Physics-Informed Neural Networks for Excited Liquid Sloshing with Beating Response in a Rectangular Tank

1. Introduction

Liquid sloshing in partially filled containers arises in naval architecture, aerospace engineering, and civil engineering, where it affects the structural integrity of LNG carriers, spacecraft fuel tanks, and seismically loaded storage vessels [1,2]. Conventional computational approaches—finite-element methods [3], finite-difference methods [4], and boundary element methods—require mesh generation and time integration; for long transient simulations spanning tens of natural periods, these costs become considerable. Reduced-order models alleviate part of the burden but are typically restricted to a pre-selected set of modes.

Physics-informed neural networks (PINNs), introduced by Raissi et al. [5], embed the governing PDEs directly into the training loss, so the network is trained without labeled solution data. Reviews by Cai et al. [6] and Cuomo et al. [7] cover applications from Navier–Stokes flows to heat transfer; the open-source library DeepXDE [8] has further broadened accessibility.

Training PINNs for oscillatory problems is difficult in practice. Wang et al. [9] showed that imbalanced multi-component losses cause gradient pathologies in which the PDE residual dominates and suppresses boundary and initial condition learning; they proposed gradient-based adaptive weighting as a remedy. Jagtap et al. [10] addressed slow convergence through adaptive activation functions. A more fundamental obstacle is spectral bias—the preference of standard fully connected networks for low-frequency components of the target function, first analyzed by Rahaman et al. [11]. For oscillatory PDEs covering many cycles this can lead to severe temporal under-resolution. Tancik et al. [12] showed that mapping inputs through sinusoidal Fourier features enables networks to learn high-frequency functions; Wang et al. [13] extended the analysis to PINNs via Neural Tangent Kernel theory. Wang et al. [14] further proposed PirateNets with causal training to enforce temporal causality and prevent PINNs from using future information to predict past states, a common pathology in time-dependent problems.

For free surface and water wave problems, several strands of work have appeared in recent years. Sheikholeslami et al. [15] examined soft and hard boundary enforcement for linear Laplace-based waves, achieving velocity errors below 3% with a periodic hard boundary layer; their trial-function approach for the bottom boundary condition, however, restricted the solution space and increased errors by up to 15 times. Li et al. [16] treated fully nonlinear wave propagation using a quasi-$\sigma$ coordinate transformation that maps the time-dependent free surface to a fixed domain, with a two-stage Adam/L-BFGS optimizer. Chen et al. [17] reconstructed rotational flow beneath periodic waves from sparse gauge data using WaveNets. More recently, Ehlers et al. [18] applied PINNs to the fully nonlinear potential flow equations for phase-resolved ocean wave assimilation and prediction, inferring the velocity potential throughout the fluid volume from surface elevation measurements alone.

Closer to the sloshing problem, Wessels et al. [19] proposed the Neural Particle Method, a Lagrangian PINN that tracks particle paths to handle the moving free surface; they simulated both small- and large-amplitude sloshing in a container and the dam-break problem, with good agreement with experimental data. Huang et al. [20] developed a boundary-and-initial-conditions-free PINN (bif-PINN) in Lagrangian coordinates that hardwires all boundary and initial conditions into the formulation via nonlinear residual connections; their approach solved the small- and large-amplitude sloshing cases accurately, whereas the original PINN failed entirely for all free-surface configurations. Shao et al. [21] combined PINNs with an improved neural particle method for violent sloshing simulation, validated against SPH results.

Despite this progress, all the above free-surface PINN studies address either periodic propagating waves over a few wavelengths or free sloshing (unforced) over a handful of cycles. The forced sloshing problem—where an external excitation near resonance produces a multi-frequency beating response spanning tens of natural periods—poses additional challenges that have not been addressed: the simulation must cover an order of magnitude more periods than a single natural cycle; the weak forcing creates free surface residuals orders of magnitude smaller than the Laplace residual, opening a zero-solution trap; and the multi-frequency temporal content exceeds the effective bandwidth of a standard MLP. The present work tackles these challenges for laterally forced sloshing at ${\omega }_e=0.9{\omega }_1$, where beating produces a response envelope over approximately $10T_1$. Four configurations (C–F) progressively address the above difficulties, reducing ${\epsilon }_u$ from 12.46% to 0.84%. The methodology is then extended to a three-dimensional rectangular tank with bi-directional excitation, demonstrating that the same techniques—Fourier embedding, hard wall boundary conditions, and gradient-enhanced regularisation—transfer effectively to the higher-dimensional setting.

2. Mathematical Formulation

2.1. Problem Setup

Consider a two-dimensional rectangular tank of length L and still-water depth h. The fluid is inviscid, incompressible, and irrotational, so the flow is described by a velocity potential $\phi (x,z,t)$. The parameters are:

$$L=2.0\mathrm{m},h=1.0\mathrm{m},n=1(\mathrm{first}\mathrm{sloshing}\mathrm{mode}).$$

(1)

The wave number and natural frequency are:

$$k_1={\displaystyle \frac{\pi }{L}}=1.5708\mathrm{rad}/\mathrm{m},$$

(2)

$${\omega }_1=\sqrt{g{k}_{1}tanh\left(k_{1}h\right)}=3.7594\mathrm{rad}/\mathrm{s},$$

(3)

with $g=9.81$ m/s². The natural period is $T_1=2\pi /{\omega }_1=1.6713$ s.

The tank is subjected to a lateral harmonic acceleration:

$$a\left(t\right)=a_{0}sin\left({\omega }_{e}t\right),a_0=0.002{\mathrm{m}/\mathrm{s}}^2,{\omega }_e=0.9{\omega }_1=3.3834\mathrm{rad}/\mathrm{s}.$$

(4)

The excitation frequency is chosen near the first natural frequency to induce a beating response with period:

$$T_{\mathrm{beat}}={\displaystyle \frac{2\pi }{|{\omega }_1-{\omega }_e|}}=16.71\mathrm{s}\approx 10T_1.$$

(5)

The simulation domain spans one full beating cycle: $t\in [0,T_{\mathrm{beat}}]$.

2.2. Governing Equation

In the non-inertial reference frame attached to the tank, the Euler equations for an inviscid, incompressible fluid read:

$${\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+(\mathbf{v}·\nabla )\mathbf{v}=-{\displaystyle \frac{1}{\rho }}\nabla p-g\widehat{\mathbf{z}}-a\left(t\right)\widehat{\mathbf{x}},$$

(6)

where $\mathbf{v}=(u,w)$ is the velocity vector, p is the pressure, $\rho$ is the (constant) fluid density, g is gravitational acceleration, and $a\left(t\right)\widehat{\mathbf{x}}$ is the pseudo-force due to the non-inertial frame. For irrotational flow ($\nabla \times \mathbf{v}=\mathbf{0}$), a velocity potential $\phi$ exists such that $\mathbf{v}=\nabla \phi$. The incompressibility constraint $\nabla ·\mathbf{v}=0$ then yields the Laplace equation:

$${\nabla }^2\phi ={\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\phi }{\partial z^2}}=0,0<x<L,-h<z<0.$$

(7)

Substituting $\mathbf{v}=\nabla \phi$ into the Euler equation (6) and using the vector identity $(\mathbf{v}·\nabla )\mathbf{v}=\nabla \left(\right|\mathbf{v}{|}^2/2)$ (valid for irrotational flow), one obtains the unsteady Bernoulli equation:

$${\displaystyle \frac{\partial \phi }{\partial t}}+{\displaystyle \frac{1}{2}}{|\nabla \phi |}^2+{\displaystyle \frac{p}{\rho }}+gz+a\left(t\right)x=C\left(t\right),$$

(8)

where $C\left(t\right)$ is a spatially uniform function absorbed into $\phi$ by redefinition. This Bernoulli equation is the starting point for deriving the dynamic free surface boundary condition.

The velocity field is recovered from $\mathbf{v}=\nabla \phi =(u,w)$.

2.3. Boundary Conditions

KBBC (impermeable bottom):

$${\displaystyle \frac{\partial \phi }{\partial z}}=0\mathrm{at}z=-h.$$

(9)

KWBC (impermeable walls):

$${\displaystyle \frac{\partial \phi }{\partial x}}=0\mathrm{at}x=0\mathrm{and}x=L.$$

(10)

KFSBC (kinematic free surface):   The free surface is defined implicitly by $F(x,z,t)=z-\eta (x,t)=0$. The kinematic condition requires that the free surface is a material surface, i.e., $\mathrm{D}F/\mathrm{D}t=0$:

$${\displaystyle \frac{\mathrm{D}F}{\mathrm{D}t}}={\displaystyle \frac{\partial F}{\partial t}}+u{\displaystyle \frac{\partial F}{\partial x}}+w{\displaystyle \frac{\partial F}{\partial z}}=-{\displaystyle \frac{\partial \eta }{\partial t}}-{\displaystyle \frac{\partial \phi }{\partial x}}{\displaystyle \frac{\partial \eta }{\partial x}}+{\displaystyle \frac{\partial \phi }{\partial z}}=0\mathrm{at}z=\eta (x,t),$$

(11)

which yields the fully nonlinear KFSBC: $\partial \eta /\partial t+(\partial \phi /\partial x)(\partial \eta /\partial x)=\partial \phi /\partial z$ at $z=\eta$. Under the small-amplitude assumption ($\eta /h=\mathcal{O}\left(ϵ\right)$, $ϵ\ll 1$), we perform a Taylor expansion about the mean free surface $z=0$:

$${\displaystyle \frac{\partial \phi }{\partial z}}{|}_{z=\eta }={\displaystyle \frac{\partial \phi }{\partial z}}{|}_{z=0}+\eta {\displaystyle \frac{{\partial }^2\phi }{\partial z^2}}{|}_{z=0}+\mathcal{O}\left({ϵ}^2\right).$$

(12)

Retaining only $\mathcal{O}\left(ϵ\right)$ terms and noting that $(\partial \phi /\partial x)(\partial \eta /\partial x)=\mathcal{O}\left({ϵ}^2\right)$, the linearized KFSBC becomes:

$${\displaystyle \frac{\partial \phi }{\partial z}}={\displaystyle \frac{\partial \eta }{\partial t}}\mathrm{at}z=0.$$

(13)

DFSBC (dynamic free surface):   Setting $p=p_{\mathrm{atm}}$ (atmospheric pressure, taken as the reference) at $z=\eta$ in the Bernoulli equation (8):

$${\displaystyle \frac{\partial \phi }{\partial t}}{|}_{z=\eta }+{\displaystyle \frac{1}{2}}{\left[{\left({\displaystyle \frac{\partial \phi }{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial \phi }{\partial z}}\right)}^2\right]}_{z=\eta }+g\eta +a\left(t\right)x=0.$$

(14)

Again applying the Taylor expansion about $z=0$:

$${\displaystyle \frac{\partial \phi }{\partial t}}{|}_{z=\eta }={\displaystyle \frac{\partial \phi }{\partial t}}{|}_{z=0}+\eta {\displaystyle \frac{{\partial }^2\phi }{\partial t\partial z}}{|}_{z=0}+\mathcal{O}\left({ϵ}^2\right),$$

(15)

and noting that ${|\nabla \phi |}^2=\mathcal{O}\left({ϵ}^2\right)$, the linearized DFSBC reads:

$${\displaystyle \frac{\partial \phi }{\partial t}}+g\eta +a\left(t\right)x=0\mathrm{at}z=0,$$

(16)

where $a\left(t\right)x$ is the pseudo-pressure term arising from the non-inertial reference frame.

Combined free surface condition.   Eliminating $\eta$ from the linearized KFSBC (13) and DFSBC (16) by differentiating the DFSBC with respect to t and substituting $\partial \eta /\partial t$ from the KFSBC yields:

$${\displaystyle \frac{{\partial }^2\phi }{\partial t^2}}+g{\displaystyle \frac{\partial \phi }{\partial z}}+\dot{a}\left(t\right)x=0\mathrm{at}z=0,$$

(17)

where $\dot{a}\left(t\right)=a_0{\omega }_{e}cos\left({\omega }_{e}t\right)$. This combined condition, while not used directly in the PINN loss, provides physical insight: the free surface acts as a dispersive wave guide with the forcing $\dot{a}\left(t\right)x$ injecting energy into the system.

2.4. Initial Condition

Since the tank starts from rest, $\phi (x,z,0)=0$ and $\partial \phi /\partial t(x,z,0)=0$. No terminal condition is imposed because $\phi$ does not vanish at $t=T_{\mathrm{beat}}$; the two initial conditions together serve as temporal anchors.

3. PINN Methodology

3.1. Network Architecture

The neural network approximation ${\widehat{\phi }}_{\theta }(x,z,t)$ is a fully connected feed-forward network (multilayer perceptron) defined by the recursive composition:

$${\mathbf{h}}^{\left(0\right)}={\mathbf{x}}_{\mathrm{input}},{\mathbf{h}}^{(\ell )}=tanh\left({\mathbf{W}}^{(\ell )}{\mathbf{h}}^{(\ell -1)}+{\mathbf{b}}^{(\ell )}\right),\widehat{\phi }={\mathbf{W}}^{(L+1)}{\mathbf{h}}^{\left(L\right)}+b^{(L+1)},$$

(18)

where ${\mathbf{W}}^{(\ell )}\in {\mathbb{R}}^{n_{\ell }\times n_{\ell -1}}$ and ${\mathbf{b}}^{(\ell )}\in {\mathbb{R}}^{n_{\ell }}$ are the weight matrix and bias vector for layer ℓ. The collection of all trainable parameters is $\theta ={\{{\mathbf{W}}^{(\ell )},{\mathbf{b}}^{(\ell )}\}}_{\ell =1}^{L+1}$. The network uses $L=16$ hidden layers with $n_{\ell }=64$ neurons, Xavier normal initialization, and a single scalar output. Velocities follow from automatic differentiation:

$$u={\displaystyle \frac{\partial \widehat{\phi }}{\partial x}},w={\displaystyle \frac{\partial \widehat{\phi }}{\partial z}}.$$

(19)

In the enhanced configuration (Config E), the input is augmented with Fourier time features, increasing the input dimension from 3 to 9 (see Section 3.6). The configuration is summarized in Table 1.

3.2. Loss Function and Automatic Differentiation

All required derivatives of the network output $\widehat{\phi }=f_{\theta }\left(\gamma \left(\mathbf{x}\right)\right)$ are computed via automatic differentiation (AD). By the chain rule through the Fourier embedding, the physical derivatives involve contributions from both the normalized coordinates and the sinusoidal features. For instance:

$${\displaystyle \frac{\partial \widehat{\phi }}{\partial t}}={\displaystyle \frac{\partial f_{\theta }}{\partial \tilde{t}}}·{\displaystyle \frac{1}{T_{\mathrm{sim}}}}+\sum _{i=1}^3\left({\displaystyle \frac{\partial f_{\theta }}{\partial s_i}}{\omega }_{i}cos{\omega }_{i}t-{\displaystyle \frac{\partial f_{\theta }}{\partial c_i}}{\omega }_{i}sin{\omega }_{i}t\right),$$

(20)

where $s_i=sin{\omega }_{i}t$, $c_i=cos{\omega }_{i}t$, and ${\omega }_i\in \{{\omega }_1,{\omega }_e,\Delta \omega \}$. The spatial derivatives are simpler: $\partial \widehat{\phi }/\partial x=(\partial f_{\theta }/\partial \tilde{x})/L$ and $\partial \widehat{\phi }/\partial z=(\partial f_{\theta }/\partial \tilde{z})/h$. The second-order derivatives needed for the Laplace residual are:

$${\displaystyle \frac{{\partial }^2\widehat{\phi }}{\partial x^2}}={\displaystyle \frac{1}{L^2}}{\displaystyle \frac{{\partial }^2{f}_{\theta }}{\partial {\tilde{x}}^2}},{\displaystyle \frac{{\partial }^2\widehat{\phi }}{\partial z^2}}={\displaystyle \frac{1}{h^2}}{\displaystyle \frac{{\partial }^2{f}_{\theta }}{\partial {\tilde{z}}^2}}.$$

(21)

In PyTorch, these are obtained via nested calls to torch.autograd.grad with create_graph=True.

The analytically derived $\eta (x,t)$ from Eq. (44) and its time derivative $\partial \eta /\partial t$ are directly substituted into the KFSBC and DFSBC as known boundary data:

$$r_{\mathrm{KFSBC}}\left(\mathbf{p}\right):={\displaystyle \frac{\partial \widehat{\phi }}{\partial z}}{|}_{z=0}-{\displaystyle \frac{\partial {\eta }_{\mathrm{exact}}}{\partial t}},r_{\mathrm{DFSBC}}\left(\mathbf{p}\right):={\displaystyle \frac{\partial \widehat{\phi }}{\partial t}}{|}_{z=0}+g{\eta }_{\mathrm{exact}}+a\left(t\right)x.$$

(22)

This converts the free surface conditions from coupled evolution equations into deterministic Neumann–Dirichlet boundary data for the Laplace equation. The remaining residual operators are:

$$r_{\mathrm{GE}}\left(\mathbf{p}\right):={\displaystyle \frac{{\partial }^2\widehat{\phi }}{\partial x^2}}+{\displaystyle \frac{{\partial }^2\widehat{\phi }}{\partial z^2}},r_{\mathrm{KBBC}}\left(\mathbf{p}\right):={\displaystyle \frac{\partial \widehat{\phi }}{\partial z}}{|}_{z=-h}.$$

(23)

The lateral wall boundary condition $\partial \widehat{\phi }/\partial x=0$ at $x\in \{0,B\}$ is enforced exactly by the network architecture rather than as a soft penalty. The spatial input is transformed as $x_n=cos(\pi x/B)$, so that the network output $\widehat{\phi }$ depends on x only through $cos(\pi x/B)$. By the chain rule:

$${\displaystyle \frac{\partial \widehat{\phi }}{\partial x}}={\displaystyle \frac{\partial f_{\theta }}{\partial x_n}}·\left(-{\displaystyle \frac{\pi }{B}}sin{\displaystyle \frac{\pi x}{B}}\right),$$

(24)

which vanishes identically at $x=0$ and $x=B$. This hard boundary condition removes the need for wall collocation points and reduces the number of loss terms.

A gradient-enhanced Laplace regulariser is added to enforce velocity smoothness. Since ${\nabla }^2\phi =0$ implies ${\nabla }^{2}u=0$ and ${\nabla }^{2}w=0$ for harmonic functions, we penalise the spatial gradients of the Laplace residual:

$$r_{\mathrm{GL}}\left(\mathbf{p}\right):={∥\nabla \left({\nabla }^2\widehat{\phi }\right)∥}^2={\left({\displaystyle \frac{\partial }{\partial x}}{\nabla }^2\widehat{\phi }\right)}^2+{\left({\displaystyle \frac{\partial }{\partial z}}{\nabla }^2\widehat{\phi }\right)}^2.$$

(25)

This requires third-order automatic differentiation and is evaluated on a randomly subsampled subset of 4,096 interior points per step to limit GPU memory consumption.

The remaining residuals are:

$$r_{\mathrm{IC}}\left(\mathbf{p}\right):=\widehat{\phi }(x,z,0),r_{\mathrm{IC}}\mathrm{v}\left(\mathbf{p}\right):={\displaystyle \frac{\partial \widehat{\phi }}{\partial t}}(x,z,0).$$

(26)

Each individual loss component is the mean squared residual over the corresponding point set $S_m$:

$${\mathcal{J}}_m={\displaystyle \frac{1}{|S_m|}}\sum _{\mathbf{p}\in S_m}{\left[r_m\left(\mathbf{p}\right)\right]}^2.$$

(27)

Residual normalization. The main difficulty in the forced sloshing problem is the disparity in residual magnitudes. The GE residual involves second-order spatial derivatives ($\mathcal{O}\left(1\right)$), while the DFSBC contains $\partial \widehat{\phi }/\partial t\sim \Phi {\omega }_1$ and the KFSBC contains $\partial \widehat{\phi }{/\partial z|}_{z=0}\sim \Phi k_1$, where $\Phi$ is the characteristic potential amplitude. To equalize their contributions, we define the characteristic scale:

$$\Phi =\underset{x,t}{max}\left|{\phi }_{\mathrm{exact}}(x,0,t)\right|,$$

(28)

and normalize the free surface residuals:

$${\tilde{r}}_{\mathrm{DFSBC}}\left(\mathbf{p}\right)={\displaystyle \frac{r_{\mathrm{DFSBC}}\left(\mathbf{p}\right)}{\Phi {\omega }_1}},{\tilde{r}}_{\mathrm{KFSBC}}\left(\mathbf{p}\right)={\displaystyle \frac{r_{\mathrm{KFSBC}}\left(\mathbf{p}\right)}{\Phi k_1}}.$$

(29)

The normalized MSE losses are then:

$${\tilde{\mathcal{J}}}_{\mathrm{DFSBC}}={\displaystyle \frac{1}{|S_{\mathrm{FS}}|}}\sum _{\mathbf{p}\in S_{\mathrm{FS}}}{\left[{\tilde{r}}_{\mathrm{DFSBC}}\left(\mathbf{p}\right)\right]}^2,{\tilde{\mathcal{J}}}_{\mathrm{KFSBC}}={\displaystyle \frac{1}{|S_{\mathrm{FS}}|}}\sum _{\mathbf{p}\in S_{\mathrm{FS}}}{\left[{\tilde{r}}_{\mathrm{KFSBC}}\left(\mathbf{p}\right)\right]}^2.$$

(30)

The total weighted loss is:

$${\mathcal{J}}_{\mathrm{total}}=\underset{\mathrm{Laplace}}{\underbrace{{\mathcal{J}}_{\mathrm{GE}}}}+\underset{\mathrm{grad}-\mathrm{enhanced}}{\underbrace{{\lambda }_{\mathrm{GL}}{\mathcal{J}}_{\mathrm{GL}}}}+\underset{\mathrm{bottom}}{\underbrace{{\mathcal{J}}_{\mathrm{KBBC}}}}+\underset{\mathrm{free}}{\underbrace{{\lambda }_D{\tilde{\mathcal{J}}}_{\mathrm{DFSBC}}+{\lambda }_K{\tilde{\mathcal{J}}}_{\mathrm{KFSBC}}}}\mathrm{surface}\left(\mathrm{weighted}\right)+\underset{\mathrm{initial}}{\underbrace{{\lambda }_I\left({\mathcal{J}}_{\mathrm{IC}}+{\mathcal{J}}_{\mathrm{IC}}\mathrm{v}\right)}}\mathrm{conditions},$$

(31)

where ${\lambda }_D={\lambda }_K=100$, ${\lambda }_{\mathrm{GL}}=0.1$, and ${\lambda }_I=10$. The lateral wall condition is enforced exactly by the $cos(\pi x/B)$ embedding and does not appear as a loss term.

3.3. Collocation Points

Points are sampled on a structured grid of $N_x\times N_z\times N_t$ interior points in $(x,z,t)\in [0,B]\times [-h,0]\times [0,T_{\mathrm{beat}}]$, with the vertical coordinate clustered toward the free surface via $z=-h(1-{\xi }^{0.6})$ where $\xi \in [0,1]$ is uniform. An additional 20 dense temporal points near $t=0$ resolve the initial transient. Boundary points are sampled on the bottom ($N_x\times N_t$) and free surface ($N_x\times N_t$). No wall collocation points are needed due to the hard wall embedding. The initial condition is enforced on $N_x\times N_z$ points at $t=0$. The default grid uses $N_x=31$, $N_z=31$, $N_t=161$, yielding ∼173,000 interior points per window.

3.4. Training Strategy

Training proceeds in two stages with stochastic mini-batch sampling to manage GPU memory. At each Adam step, 15% of the interior, free-surface, and bottom points are randomly subsampled and transferred to GPU (all initial-condition points are retained). For the L-BFGS phase, a fixed 30% subsample is used to maintain the deterministic closure required by the quasi-Newton algorithm. This keeps peak GPU memory below 8 GB regardless of the total grid size.

1.

Adam (first-order): For the single-beat case, 8,000 epochs are used; for the three-beat case, 20,000 epochs with a 500-step linear warmup followed by cosine annealing:

$${\alpha }_k={\alpha }_{\min }+{\displaystyle \frac{1}{2}}({\alpha }_{\max }-{\alpha }_{\min })\left(1+cos{\displaystyle \frac{(k-k_{\mathrm{warm}})\pi }{K_{\mathrm{Adam}}-k_{\mathrm{warm}}}}\right),$$

(32)

where ${\alpha }_{\max }={10}^{-3}$, ${\alpha }_{\min }={10}^{-5}$, and $k_{\mathrm{warm}}=500$. Gradient clipping ($\parallel {\nabla }_{\theta }\mathcal{J}\parallel \le 1$) is applied throughout.

2.

L-BFGS (quasi-Newton): Up to 2,000 iterations (single-beat) or 3,000 iterations (three-beat) with strong Wolfe line search ($c_1={10}^{-4}$, $c_2=0.9$), gradient tolerance ${10}^{-12}$, function change tolerance ${10}^{-14}$, and history size $m=100$. For the three-beat mini-batch case, the batch is rotated every five outer L-BFGS steps to improve domain coverage while maintaining the deterministic closure required by the algorithm. The Hessian approximation follows the standard BFGS update:

$$H_k\approx B_k^{-1},B_{k+1}=B_k-{\displaystyle \frac{B_k{s}_k{s}_k^{\top }{B}_k}{s_k^{\top }{B}_k{s}_k}}+{\displaystyle \frac{y_k{y}_k^{\top }}{y_k^{\top }{s}_k}},$$

(33)

where $s_k={\theta }_{k+1}-{\theta }_k$ and $y_k={\nabla }_{\theta }{\mathcal{J}}_{k+1}-{\nabla }_{\theta }{\mathcal{J}}_k$.

3.5. Error Metric

The accuracy is quantified by the normalized mean absolute error (NMAE):

$${\epsilon }_q[\%]={\displaystyle \frac{{\displaystyle \sum _{i=1}^{\left|E\right|}|q_{\mathrm{PINN}}\left({\mathbf{p}}_i\right)-q_{\mathrm{exact}}\left({\mathbf{p}}_i\right)|}}{{\displaystyle \left|E\right|\underset{\mathbf{p}\in E}{max}\left|q_{\mathrm{exact}}\left(\mathbf{p}\right)\right|}}}\times 100,$$

(34)

where $q\in \{u,w,\eta \}$ is any field quantity, and $E=\left\{(x_i,z_j,t_k)\right\}$ is a separate evaluation grid ($N_x^E\times N_z^E\times N_t^E$) that does not coincide with the training collocation points. We also report the relative $L^2$ error for comparison:

$${\epsilon }_q^{L^2}[\%]={\displaystyle \frac{{\left({\displaystyle \sum _{i=1}^{\left|E\right|}{[q_{\mathrm{PINN}}\left({\mathbf{p}}_i\right)-q_{\mathrm{exact}}\left({\mathbf{p}}_i\right)]}^2}\right)}^{1/2}}{{\left({\displaystyle \sum _{i=1}^{\left|E\right|}{\left[q_{\mathrm{exact}}\left({\mathbf{p}}_i\right)\right]}^2}\right)}^{1/2}}}\times 100.$$

(35)

3.6. Fourier Time Embedding

Standard MLPs with tanh activations exhibit spectral bias [11]: the normalized input $t/T_{\mathrm{sim}}\in [0,1]$ must encode oscillations at ${\omega }_1{T}_{\mathrm{sim}}\approx 63$ rad, far beyond the effective bandwidth of a standard MLP. Following Tancik et al. [12] and Wang et al. [13], we augment the input with sinusoidal functions at the three physical frequencies of the problem:

$$\gamma \left(\mathbf{x}\right)=\left(\tilde{x},\tilde{z},\tilde{t},sin\left(\omega t\right),cos\left(\omega t\right)\right)\in {\mathbb{R}}^9,$$

(36)

where $\tilde{x}=x/L$, $\tilde{z}=(z+h)/h$, $\tilde{t}=t/T_{\mathrm{sim}}$ are normalised coordinates and $\omega ={({\omega }_1,{\omega }_e,\Delta \omega )}^{\top }$ with $\Delta \omega ={\omega }_1-{\omega }_e$. Explicitly:

$${\mathbf{x}}_{\mathrm{input}}=\gamma \left(\mathbf{x}\right)=\left({\displaystyle \frac{x}{L}},{\displaystyle \frac{z+h}{h}},{\displaystyle \frac{t}{T_{\mathrm{sim}}}},sin{\omega }_{1}t,cos{\omega }_{1}t,sin{\omega }_{e}t,cos{\omega }_{e}t,sin\Delta \omega t,cos\Delta \omega t\right).$$

(37)

This can be written compactly as $\gamma \left(\mathbf{x}\right)={({\tilde{\mathbf{x}}}^{\top },sin{\left(\mathbf{B}t\right)}^{\top },cos{\left(\mathbf{B}t\right)}^{\top })}^{\top }$ where $\mathbf{B}={({\omega }_1,{\omega }_e,\Delta \omega )}^{\top }$.

Unlike the random Fourier features of [13], the frequencies here are known a priori from the problem physics. The embedding increases the input dimension from 3 to 9, adding only $6\times 64=384$ parameters ($<1\%$ of total).

3.7. PINN Training Algorithm

The complete training procedure is summarized in Algorithm 1.

Algorithm 1:PINN Training for Excited Sloshing with Fourier Embedding

Require:  Tank parameters ($L,h,a_0,{\omega }_e$), network config ($N_L,N_n$), grid ($N_x,N_z,N_t$) Ensure:  Trained network ${\widehat{\phi }}_{\theta }(x,z,t)$ 1: Precompute:${\omega }_1,k_n,A_n\left(t\right),{\eta }_{\mathrm{exact}}(x,t),\partial \eta /\partial t$ 2: Generate collocation points: interior, boundary, IC on structured grid 3: Initialize: MLP with Xavier normal, ${\theta }_0$ 4: Construct Fourier input: Augment $(x,z,t)$ with $sin\left({\omega }_{1}t\right),cos\left({\omega }_{1}t\right),\dots$$\to {\mathbb{R}}^9$ 5: for epoch $=1$ to 8000 (Adam) do 6:    Forward pass: $\widehat{\phi }={\mathrm{MLP}}_{\theta }\left({\mathbf{x}}_{\mathrm{input}}\right)$ 7:    Compute residuals via auto-diff: $r_{\mathrm{GE}},r_{\mathrm{KBBC}},r_{\mathrm{KWBC}},r_{\mathrm{KFSBC}},r_{\mathrm{DFSBC}},r_{\mathrm{IC}},r_{\mathrm{IC}}\mathrm{v}$ 8:    Normalize: ${\tilde{r}}_{\mathrm{DFSBC}}\leftarrow r_{\mathrm{DFSBC}}/\left(\Phi {\omega }_1\right)$,    ${\tilde{r}}_{\mathrm{KFSBC}}\leftarrow r_{\mathrm{KFSBC}}/\left(\Phi k_1\right)$ 9:    Compute weighted total loss: $\mathcal{J}={\mathcal{J}}_{\mathrm{GE}}+\dots +{\lambda }_D{\tilde{\mathcal{J}}}_{\mathrm{DFSBC}}+{\lambda }_K{\tilde{\mathcal{J}}}_{\mathrm{KFSBC}}+{\lambda }_I({\mathcal{J}}_{\mathrm{IC}}+{\mathcal{J}}_{\mathrm{IC}}\mathrm{v})$ 10:    Clip gradients: $\parallel {\nabla }_{\theta }\mathcal{J}\parallel \le 1.0$ 11:    Update $\theta$ with Adam (cosine LR: ${10}^{-3}\to {10}^{-5}$) 12: end for 13: for iter $=1$ to 2000 (L-BFGS) do 14:    Compute $\mathcal{J}$ and ${\nabla }_{\theta }\mathcal{J}$ (full batch) 15:    Update $\theta$ with L-BFGS (strong Wolfe, history=100) 16: end for 17: Evaluate: Compute $u,w,\eta$ on separate evaluation grid 18: Return:${\widehat{\phi }}_{\theta }$, error metrics ${\epsilon }_u,{\epsilon }_w,{\epsilon }_{\eta }$

3.8. PINN Framework Overview

Figure 1 shows the overall architecture. Coordinates $(x,z,t)$ are first mapped through the Fourier embedding into a 9-dimensional input; this feeds a 16-layer MLP with 64 neurons per layer and tanh activations. All PDE and boundary residuals are formed by automatic differentiation of the scalar output $\widehat{\phi }$, assembled into the weighted loss, and minimised by backpropagation.

4. Results and Discussion

4.1. Analytical Solution

PINN predictions are validated against the closed-form solution of the linearised problem. The wall and Laplace equations admit spatial eigenfunctions ${\psi }_n(x,z)=cos\left(k_{n}x\right)cosh\left[k_n(h+z)\right]/cosh\left(k_{n}h\right)$ with $k_n=n\pi /L$ for odd n. Substituting the modal ansatz

$$\phi =\sum _{n\mathrm{odd}}{A}_n\left(t\right){\psi }_n(x,z)$$

(38)

into the combined free surface condition (17) and projecting onto each mode yields the forced harmonic oscillator ODE for the modal amplitude:

$$A_n^{\prime \prime }+{\omega }_n^2{A}_n=-\dot{a}\left(t\right)c_n,\dot{a}\left(t\right)=a_0{\omega }_{e}cos\left({\omega }_{e}t\right),$$

(39)

where the natural frequency, forcing coefficient, and dispersion relation are:

$${\omega }_n=\sqrt{g{k}_{n}tanh\left(k_{n}h\right)},c_n={\displaystyle \frac{-4L}{{\left(n\pi \right)}^2}},k_n={\displaystyle \frac{n\pi }{L}}.$$

(40)

With zero initial conditions $A_n\left(0\right)=0$, $A_n^{\prime }\left(0\right)=0$, the particular solution is obtained via the method of undetermined coefficients:

$$A_n\left(t\right)=C_n\left[cos\left({\omega }_{e}t\right)-cos\left({\omega }_{n}t\right)\right],C_n={\displaystyle \frac{-a_0{\omega }_e{c}_n}{{\omega }_n^2-{\omega }_e^2}}.$$

(41)

The general solution of Eq. (39) is $A_n\left(t\right)={\alpha }_{n}cos\left({\omega }_{n}t\right)+{\beta }_{n}sin\left({\omega }_{n}t\right)+A_n^{\left(p\right)}\left(t\right)$, where the particular solution $A_n^{\left(p\right)}\left(t\right)=C_{n}cos\left({\omega }_{e}t\right)$. Applying the initial conditions determines ${\alpha }_n=-C_n$, ${\beta }_n=0$, recovering Eq. (41).

For the first mode ($n=1$), the beating structure is explicit via the product-to-sum identity:

$$cos\left({\omega }_{e}t\right)-cos\left({\omega }_{1}t\right)=2sin\left[{\displaystyle \frac{({\omega }_1-{\omega }_e)t}{2}}\right]sin\left[{\displaystyle \frac{({\omega }_1+{\omega }_e)t}{2}}\right],$$

(42)

exhibiting a slow envelope at the beat frequency $\Delta \omega ={\omega }_1-{\omega }_e$ and a fast carrier at the mean frequency $\overline{\omega }=({\omega }_1+{\omega }_e)/2$.

The complete velocity potential and free surface elevation are:

$$\phi (x,z,t)={\phi }_0\left(t\right)+\sum _{n\mathrm{odd}}{A}_n\left(t\right)cos\left(k_{n}x\right){\displaystyle \frac{cosh\left[k_n(h+z)\right]}{cosh\left(k_{n}h\right)}},$$

(43)

$$\eta (x,t)=-{\displaystyle \frac{1}{g}}\left[\sum _{n\mathrm{odd}}{A}_n^{\prime }\left(t\right)cos\left(k_{n}x\right)+a\left(t\right)\left(x-{\displaystyle \frac{L}{2}}\right)\right],$$

(44)

where ${\phi }_0\left(t\right)=a_{0}L/\left(2{\omega }_e\right)[cos\left({\omega }_{e}t\right)-1]$ ensures volume conservation (${\int }_0^L\eta \mathrm{d}x=0$), and the first five odd modes ($n=1,3,5,7,9$) are retained.

The analytical velocity components are obtained by differentiating Eq. (43):

$$u(x,z,t)={\displaystyle \frac{\partial \phi }{\partial x}}=-\sum _{n\mathrm{odd}}{k}_n{A}_n\left(t\right)sin\left(k_{n}x\right){\displaystyle \frac{cosh\left[k_n(h+z)\right]}{cosh\left(k_{n}h\right)}},$$

(45)

$$w(x,z,t)={\displaystyle \frac{\partial \phi }{\partial z}}=\sum _{n\mathrm{odd}}{k}_n{A}_n\left(t\right)cos\left(k_{n}x\right){\displaystyle \frac{sinh\left[k_n(h+z)\right]}{cosh\left(k_{n}h\right)}}.$$

(46)

The time derivative of the modal amplitude, needed for $\eta$, is:

$$A_n^{\prime }\left(t\right)=C_n\left[-{\omega }_{e}sin\left({\omega }_{e}t\right)+{\omega }_{n}sin\left({\omega }_{n}t\right)\right].$$

(47)

The wave elevation is recovered from the DFSBC:

$${\eta }_{\mathrm{PINN}}=-{\displaystyle \frac{1}{g}}\left[{\displaystyle \frac{\partial \widehat{\phi }}{\partial t}}{|}_{z=0}+a\left(t\right)x\right].$$

(48)

4.2. Loss Landscape Analysis

This subsection explains why standard PINN training fails for forced sloshing and how the proposed remedies work.

The zero-solution trap. The trivial solution $\widehat{\phi }=0$ satisfies the Laplace equation, all wall/bottom boundary conditions, and both initial conditions exactly. The only nonzero residuals are on the free surface:

$$r_{\mathrm{DFSBC}}{|}_{\widehat{\phi }=0}=g{\eta }_{\mathrm{exact}}+a\left(t\right)x=-{\displaystyle \frac{\partial {\phi }_{\mathrm{exact}}}{\partial t}}{|}_{z=0},$$

(49)

$$r_{\mathrm{KFSBC}}{|}_{\widehat{\phi }=0}=-{\displaystyle \frac{\partial {\eta }_{\mathrm{exact}}}{\partial t}}.$$

(50)

Since the first mode dominates ($|C_3/C_1|<1\%$), the DFSBC loss at $\widehat{\phi }=0$ can be evaluated analytically retaining only $n=1$:

$${\mathcal{J}}_{\mathrm{DFSBC}}{|}_{\widehat{\phi }=0}\approx {\displaystyle \frac{C_1^2({\omega }_e^2+{\omega }_1^2)}{4}}+{\displaystyle \frac{a_0^2{L}^2}{8}}=2.87\times {10}^{-5},$$

(51)

using $\langle {sin}^2\omega t\rangle \approx 1/2$ and $\langle {cos}^2{k}_{1}x\rangle =1/2$ (verified against numerical quadrature to 0.2% accuracy). A randomly initialized network produces ${\mathcal{J}}_{\mathrm{GE}}^{(}\mathrm{init})\sim \mathcal{O}\left({10}^{-1}\right)$, yielding the ratio

$${\displaystyle \frac{{\mathcal{J}}_{\mathrm{DFSBC}}{|}_{\widehat{\phi }=0}}{{\mathcal{J}}_{\mathrm{GE}}^{(}\mathrm{init})}}\sim {10}^{-4},$$

(52)

so the DFSBC penalty for the zero solution is negligible—the optimizer has no incentive to deviate from $\widehat{\phi }=0$.

Breaking the trap. Two mechanisms are combined to escape the zero-solution trap. First, the DFSBC and KFSBC residuals are divided by their characteristic scales $\Phi {\omega }_1$ and $\Phi k_1$, which amplifies the normalized zero-solution loss from $2.87\times {10}^{-5}$ to ${\tilde{\mathcal{J}}}_{\mathrm{DFSBC}}{|}_{\widehat{\phi }=0}=0.121$. Second, applying a weight ${\lambda }_D=100$ further elevates this to ${\lambda }_D{\tilde{\mathcal{J}}}_{\mathrm{DFSBC}}{|}_{\widehat{\phi }=0}=12.1$, which exceeds ${\mathcal{J}}_{\mathrm{GE}}^{(}\mathrm{init})$ by two orders of magnitude and reshapes the loss landscape so that $\widehat{\phi }=0$ is no longer a local minimum.

4.3. External Validation Against Chen & Nokes (2005)

The PINN was re-run with the parameters of Chen & Nokes [4] ($B=0.9$ m, $d_0=0.6$ m, $X_0/d_0=0.00186$, ${\omega }_e/{\omega }_1=0.9$, Fig. 5c therein). As shown in Figure 2, the PINN agrees with the linear analytical solution to within ${\epsilon }_{\eta }=0.15\%$ over the full beating cycle ($T_{\mathrm{beat}}=10T_1$). The slightly lower peak amplitude in Chen’s viscous nonlinear FDM is attributable to viscous dissipation and nonlinear effects absent from the present inviscid linearised formulation.

Ablation of remaining techniques.Table 2 isolates the effect of temporal resolution and Fourier embedding, starting from the configuration where the zero-solution trap has already been broken (normalized + weighted, ${\lambda }_D={\lambda }_K=100$).

Config C ($N_t=51\approx 5$ pts/$T_1$) gives ${\epsilon }_{\eta }=0.30\%$ but ${\epsilon }_w=21.36\%$ due to temporal under-resolution. Increasing to $N_t=201$ ($\approx 20$ pts/$T_1$, Config D) reduces ${\epsilon }_w$ to 15.40%. The Fourier embedding (Config E) alleviates spectral bias: ${\epsilon }_{\eta }$ drops to 0.18% and ${\epsilon }_w$ to 4.68%. Config F adds the hard wall boundary condition (removing wall collocation points) and the gradient-enhanced Laplace regulariser, yielding ${\epsilon }_u=0.84\%$ and ${\epsilon }_w=1.65\%$—a further $4.5\times$ improvement in velocity accuracy. Config F is used for all subsequent simulations.

4.4. Velocity Fields

Figure 3 compares the analytical and PINN velocity fields at $t=T_{\mathrm{beat}}/4$ (within Window 0 of the 3-beat decomposition), together with the pointwise error. The horizontal velocity u in Window 0 has ${\epsilon }_u=0.83\%$; the PINN reproduces the sinusoidal spatial variation and the cosh-dependent depth profile. The largest u-errors appear near the tank bottom ($z/h<-0.6$), where velocity magnitudes are small and relative errors are amplified.

The vertical velocity w in Window 0 has ${\epsilon }_w=1.50\%$, with errors concentrated near the bottom corners where the sinh depth profile attenuates the signal to near-zero values. In the free surface region ($z/h>-0.3$), both velocity components have errors below 2%. The overall 3-beat errors (${\epsilon }_u=3.26\%$, ${\epsilon }_w=3.87\%$) are higher because of the later windows (see Table 4).

5. Simulation Results

The single-beat case ($T_{\mathrm{sim}}=T_{\mathrm{beat}}=10T_1$) is examined first through wave elevation, velocity potential, and flow fields; the method is then extended to a three-beat simulation ($3T_{\mathrm{beat}}=30T_1$) to test long-time accuracy.

5.1. Wave Elevation

Figure 4 shows the wave elevation at the left wall ($x=0$) and right wall ($x=L$) over three beating cycles ($3T_{\mathrm{beat}}=30T_1$). Both walls are antinodes of the first sloshing mode, with peak amplitudes of $\pm 1.65$ mm. The PINN predictions track the analytical solution throughout all three beats (${\epsilon }_{\eta }=0.30\%$), reproducing the growth, peak, and decay phases of each beating envelope without visible drift.

5.2. Velocity Potential at Free Surface

Figure 5 shows the velocity potential $\phi$ at the free surface ($z=0$) as a function of $x/B$ and $t/T_1$ over three beating cycles. The analytical field exhibits alternating positive and negative bands whose amplitude follows the beating envelope. The PINN output is visually identical; the maximum absolute error is $9.5\times {10}^{-4}$ m²/s against a potential amplitude of $2.5\times {10}^{-2}$ m²/s. Errors are largest in Window 2, consistent with the per-window velocity errors in Table 4. No discontinuities appear at the window boundaries (dashed lines).

5.3. Frequency Parameter Study

To test how the method performs across different excitation conditions, five frequency ratios are selected: ${\omega }_e/{\omega }_1=0.5$ (far from resonance), $0.8$ (moderately close), $0.9$ (baseline), $0.95$ (near-resonance), and $1.1$ (above resonance). All simulations use the Chen & Nokes [4] tank geometry ($B=0.9$ m, $d_0=0.6$ m, $X_0/d_0=0.00186$) and the same network architecture (16 layers × 64 neurons with Fourier embedding).

Figure 6 shows the left-wall wave elevation for all five cases. The beating period $T_{\mathrm{beat}}=2\pi /|{\omega }_1-{\omega }_e|$ ranges from $2.0T_1$ (${\omega }_e/{\omega }_1=0.5$) to $20.0T_1$ (${\omega }_e/{\omega }_1=0.95$), and the peak amplitude grows sharply as ${\omega }_e/{\omega }_1\to 1$. In each case the PINN reproduces both the carrier oscillations and the slowly varying envelope.

Table 3 lists the error metrics for all frequency ratios. The wave elevation error ${\epsilon }_{\eta }$ stays below $0.25\%$ in every case. For excitation frequencies near the natural frequency (${\omega }_e/{\omega }_1=0.8$–$1.1$), velocity errors are low (${\epsilon }_u<2.3\%$, ${\epsilon }_w<2.9\%$). The far-from-resonance case (${\omega }_e/{\omega }_1=0.5$) shows large relative velocity errors because the response amplitude is an order of magnitude smaller, so the normalisation denominator shrinks; absolute velocity errors are comparable. Training time is about 1900 s per beating cycle at all frequency ratios.

Figure 7 and Figure 8 compare the PINN velocity fields at $t=T_{\mathrm{beat}}/4$ for ${\omega }_e/{\omega }_1=0.90$ (baseline) and $0.95$ (near-resonance). Each figure shows the analytical solution (top), PINN prediction (middle), and pointwise error (bottom) for both u and w. Errors concentrate near the bottom corners where velocity gradients are steepest. At ${\omega }_e/{\omega }_1=0.95$ the response amplitude is the largest among all cases, yet ${\epsilon }_u=1.05\%$—the Fourier embedding handles the beating dynamics at all tested frequency ratios.

5.4. Long-Time Multi-Beat Simulation via Time-Domain Decomposition

A single PINN trained over a long time window runs into a capacity bottleneck: ${\epsilon }_{\eta }$ stays low, but velocity errors grow with the temporal domain (e.g. ${\epsilon }_u=15.79\%$ for $3T_{\mathrm{beat}}$ with one network, versus $0.84\%$ for $1T_{\mathrm{beat}}$). We therefore partition the simulation into one-beat windows.

Method. The simulation domain $[0,N_w{T}_{\mathrm{beat}}]$ is partitioned into $N_w$ windows, each covering one beating period. For window k:

• The temporal domain is $t\in [k{T}_{\mathrm{beat}},(k+1)T_{\mathrm{beat}}]$ with local normalisation $\widehat{t}=(t-k{T}_{\mathrm{beat}})/T_{\mathrm{beat}}$.
• The initial condition at $t=k{T}_{\mathrm{beat}}$ is prescribed from the analytical solution: $\phi ={\phi }_{\mathrm{exact}}$, $\partial \phi /\partial t={(\partial \phi /\partial t)}_{\mathrm{exact}}$.
• Transfer learning: window k is initialised with the trained weights from window $k-1$, providing a warm start that reduces training time by approximately $30\%$ compared to training from scratch.

Each window is trained independently (8,000 Adam + 2,000 L-BFGS steps) with the same collocation density as the single-beat case, so the per-window difficulty is the same as the already validated single-beat problem.

Figure 9 shows the wave elevation over five consecutive beating cycles ($5T_{\mathrm{beat}}=50T_1$). The predictions from all five windows are stitched together and compared with the analytical solution; no visible degradation appears at the window boundaries.

Table 4 compares the per-window and overall errors for the 3-beat and 5-beat simulations. For the 3-beat case, the decomposition gives ${\epsilon }_u=3.26\%$ and ${\epsilon }_w=3.87\%$ overall, whereas the single-network approach yields ${\epsilon }_u=15.79\%$ for the same domain—a reduction by nearly a factor of five.

Extending to 5 beats ($50T_1$), the overall errors are ${\epsilon }_u=3.43\%$, ${\epsilon }_w=4.07\%$, and ${\epsilon }_{\eta }=0.30\%$. The per-window errors do not accumulate monotonically: Window 3 (${\epsilon }_u=2.60\%$) is more accurate than Window 2 ($6.44\%$), indicating that accuracy fluctuates between windows rather than systematically degrading. The total training time for 5 windows is 9,760 s ($\approx 2.7$ h), scaling linearly with the number of windows.

Figure 10 shows the PINN velocity field at $t=3.5T_{\mathrm{beat}}$ (middle of the fourth window in the 5-beat simulation). The spatial patterns match the single-beat results, with no sign of progressive degradation as more windows are added.

6. Extension to Three Dimensions

The preceding sections demonstrated that PINNs with Fourier embedding and hard wall boundary conditions can accurately capture multi-frequency beating in a 2-D tank. A natural question is whether these techniques extend to a three-dimensional configuration, where the additional spatial dimension introduces transverse modes and substantially increases the computational domain. This section presents the 3-D formulation, the required modifications to the network architecture, and preliminary results.

6.1. 3-D Problem Setup

Consider a three-dimensional rectangular tank of length B, width W, and still-water depth $d_0$. The physical parameters are taken from the 2-D Chen & Nokes configuration with an added transverse dimension:

$$B=0.9\mathrm{m},W=0.6\mathrm{m},d_0=0.6\mathrm{m}.$$

(53)

The first-mode wave numbers and natural frequencies in each horizontal direction are: `small

$$k_{1x}={\displaystyle \frac{\pi }{B}},{\omega }_{1x}=\sqrt{g{k}_{1x}tanh\left(k_{1x}{d}_0\right)}=5.764\mathrm{rad}/\mathrm{s},$$

(54)

$$k_{1y}={\displaystyle \frac{\pi }{W}},{\omega }_{1y}=\sqrt{g{k}_{1y}tanh\left(k_{1y}{d}_0\right)}=7.154\mathrm{rad}/\mathrm{s}.$$

(55)

The tank undergoes bi-directional lateral excitation:

$$a_x\left(t\right)=a_{0x}sin\left({\omega }_{ex}t\right),a_y\left(t\right)=a_{0y}sin\left({\omega }_{ey}t\right),$$

(56)

where ${\omega }_{ex}=0.9{\omega }_{1x}$ and ${\omega }_{ey}=0.9{\omega }_{1y}$, with $a_{0x}=X_0{\omega }_{ex}^2$ ($X_0/d_0=0.00186$) and $a_{0y}=Y_0{\omega }_{ey}^2$ ($Y_0/d_0=0.00120$). This produces two independent beating responses with periods $T_{\mathrm{beat},x}=2\pi /|{\omega }_{1x}-{\omega }_{ex}|\approx 10T_{1x}$ and $T_{\mathrm{beat},y}=2\pi /|{\omega }_{1y}-{\omega }_{ey}|\approx 10T_{1y}$.

6.2. 3-D Governing Equations

The velocity potential $\phi (x,y,z,t)$ satisfies the three-dimensional Laplace equation:

$${\displaystyle \frac{{\partial }^2\phi }{\partial x^2}}+{\displaystyle \frac{{\partial }^2\phi }{\partial y^2}}+{\displaystyle \frac{{\partial }^2\phi }{\partial z^2}}=0,0<x<B,0<y<W,-d_0<z<0.$$

(57)

The boundary conditions are extended as follows:

• KBBC ($z=-d_0$): $\partial \phi /\partial z=0$.
• KWBC (x-walls): $\partial \phi /\partial x=0$ at $x=0$ and $x=B$.
• KWBC (y-walls): $\partial \phi /\partial y=0$ at $y=0$ and $y=W$.
• KFSBC ($z=0$): $\partial \phi /\partial z=\partial \eta /\partial t$.
• DFSBC ($z=0$): $\partial \phi /\partial t+g\eta +a_x\left(t\right)x+a_y\left(t\right)y=0$.
• IC ($t=0$): $\phi =0$, $\partial \phi /\partial t=0$.

The DFSBC now includes a y-direction forcing term $a_y\left(t\right)y$, and the new KWBC on the y-walls is the only structurally new equation compared to the 2-D case.

6.3. 3-D Analytical Solution

Under linearised potential flow, the x- and y-excitations decouple into independent modal families. The x-excitation drives $(m,0)$ modes (odd m) and the y-excitation drives $(0,n)$ modes (odd n):

$$\phi ={\phi }_{0x}\left(t\right)+\sum _{m\mathrm{odd}}{A}_m\left(t\right)cos\left(k_{mx}x\right){\displaystyle \frac{cosh\left[k_{mx}(d_0+z)\right]}{cosh\left(k_{mx}{d}_0\right)}}+{\phi }_{0y}\left(t\right)+\sum _{n\mathrm{odd}}{B}_n\left(t\right)cos\left(k_{ny}y\right){\displaystyle \frac{cosh\left[k_{ny}(d_0+z)\right]}{cosh\left(k_{ny}{d}_0\right)}},$$

(58)

where $A_m\left(t\right)=C_m^x[cos\left({\omega }_{ex}t\right)-cos\left({\omega }_{m}t\right)]$ and $B_n\left(t\right)=C_n^y[cos\left({\omega }_{ey}t\right)-cos\left({\omega }_{n}t\right)]$ follow the same forced-oscillator structure as the 2-D case. The free-surface elevation is:

$$\eta (x,y,t)=-{\displaystyle \frac{1}{g}}\left[\sum _{m\mathrm{odd}}{A}_m^{\prime }\left(t\right)cos\left(k_{mx}x\right)+a_x\left(t\right)\left(x-\frac{B}{2}\right)+\sum _{n\mathrm{odd}}{B}_n^{\prime }\left(t\right)cos\left(k_{ny}y\right)+a_y\left(t\right)\left(y-\frac{W}{2}\right)\right].$$

(59)

The velocity components are $u=\partial \phi /\partial x$, $v=\partial \phi /\partial y$, $w=\partial \phi /\partial z$. Numerical verification confirms that this solution satisfies all governing equations and boundary conditions to within finite-difference truncation error ($\sim {10}^{-5}$). Importantly, the y-modes contribute approximately 42% of the velocity potential at the tank corners, confirming that the 3-D problem is not a trivial extrusion of the 2-D case.

6.4. 3-D Network Architecture

The 3-D network retains the same MLP structure (16 hidden layers, 64 neurons, tanh) with the following modifications:

Hard wall BC in both directions. The spatial inputs are embedded as $x_n=cos(\pi x/B)$ and $y_n=cos(\pi y/W)$. By the chain rule:

$${\displaystyle \frac{\partial \widehat{\phi }}{\partial x}}={\displaystyle \frac{\partial f_{\theta }}{\partial x_n}}·\left(-{\displaystyle \frac{\pi }{B}}sin{\displaystyle \frac{\pi x}{B}}\right),{\displaystyle \frac{\partial \widehat{\phi }}{\partial y}}={\displaystyle \frac{\partial f_{\theta }}{\partial y_n}}·\left(-{\displaystyle \frac{\pi }{W}}sin{\displaystyle \frac{\pi y}{W}}\right),$$

(60)

which vanish identically at $x\in \{0,B\}$ and $y\in \{0,W\}$ respectively. This eliminates all wall collocation points in both horizontal directions.

Extended Fourier embedding. The input is augmented from 9 to 16 dimensions to accommodate six physical frequencies:

$${\gamma }_{3\mathrm{D}}\left(\mathbf{x}\right)=\left(cos\frac{\pi x}{B},cos\frac{\pi y}{W},\frac{z+d_0}{d_0},\frac{t}{T_w},sin{\omega }_{1x}t,cos{\omega }_{1x}t,sin{\omega }_{ex}t,cos{\omega }_{ex}t,sin\Delta {\omega }_{x}t,cos\Delta {\omega }_{x}t,sin{\omega }_{1y}t,cos{\omega }_{1y}t,sin{\omega }_{ey}t,cos{\omega }_{ey}t,sin\Delta {\omega }_{y}t,cos\Delta {\omega }_{y}t\right),$$

(61)

where $\Delta {\omega }_x={\omega }_{1x}-{\omega }_{ex}$ and $\Delta {\omega }_y={\omega }_{1y}-{\omega }_{ey}$. The additional 7 input features add only $(16-9)\times 64=448$ parameters ($<1\%$ of total). Table 5 compares the 2-D and 3-D architectures.

6.5. 3-D Loss Function

The loss function retains the same structure and weights as the 2-D case, with the following modifications:

• 3-D Laplace:$r_{\mathrm{GE}}={\widehat{\phi }}_{xx}+{\widehat{\phi }}_{yy}+{\widehat{\phi }}_{zz}$.
• 3-D gradient-enhanced Laplace:$r_{\mathrm{GL}}={\left({\partial }_x{\nabla }^2\widehat{\phi }\right)}^2+{\left({\partial }_y{\nabla }^2\widehat{\phi }\right)}^2+{\left({\partial }_z{\nabla }^2\widehat{\phi }\right)}^2$.
• DFSBC:$r_{\mathrm{DFSBC}}={\partial }_t\widehat{\phi }+g\eta +a_x\left(t\right)x+a_y\left(t\right)y$ at $z=0$.

The loss weights (${\lambda }_D={\lambda }_K=100$, ${\lambda }_{\mathrm{GL}}=0.1$, ${\lambda }_I=10$) and training strategy (8,000 Adam + 2,000 L-BFGS) are unchanged from the 2-D case.

6.6. 3-D Collocation Grid

A structured grid of $N_x\times N_y\times N_z\times N_t$ interior points is used, with the same vertical clustering ($z=-d_0(1-{\xi }^{0.6})$) and early-time densification as the 2-D case. The default 3-D grid uses $N_x=21$, $N_y=11$, $N_z=21$, $N_t=101$, yielding approximately 490,000 interior points per window. The 15% mini-batch strategy keeps per-step GPU memory comparable to the 2-D case.

6.7. 3-D Results

The 3-D PINN was trained on a single NVIDIA RTX 5090 GPU. Training for one beating cycle took 2,499 s ($\approx 42$ min), comparable to the 2-D training time of ∼1,920 s on an A100.

6.7.1. Wave Elevation

Figure 11 shows the wave elevation time history at two tank corners: $(x,y)=(0,0)$ and $(x,y)=(B,0)$, where both x- and y-modal families contribute. The PINN tracks the analytical solution over one beating cycle with ${\epsilon }_{\eta }=0.24\%$.

6.7.2. Free-Surface Pattern

Figure 12 shows the three-dimensional free-surface elevation $\eta (x,y)$ at four time instants spanning one beating cycle ($t/T_{\mathrm{beat}}=0.10,0.25,0.40,0.50$). The analytical solution is rendered as a coloured wireframe mesh, while the PINN predictions are plotted as scatter points at the same grid nodes. The points sit precisely on the mesh surface at all four instants, confirming sub-percent accuracy. The bi-directional excitation produces a genuinely two-dimensional wave pattern: the free surface tilts in both x and y simultaneously, with peak amplitudes of $\sim 12$ mm at the corners growing from $\sim 5$ mm at $t=0.10T_{\mathrm{beat}}$ to the maximum at $t=0.25T_{\mathrm{beat}}$, then decreasing as the beating envelope passes through zero.

6.7.3. Velocity Fields

Figure 13 compares the velocity fields at the $y=W/2$ mid-plane at $t=T_{\mathrm{beat}}/4$. All three velocity components (u, v, w) are shown. The transverse velocity $v=\partial \phi /\partial y$ is a purely 3-D quantity with no 2-D counterpart; at $y=W/2$ the transverse modes have $sin(n\pi y/W)=\pm 1$ (odd n), so v reaches its maximum magnitude at this plane, providing a stringent test of the network’s ability to learn the y-modal contributions.

6.7.4. Error Metrics

Table 6 summarises the 3-D PINN accuracy over one beating cycle. The errors are computed over the full 3-D volume, not on a single slice.

Comparing with the 2-D results (${\epsilon }_{\eta }=0.15\%$, ${\epsilon }_u=0.84\%$, ${\epsilon }_w=1.65\%$), the 3-D errors are expected to be moderately higher due to the increased input dimensionality and the additional PDE terms. The transverse velocity error ${\epsilon }_v$ is reported for the first time and reflects the network’s ability to learn the y-modal family.

7. Conclusions

This paper has applied PINNs to near-resonant excited sloshing with a beating response in both two-dimensional and three-dimensional rectangular tanks. For the 2-D case, two training obstacles specific to forced sloshing are identified and addressed. Residual normalization and loss weighting (${\lambda }_D={\lambda }_K=100$) eliminate the zero-solution trap; Fourier time embedding at the three physical frequencies (${\omega }_1$, ${\omega }_e$, $\Delta \omega$) resolves spectral bias. A hard wall boundary condition via the $cos(\pi x/B)$ embedding enforces the lateral wall condition exactly, and a gradient-enhanced Laplace regulariser ($\nabla \left({\nabla }^2\phi \right)=\mathbf{0}$) improves velocity accuracy. Together, these techniques yield ${\epsilon }_{\eta }=0.15\%$, ${\epsilon }_u=0.84\%$, and ${\epsilon }_w=1.65\%$ over one beating cycle, validated against the viscous finite-difference results of Chen & Nokes [4].

A frequency parameter study covers ${\omega }_e/{\omega }_1=0.5$–$1.1$, including the near-resonance regime. A time-domain decomposition strategy with transfer learning partitions the domain into one-beat windows. For the 3-beat case ($30T_1$), decomposition reduces velocity errors from ${\epsilon }_u=15.79\%$ (single network) to $3.26\%$. Extending to 5 beats ($50T_1$), overall accuracy is maintained at ${\epsilon }_u=3.43\%$ and ${\epsilon }_{\eta }=0.30\%$ with no monotonic error accumulation.

The extension to three dimensions demonstrates that the same techniques—Fourier embedding, hard wall boundary conditions, and gradient-enhanced regularisation—transfer effectively to higher-dimensional sloshing. The 3-D formulation introduces bi-directional excitation that excites both $(m,0)$ and $(0,n)$ modal families, producing a genuinely three-dimensional free-surface response with y-modes contributing approximately 42% of the velocity potential at the tank corners. The 16-dimensional Fourier embedding accommodating six physical frequencies, combined with hard wall enforcement in both x and y directions, enables the network to capture the additional modal family with errors of ${\epsilon }_{\eta }=0.24\%$, ${\epsilon }_u=1.31\%$, ${\epsilon }_v=1.78\%$, and ${\epsilon }_w=2.32\%$ over one beating cycle.

Future work will address nonlinear large-amplitude sloshing, complex tank geometries, broadband excitation spectra, and the coupling of time-domain decomposition with the 3-D formulation for long-time multi-beat simulations.