Numerical simulations of shallow-water sloshing coupled to horizontal vessel motion in the presence of a time-dependent porous baffle

– Shallow-water fluid sloshing in the Lagrangian Particle Path formulation, with the addition of an energy-extracting porous baffle, is simulated numerically using a symplectic numerical scheme which captures, in an essential way, the energy exchange. The fluid motion in a rectangular vessel is dynamically coupled to a surface-piercing porous baffle. The fluid transmission through the baffle is characterized by a nonlinear Darcy-Forchheimer model equation. The numerical scheme is symplectic, based on the implicit-midpoint rule, and thus is strategically designed to maintain the energy partition between the fluid and vessel throughout numerous time steps. Our results demonstrate the non-conservative nature of the system, with the porous baffle effectively dissipating energy from the overall system. Furthermore, we present findings that demonstrate the role of time-periodic variations in baffle porosity on energy dissipation. By manipulating the frequency and magnitude of this time-dependent variability, it is established that a greater amount of energy can be extracted from the system compared with the optimal fixed porosity baffle. These results shed new light on potential strategies for enhancing energy dissipation in such configurations.


Introduction
Research into fluid sloshing motions in vessels forced to move in prescribed manners is a vast area of research in both science and technology (see [1] and [2] and references therein).The motion of the fluid is typically complex, and as the fluid interacts with the walls of the vessel it generates forces and moments on these walls which can have unintended consequences, such as causing faults in fuel tanks which have been excited by phenomena such as earthquakes [3].Other possible scenarios where fluid sloshing may occur might include maritime and terrestrial fuel transport [4].If the vessel is free to move (either fully, or in some constrained way) then those fluid forces and moments can generate movement of the vessel, which in turn feedback into the motion of the fluid.This coupled motion is very complex, and can lead to devastating destabilization consequences, such as capsizing King crab boats [5,6] and the spilling of coffee from a mug while walking [7].In order to mitigate against destabilization it is important to recognize when destabilization is occurring and put measures in place to control the system.For example, a 'carry cradle' has been shown to be very successful in reducing fluid sloshing in carried coffee mugs [8].The novel feature of this paper is we consider a time-dependent porous baffle as a mechanism to control the system.
Baffles that are non-porous have been utilized in sloshing systems to reduce the negative consequences of the sloshing fluid, essentially by blocking the flow [9,1].This blocking of the flow leads to a change in the natural frequency of the system which can alter the forces generated inside the vessel, making them less severe.Such a scenario was considered by [10] who considered a rectangular vessel with N surface-piercing baffles and found the natural frequency of the system increased as N was increased.Porous baffles can also modify the system's natural frequency [11], but they play a more significant role by extracting energy from the system and damping the motion [12,13].There have been significant experimental and numerical studies on the effect of porous baffles on forced systems (not dynamically coupled), in particular focusing on their position, size and material construction [14,15,16,17,18,19,20,21,22,23].
The key features of the baffles turn out to be their size and construction.The baffle construction, i.e is it a randomly drilled porous block or perforated plate, is particularly significant to its overall damping properties.In numerical simulations of the baffle systems in the literature, the flow though the baffle is not explicitly calculated, and instead fluid transmission through the baffle is modelled via a model equation.Such an equation can include using Darcy's law in linear and weakly nonlinear scenarios [24,25,26,27], using a nonlinear pressure drop condition for highly nonlinear scenarios [28,29] or using a Darcy-Forchheimer condition [30] for intermediate conditions.The Darcy-Forchheimer condition has proved particularly successful in modelling the sloshing motion of a fluid in a tuned liquid damper filled with a porous media [31,32].Turner [11] examined the linear normal modes for the coupled sloshing problem, and investigated how the frequency of these modes varied with the system parameters.
In the current paper we perform, to the best of our knowledge, the first numerical simulations of a dynamically coupled system, to include porous baffles.To model the system, we approximate the fluid as a shallow-water layer, and use the Lagrangian Particle Path (LPP) approach, developed for dynamic sloshing systems by [33].This approach has been widely implemented in similar dynamic sloshing problems [34,35,36] due to the fact that the scheme has a symplectic structure and, hence, can be integrated by a geometric integration scheme such as the implicitmidpoint rule [37,34,38].Such a scheme preserves the partition of energy between the vessel and the fluid, and bounds any numerical leakage of energy from the fluid to the vessel.This is particularly important when we allow the porosity of the baffle to vary with time.
The current paper is laid out as follows.In §2 we formulate the governing equations for the nonlinear system and then take the shallow-water approximation of these equations.In §3 we consider the linear solution to these equations for scheme validation, while §4 formulates the numerical scheme in both the nonlinear and linear regimes.Results of the scheme are presented in §5.Linear validation results and weakly nonlinear simulations are found in §5.1, while a mode switching phenomena in the limit as the fluid layer depth goes to zero is presented in §5.2.Section 5.3 investigates the control properties on the system when a baffle with timedependent porosity is considered.Concluding remarks can be found in §6.
2 Formulation of Governing Equations

Finite Depth Fluid Equations
The problem we consider is that of a two-dimensional vessel with rectangular cross-section.The vessel is of length L and contains a surface piercing porous baffle at x = L 1 , while the impermeable rigid side walls of the baffle lie at x = 0 and x = L = L 1 + L 2 .The vessel is partially filled with an inviscid, incompressible, irrotational fluid of constant density ρ, which has a mean depth of H, meaning that the vessel is filled with a fluid of mass m f = ρHL.Figure 1 depicts a schematic of the problem of interest, where variables are indexed by subscripts 1 and 2 in the left and right compartments respectively.
Figure 1: Schematic diagram of the problem under investigation.The rectangular vessel undergoes rectilinear motion in the X-direction and is partially filled with an inviscid, incompressible fluid which is allowed to pass through the porous baffle at x = L 1 .The vessel restoring force is given by the linear spring connecting the vessel to X = 0.
In the spatial frame the coordinates are X = (X, Y ), while in the body frame fixed to the vessel, the rectangular coordinates are x = (x, y) as shown in figure 1.The vessel undergoes unidirectional motion parallel to the x-direction, of amplitude q(t) and hence the two reference frames are related via where (X 0 , Y 0 ) is a constant displacement.The vessel is attached to a wall at X = 0 by a linear spring, which acts as a restoring force on the vessel, hence the function q(t) is determined as part of the solution.In each compartment the fluid occupies the region for j = 1, 2, where x 0 = 0 and x 2 = L.We represent the fluid velocity field via the two-dimensional vector (u j (x, y, t), v j (x, y, t)) and the pressure field via the scalar quantity p j (x, y, t).Thus the Eulerian representation of the fluid momentum equations and the conservation of mass equation in each compartment, relative to the body coordinates, are given by the Euler equations for j = 1, 2.Here g is the gravitational constant, is the twodimensional material derivative and the overdots denote derivatives with respect to t.
On the free surface of the fluid in each compartment, the pressure and velocity fields satisfy the dynamic and kinematic boundary conditions, namely p j = 0, and ∂h j ∂t + u j ∂h j ∂x = v j , at y = h j (x, t) , for j = 1, 2 . (2. 2) The boundary conditions on the rigid side-walls and vessel bottom are while the baffle at x = x 1 = L 1 is porous and thus fluid can be exchanged between the two compartments.We do not explicitly calculate the flow through the baffle, and instead we model the transmission of the fluid through the baffle with a model equation.We assume that at the baffle the fluid velocities and the fluid pressure satisfy a Darcy-Forchheimer model [30] and hence these quantities are related via (2.4) The coefficient β ∈ R is a measure of the porosity of the baffle and γ ∈ R is known as the internal permeability.If β or γ are zero, with the other nonzero, then we have an impermeable baffle [34], while if β and γ are infinite then the baffle is no longer present, and we now have a single compartment vessel [33].
Using (2.2) and the second equation of (2.1), we can modify the right-hand-side of the transmission boundary condition (2.4) at the baffle such that it is written in terms of the surface elevations h j rather than the pressures p j .To do this we integrate the second equation of (2.1) with respect to y to give using the dynamic boundary condition.Thus This is an exact equation for the pressure difference between the two compartments, which will simplify further once we make assumptions about the depth of the fluid layer in §2.2.
Rather than using the velocity field throughout the whole fluid, we instead choose to express the momentum equations (2.1) in terms of the surface velocity field, defined by This transformation is laid out in detail in [33] for the single compartment vessel, and so here we just quote the result, which is that the surface velocity components in each compartment satisfy the two exact equations at y = h j (x, t) for j = 1, 2, where (2.5) has been used to eliminate the pressure.
The motion of the vessel is determined by a forced pendulum equation and is derived using a variational principle.This derivation was laid out clearly for the impermeable baffle problem in [34], and hence we just state the equation as where m f is the fluid mass.Note, in this work we only consider a linear spring equation, but a nonlinear spring can easily be incorporated by including a −ν 2 q 3 term on the LHS of (2.9) [34].

Shallow-Water Approximation to Governing Equations
The exact set of equations for the motion of the system (2.8) and (2.9) are not closed because the terms on the right-hand sides contain terms involving v j (x, y, t), u j (x, y, t) and V j (x, t).In order to close this set of equations we need to neglect the right-hand sides of (2.8) and (2.9), which can be achieved by making the shallow-water assumption for our problem.The principle assumptions in this limit are then for j = 1, 2, where U 0 is an order one reference velocity.These assumptions are the usual shallow-water fluid theory assumptions.Namely, that the vertical Lagrangian acceleration at the free-surface is smaller than the gravitational restoring force, and that the horizontal velocity is independent of the vertical coordinate y.The first condition in (2.10) is a special case of the more common shallow-water assumption that the Lagrangian acceleration in the vertical direction is small everywhere [39], while if we assume ∂u j /∂y = 0 in each compartment then integrating the continuity equation leads to v j + y∂u j /∂x = 0 which leads to the second condition when evaluated on the free-surface.The third condition is then a consequence of this assumption.
Applying these assumptions to our governing equations from §2.1 gives the shallow-water equations for the fluid motion as for j = 1, 2, and these are coupled to the vessel motion via the forced pendulum equation which is written as d dt (2.12) In the shallow-water regime the rigid side-wall boundary conditions and the porous baffle transmission condition become The governing shallow-water equations above can be determined by the Euler-Lagrange equations of the Lagrangian functional where λ i are Lagrange multipliers.This Lagrangian is of the form where the kinetic (KE) and potential (P E) energies are with the constraints imposed.It is also possible to show that the total energy of the system decays with time, see appendix A.
In the next section we consider the solution of the shallow-water equations in the linear regime of small amplitude motions.

Linear Coupled Problem
The solution of the shallow-water problem in the linear regime of small amplitude motions allows us to determine the natural frequencies of the system, i.e. those frequencies at which the system wants to naturally oscillate.To consider the small amplitude motions of the system, we seek small velocity perturbations about a quiescent free-surface at y = H and hence the governing shallow-water equations in each compartment reduce to for j = 1, 2, with the boundary conditions given by (2.13) and (2.14).The linearised form of (2.14) for small velocities about a quiescent system is while the coupled linear vessel equation, linearised about q = 0 reduces to d dt To identify the characteristic equation for this system of equations, we consider a normal mode decomposition of the perturbation quantities, hence we let where c.c denotes the complex conjugate.This leads to the linearised set of equations iω h j + H U jx = 0, iω U j + g h jx = ω 2 q, for j = 1, 2, where here the subscripts denote partial derivatives.Eliminating h j from each pair of equations leads to a governing equation for each U j as where α 2 = ω 2 gH .
Solving this expression for j = 1, 2 gives for unknown constants A j , B j .The two side-wall conditions from (2.13) imply that while the second two equations lead to the pair of equations where L 2 = L − L 1 is the length of the second compartment.
The case when sin αL 1 = sin αL 2 = 0 was examined in [11] and this case leads to symmetric sloshing modes in a stationary vessel, i.e. modes which do not link to the vessel motion.For the case when these quantities are non-zero the above pair of equations are solvable, giving The linearised vessel equation gives a third equation linking B 1 , B 2 and q.To evaluate this equation we use the fact that giving the third equation as Writing these three equations as a matrix system and solving for non-trivial solutions leads to a characteristic equation of the form The above formulation gives the natural frequencies of the system.These natural frequencies give the eigenbasis for the nonlinear formulation of the problem, i.e. nonlinear solutions are superpositions of these linear eigenmodes.These linear solutions are also vital for providing initial conditions allowing us to validate the numerical scheme against known results.We discuss the solutions to (3.22) in §5.1 to which there are an infinite number of modes, the lowest frequency ones of which are plotted in figure 2.

Numerical Scheme for Nonlinear Simulations
The numerical scheme used to solve the governing shallow-water equations is based on a Lagrangian Particle Path (LPP) formulation [33], and is very similar to that laid out for the impermeable baffle problem in [34].The reader is referred to [34] for more details of the formulation.
The LPP formulation comes from converting from Eulerian to Lagrangian coordinates via the non-degenerate mapping which when substituted into (2.11) and (2.12) gives for 0 < a < a B for j = 1 and a B < a < L for j = 2.Here a B is the unknown Lagrangian coordinate which fixes the position of the baffle in the Eulerian framework and needs to be determined as part of the solution.The vessel equation becomes where comes from the mass conservation equation and a 0 = 0, a 1 = a B and a 2 = L.Here χ (j) (a) is independent of τ and hence is fixed at τ = 0. From this equation we are able to determine the free-surface elevation via h j = χ (j) /x a .In this LPP formulation the unknowns are x(a, τ ) and q(τ ).
The above equations (4.24) and (4.25) in Lagrangian coordinates can be determined from a variational principle from the Lagrangian (2.15) where we use (4.23) to convert the Lagrangian to these Lagrangian coordinates.Converting the Lagrangian and performing a Legendre transformation leads to a non-autonomous Hamiltonian functional where the constraints hold and (x, w, q, p) are the canonical variables with momentum variables The governing equations can be written in Hamiltonian form as The impermeable side wall boundary conditions to be solved together with (4.29) are x(0, τ ) = 0, and and For the conditions at the porous baffle, the transmission conditions are and from (2.4) Note that the introduction of the condition (4.33) means that this system, unlike the impermeable baffle problem in [34], is no longer conservative and thus energy is not conserved.However, it has been shown that symplectic integration schemes work very well for symplectic equations with dissipative perturbations, even capturing the actual dissipation better than non-symplectic schemes [38].They have also been shown to be very efficient at preserving the energy budget between the fluid and the vessel [33], hence we choose to use a symplectic scheme here.
4.1 The Semi-Discretization of (4.29) The challenging numerical aspect of solving the system (4.29)-(4.34) is that the Lagrangian coordinate a B is a function of τ , and so the domains on which the governing equations are solved is not fixed.To overcome this we make the following change of variables to a fixed domain in each region, namely where ξ 1 (τ ) = a B (τ ) and ξ 2 (τ ) = L − a B (τ ).
The domains A (1) and A (2) can then be regularly discretized by setting where N and M are integers and we define ∆A (1) = 1/N and ∆A (2) = 1/M .Notationally we write x i , τ ) where j = 1, 2 denotes the region of the vessel.
The change of variables and discretization of the equations (4.29)-(4.34) is relatively straightforward, except for the w τ equation, which requires a variational discretization, where the Hamiltonian function (4.27) is first discretized and then variations of this are then taken.In this case, following [33], these equations become Thus the full semi-discretization of the transformed equations are then given by (4.37) and (4.38) together with i − 39) where i = 2, . . ., N for j = 1 and i = 2, . . ., M for j = 2, and where (1) N +1 ∆A (1) , σ (2) M +1 w (2) are the integral terms in (4.29) which have been discretized using the trapezoidal rule.
Under the change of variables (4.35) and (4.36), equation (4.26) no longer states that χ is independent of τ , and in fact the conservation of mass equation now states that we must additionally solve with the corresponding forward and backward forms of ∂x ∂a to evaluate the four end points i = 1, N + 1 for j = 1 and i = 1, M + 1 for j = 2. Now the free surface elevation is given by

Time Discretization using the Implicit Midpoint Rule
The semi-discretization of the governing equations above have the form where p = (p, w 1 , . . ., w N +1 , w 1 . . ., w M +1 ), q = (q, a B , x We discretize τ such that τ n = n∆τ where ∆τ is a fixed time-step, and introduce the notation q(τ n ) = q n .Then, equations of this form can be integrated by applying the implicit midpoint rule [38], which is second order accurate in time and discretizes the equations as (4.49) The implicit nature of this algorithm means that larger time-steps, ∆τ , can be taken compared to an explicit scheme, which helps in keeping the scheme computationally fast when weakly nonlinear simulations are considered which require a higher spatial resolution.
The system of nonlinear algebraic equations is solved using Newton's method, where iterations are continued until

Numerical Results
In this section we present numerical results to the nonlinear equations (4.39)-(4.45).One of the limitations of the shallow-water model used in this paper is, for sloshing problems the equations are only able to sustain a weak amount of nonlinearity before wave-breaking occurs.Thus we are limited to parameter regimes where wave overturning does not occur and the nonlinearities observed are weak.In cases where the scheme approaches breaking, the code breaks down due to the free-surface becoming multiply defined.Therefore, we find that at most, our results are weakly nonlinear.Hence we choose to fix γ(τ ) = 1 in the nonlinear transmission condition throughout this section, and find no quantitative difference to the results presented if this value is significantly varied.
In all the results presented in this manuscript the parameters N, M and ∆τ , which govern the resolution of the numerical scheme, were varied to ensure the presented results have converged.

Scheme Validation and Weakly Nonlinear Example
In this section we validate our numerical scheme by making direct comparisons with the linear results presented in §3.Before we compare the numerical simulation results, we first consider what frequencies of solution, ω = ω r + iω i , we expect to observe by solving the characteristic equation (3.22).Note that the value of ω i denotes the exponential decay rate of the mode as we only find solutions with ω i ≥ 0. In figure 2  with N = M = 100 and ∆τ = 10 −3 s.The results show that those modes for which ω r varies considerably as β → 0, also have the largest decay rates.The lowest frequency modes (labelled 1) are likely to be the modes observed in a physical system, as these modes typically decay the slowest.More discussion on the modal results and how the various system parameters affect these results can be found in [11], while more information on this type of eigenvalue spectrum for similar coupled fluid/vessel systems without baffles can be found in [40,41,42,43].
In figures 3-6 we compare the numerical (solid black lines) and the exact analytical results (dashed red lines) for the three cases noted above.In all cases, the exact and analytical results for the vessel displacement q(t) and the free-surface elevations h(x, t) are in excellent agreement.Note that we plot our results in terms of the Eulerian variables (x, t) as these are more intuitive than (A (j) , τ ).In this linear case t = τ and x = [ξ 1 A (1) , (a B + ξ 2 A (2) )].The mode 6 result has ω = 10.225+ 2.529i rad s −1 and so decays very quickly compared to the other two cases where ω = 1.027 + 0.032i rad s −1 and ω = 1.019 + 0.023i rad s −1 respectively.However, even in this case the agreement for q(t) and a B (t) − L 1 in figure 3(b) is excellent.For the free-surface profiles, h(x, t), in figure 5 we see some numerical round off errors in the black numerical results, but the overall agreement is still excellent.
As q is increased the problem becomes more nonlinear, however the parameter window to observe the nonlinearity is limited, because if we increase q too much, then we get wave breaking and the numerical scheme breaks down.Hence we can only realistically model weakly nonlinear situations with this scheme, which we consider in figures 7 and 8 with q = 0.05 m.In these figures the spatial resolution has been increased to N = M = 400 to give converged results.
In figure 7 we consider the evolutions of q(t) and a B (t) where we have subtracted off the equivalent linear solution functions.The results show small O( q 2 ) corrections to these values, q − q linear a B − a linear B t Figure 7: Plots of the nonlinear vessel displacement from the linear result q(t) − q linear (t) and the nonlinear Lagrangian value at the baffle difference a B (t) − a linear B (t) for the parameters in (5.50) with L 1 = L 2 = 0.5 m, β = 0.2, sm −1 and q = 0.05 m and with the resolution N = M = 400.Here t has units of s while q and a B have units of m. with higher frequency oscillations appearing for the duration of the simulation.We saw in figure 2 that these higher frequency modes typically have faster decay rates than the lowest frequency mode and so these modes decay away as time evolves.The nonlinear terms grow in magnitude over t ∈ [0, 50] s before decaying as the porous baffle removes the energy from the system.The weakly nonlinear free-surface profiles in figure 8   nonlinear oscillations, and in these appear to persist throughout the whole simulation.Now that the numerical scheme is validated, we go on to highlight an interesting mode switching phenomena in the next section, before investigating the effects on the system of incorporating a baffle with time-dependent porosity.

Mode Switching Phenomena as H → 0
In the zero-baffle problem considered by [44] they found a curious feature of the system where the dominant behaviour in the system switches to ever higher frequency modes, as the fluid depth H → 0, in order that the correct dry vessel oscillation frequency is achieved in the H = 0 limit.In the porous baffle case we would expect the same behaviour to occur, but it is less clear as to how the decay rate of the system would behave as this limit is approached.This is investigated in this section.
The problem [44] considered is slightly different to that in figure 1, and actually comprises of a vessel hung as a bifilar pendulum with string lengths l, see figure 9, where here q(t) denotes the angle of the string to the vertical.However, in the linear regime, where this mode switching phenomena is observable, both problems are identical if the linear spring coefficient is allowed to vary with m f (H) via ν = g l (m v + m f ). (5.58) Figure 9: Schematic diagram of the bifilar pendulum problem considered by [44], with the addition of a porous baffle at x = L 1 , which demonstrates the model switching phenomena examined in §5.2.In the linear regime this problem is equivalent to that in figure 1 We use (5.58) together with the system parameters ) with N = M = 100 and ∆τ = 10 −3 s, which match those parameters used by [44] in their calculation.In this section we consider an initial condition such that q(0) = q + c.c., x (1) (A (1) , 0) = ξ 1 A (1) , x (2) (A (2) , 0) = a B + ξ 2 A (2) , w (j) (A (j) , 0) = 0, for j = 1, 2 This initial condition is one which is more readily observed in an experimental setup, and comprises of a condition where the vessel is displaced from its equilibrium position, extending the spring or placing the hanging string at a non-zero angle, and then being released from rest.The solution in this case comprises of a superposition of all the modes of the system, with differing amplitudes.
We conduct a series of numerical simulations and from each of these we consider the values of ω r and ω i as a function of m f /m v (which is equivalent to varying H for fixed m v , ρ, L 1 + L 2 ) in figure 10.In this figure the five black lines signify the values of the five lowest frequency modes found by solving the characteristic equation (3.22), while the circles are estimates taken from the numerical simulations.The dashed line is the dry vessel frequency g/l, and hence is the system frequency at m f = 0. To calculate the values of the circles we compute the vessel motion q(t) from t = 0 to t = 10 s.
We then use the MATLAB subroutine findpeaks to identify the local maxima of this function, and then use each pair of successive maxima to estimate a value for ω r and ω i .We denote the numerical estimates by ω num r and ω num i respectively.Finally we average these values over all maxima in [0, 10] s to get the estimate value we plot in figure 10.
These numerical results show the frequency of the simulations switching from the mode labelled 1 to the mode labelled 2 around m f /m v ≈ 0.6.The peculiar feature of the result here is that this transition does not occur exactly at the point where the two ω i values intersect.The numerical results follow mode 1 beyond this point (m f /m v ≈ 0.7) and then switch.What this switching actually amounts to is the time-dependent amplitude of the modes switching dominance, such that mode 2 has the largest amplitude for m f /m v ≲ 0.5, as the numerical solution is a superposition of all these modes.As m f /m v is reduced below 0.5, eventually the numerical results agree with the dry vessel result given by the dashed line.The decay rate of the solution reduces as m v /m f → 0, but very close to zero (m f /m v ≈ 0.04) there is a small increase in ω i before going to zero in the dry vessel limit.This mode switching in the function q(t) can be seen in figure 11.In figure 11(a) the value of m f /m v = 1.0 and it is clear that the frequency and decay rate agree very well with the exact linear result of mode 1, given by the red dashed line.For m f /m v = 0.5 in figure 11(b) we are in the transition region between the two modes and clearly there are oscillations of a shorter wavelength (larger ω r ) emerging, but there is a modulation similar to the frequency of the mode 1 result.Then at m f /m v = 0.25 in figure 11(c) it is clear that both ω r is larger and ω i is smaller and in this figure the red dashed line gives the exact result for mode 2, showing excellent agreement.Such a switching phenomena could be significant in oscillating systems where the fluid is slowly being drained out.In such cases, being able to change the porosity of the baffle remotely might then be desirable, in order to control the system.In the next section we consider baffles with The black lines are the result of the numerical simulations, while the red dashed lines are the exact linear result for the mode labelled 1 in figure 10(a,b) and for the mode labelled 2 figure 10(c).Here t has units s and q has units m.
time-dependent porosity to try to understand how such systems behave, and whether they allow for the system to be controlled.

Time-Dependent Baffle Porosity
In this section we again consider the system from figure 1, but we investigate the effect of prescribing the porosity of the baffle to vary with time.Such a baffle could consist of a regular holed slat screen [29], where the holes are mechanically opened and closed in some prescribed manner.The interest in this problem is to try and understand how the system behaves, with the main interest in whether the decay rate of the system can be increased from the maximum decay rate for the lowest frequency mode, by periodically oscillating the baffle's porosity.
For the results in this section we again consider the system parameters (5.50) with L 1 = L 2 = 0.5 m and q = 5 × 10 −4 m.We consider five periodic forms for where Here τ has units s, and β has units sm −1 .
The forms of these functions are plotted in figure 12.The function β 1 (τ ) is chosen as it uses the whole range of value of β ∈ [0, ∞), while for β 2−5 (τ ) the time variations are a simpler periodic function with constant amplitude, and hence are more likely to be manufacturable and implementable.
In figure 13 we consider the numerical decay rate ω num i for simulations with β = β 1 (τ ), in the (Ω, T on )-plane.The values of ω num are calculated as in §5.2 except here we average over the time period τ ∈ [0, 100] s.The black contour corresponds to the maximum linear decay, ω max i = 0.03218 rad s −1 , rate for a fixed baffle with β = β max = 0.508 s m −1 .The lighter regions (green-yellow) signify parameter regions where the system decay rate is above ω max i , i.e. the system decays faster than would be seen in a system with a fixed baffle with β = β max .What this figure shows is at low frequencies Ω ≲ 8 s −1 there is some dependence on the T on value when the baffle porosity switches on, but above this frequency the results are independent of this What we believe is happening here is a resonance between β(τ ) and the frequency ω r (τ ) of the dominant mode in the system.If the two values of these are such that the velocity of the fluid is relatively low at the baffle when the baffle is closed or partially open, then very little energy is removed from the system (hence a low ω num i ), whereas if the velocity at the baffle is high at these times then the energy removed is larger (hence higher ω num i ).
In figure 14 we consider plots of ω num i (Ω) for fixed T on = 2 s for β 1 (τ ) to β 5 (τ ).For the β 1 (τ ) case in figure 14(a) this corresponds to a horizontal slice though figure 13.This figure shows clearly those values of Ω where ω num i > ω max i and it is clear that the magnitude of the peak decay rate increases with Ω.When we consider β 2 (τ ) in figure 14(b) we find that overall, the decay rate is lower, with no significant frequencies with ω i > ω max i until Ω > 10 s −1 , and even then, these regions are small.It is again true that the magnitude of the peak decay rate increases with Ω.For β 3 and β 4 in panels (c) and (d) the amplitude of the porosity function is, on the whole, smaller and larger respectively than β 2 .What is curious is that the smaller amplitude case in (c) has significant values of Ω where ω num i > ω max i , whereas the larger amplitude function in (d) has predominantly ω num i < ω max i .The reason for this could be down to the fact that β 1 (τ ) and β 3 (τ ) have a similar asymptotic behaviour close to τ − T on = 2nπ/Ω, n ∈ Z, as can be seen in figure 12. Overall, the peak values of ω num i are smaller in figure 14(c) than figure 14(a), and so opening up the baffle completely (β → ∞) is important, but the most significant aspect is the rate at which the baffle is opened from the closed (β = 0) position.In figure 14(e Ω here. While the porosity function β 1 appears to perform the best at controlling the decay rate of the coupled system, it is unlikely that such a baffle could be manufactured which has β ∈ [0, ∞), and it is more likely that a mechanical system with a more sinusoidal porosity would be possible.To this end we examine the porosity function (5.62), but we allow A to vary in figure 15 for In this case the maximum decay rate is 86% larger than the dashed line result.In figure 15(c) we consider the evolution of q(τ ) for the case when β = β max and the two amplitudes giving the maximum values of ω num i in panels (a) and (b).The results demonstrate the faster decay rate of the time-dependent baffle, in particular highlighting the significant decay of the system for Ω = 25 s −1 (black line).Panel (d) gives a zoomed in plot of the final 50 seconds of the simulation where the changes in vessel amplitude are more clearly observed.
When we consider a larger range of (Ω, A) values in figure 16, we see that the largest decay rates are concentrated in small regions of Ω around Ω = 15 s −1 and Ω = 25 s −1 with the peak amplitudes close to A = 2 s m −1 , akin to the figure 15(b) result.As well as the small regions of large decay rates, there are also larger regions, which are periodic in Ω, for amplitudes A ∈ [1, 2] m s −1 .These regions are more akin to the result in figure 15(a) where the maximum decay rate is only about 20 − 25% larger than the dashed line.What these results show is that it is entirely possible to control the decay rate of the coupled sloshing system in figure 1 using a time-dependent baffle, at least at the lower fluid velocities considered here, where (2.14) is valid.Clearly, changing the porosity value of the baffle from one constant value to another constant value is one way to control the system, but if you wish to damp out the system oscillations faster than would occur with a fixed baffle, then periodically oscillating the porosity of the baffle is the best approach, although the amount of variation and the frequency of the variation depends carefully upon the parameters of the system.

Conclusions
In this work we considered the coupled vessel-fluid dynamic system.In the system a rectangular vessel partially filled with an inviscid incompressible fluid was allowed to move in a rectilinear motion connected to a linear spring, and the fluid was allowed to slosh back and forth in the vessel.The fluid and the vessel exhibit forces on one another, giving complex coupled motions.The main focus of this paper was to model these complex motions when a vertical, surfacepiercing, porous baffle is inserted into the vessel, which was found to extract energy from the system causing the motion to decay in time.The fluid transmission condition at the baffle was modelled using a nonlinear Darcy-Forchheimer model to more accurately deal with the nonlinearity of the porous baffle in weakly nonlinear scenarios.
The complex motion was determined numerically via a fast and efficient numerical scheme which approximated the fluid layer as a shallow-water fluid and used a Lagrangian Particle Path approach.This approach was used together with a symplectic time-integration scheme based on the implicit-mid-point rule as this approach preserved the energy partition between the vessel and the fluid and effectively captured the dissipation.The implicit nature of the scheme also allowed for larger time-steps than an explicit scheme, such as Störmer-Verlet.Using an implicit scheme is beneficial because a larger time step can be used (typically 10-100 times larger than for an explicit scheme) which, despite the added iterations of the implicit scheme, generally make it faster than an explicit scheme.
Results were presented for a single baffle system, which splits the vessel into two compartments, in both the linear and weakly nonlinear regimes and excellent agreement was found with the exact linear solution for various test cases.Simulations were presented for a scenario where the baffle was assumed to mechanically open and close, such that its porosity properties varied in a time-periodic manner.It was shown that under such circumstances the average decay rate of the system over a fixed period of time, can be larger than the maximum decay rate for the lowest frequency solution for a fixed porosity baffle.This showed that time-dependent baffles are an effective mechanism to control the decay rate of the coupled system.x j−1 ρ (h j U 2 j ) x + 1 2 g(h 2 j ) x + h j q dx, = −m f q − 2 j=1 x j x j−1 ρ (h j U 2 j ) x + 1 2 g(h 2 j ) x dx, and so (2.12) can be expressed as To show that the baffles are dissipative we note that dE dt = 2 j=1 x j x j−1 1 2 h jt (U j + q) 2 + h j (U j + q)(U jt + q) + 1 2 gh j h jt ρ dx + (m v q + νq) q, = 2 j=1 x j x j−1 1 2 (h j U 3 j ) x + q2 (h j U j ) x + g(h after using (A-1), (2.11) and integrating with respect to x.Here we have used the fact that x 1 = L 1 and that U 1 (x 0 = 0, t) = U 2 (x 2 = L, t) = 0.
In order to demonstrate that the quantity on the right-hand-side is negative, we consider two limits of (2.13), the small velocity limit and the large velocity limit, with the understanding that the general case lies in between.

A.1 Small velocity case
In this case (2.13) simplifies to (3.18)This term is clearly negative definite and hence the baffle is dissipative.

A.2 Large velocity case
In this case (2.13) simplifies to Hence from this we deduce that and thus for j = 1, 2.
Inserting this into (A-2) gives which is again negative definite.

Figure 2 :
Figure 2: Plots of (a,c) normal mode frequencies ω r (β) and (b,d) normal mode decay rates ω i (β) from (3.22) for (a,b) L 1 = L 2 = 0.5 m and (c,d) L 1 = 0.25 m, L 2 = 0.75 m.The other parameters are fixed and given in (5.50).The circles at β = 0 signify modes which are symmetric modes in a stationary vessel in this limit.Here β has units sm −1 and ω r , ω i have units rad s −1 .

Figure 8 :
Figure 8:Free-surface profiles h(x, t) for the case in figure7at t = 12.5, 25.0, 37.5, 50.0, 62.5, 75.0, 87.5 and 100 s respectively and with the resolution N = M = 400.Here h and x have units of m.

2 Figure 10 :
Figure 10: Plot of (a) the normal mode frequencies ω r (m f /m v ) and (b) the normal mode decay rates ω i (m f /m v ) for the system parameters given (5.59).The black lines give the exact values from solving the characteristic equation (3.22) while the circles are results of the numerical simulations.Here m v and m f have units kg and ω r , ω i have units rad s −1 .

Figure 11 :
Figure11: Plot of the vessel displacement q(t) for the system parameters in (5.59) for (a) m f /m v = 1, (b) 0.5 and (c) 0.25.The black lines are the result of the numerical simulations, while the red dashed lines are the exact linear result for the mode labelled 1 in figure10(a,b) and for the mode labelled 2 figure10(c).Here t has units s and q has units m.

Figure 13 :
Figure 13: Contour plot of the numerical decay rate ω num i in the (Ω, T on )-plane for the system parameters (5.50) with L 1 = L 2 = 0.5 m and q = 5 × 10 −4 m and porosity function β 1 (τ ).The black contours gives the value ω i = ω max i which is the maximum linear decay rate for the system if the baffle has β = β max .Here Ω has units s −1 , T on has units s and ω i has units rad s −1 .

Figure 15 :
Figure 15: Plot of the numerical decay rate ω num i as a function of the time dependent porosity function amplitude A for the system parameter (5.50), L 1 = L 2 = 0.5 m and the porosity function β 2−5 (τ ) in (5.62) with (a) Ω = 5 s −1 and (b) Ω = 25 s −1 .The black dashed line gives the value ω i = ω max iwhich is the maximum linear decay rate for the system if the baffle has β = β max .Panels (c) and (d) give plots of the corresponding q(τ ) figures for cases (Ω, A) = (0, 0.508) -red dashed line, (5, 1.2) -blue dot-dashed line and (25, 2) -black line.Here A has units sm −1 , ω i has units rad s −1 , τ has units of s and q has units of m.

ΩAFigure 16 :
Figure 16: Contour plot of the numerical decay rate ω num i as a function of the time dependent porosity function parameters (Ω, A) for the system parameters (5.50), L 1 = L 2 = 0.5 m and T on = 2s.The black contours gives the value ω i = ω max i which is the maximum linear decay rate for the system if the baffle has β = β max .Here Ω has units s −1 , A has units sm −1 and ω i has units rad s −1 .

1 ρ 1 ρ
[h j U jt + h jt U j ] dx, [−h j (U j U jx + gh jx + q) − U j (h j U j ) x ] dx, ) which implicitly gives an equation for the unknown variable w at a = a B .Equation (4.32) essentially defines the position of the baffle a B in the Lagrangian framework.However, we find the form (4.32) troublesome to work with, as a B needs to be found via interpolation which introduces unnecessary errors.This can be overcome by forming an ODE for the value of a B (τ ) by differentiating(4.32)withrespect to τ to give ∂x ∂τ (a B , τ ) + ∂a B ∂τ ∂x ∂a (a B , τ ) = 0. Using (4.29) we can then write this equation as also show the presence of 2 j U j ) x ρ dx, (L 1 , t)U 2 (L 1 , t) U 2 (L 1 , t) 2 + q2 + gh 2 (L 1 , t) 2 U 2 (L 1 , t) (L 1 , t)U 1 (L 1 , t) U 1 (L 1 , t) 2 + q2 + gh 1 (L 1 , t) 2 U 1 (L 1 , t), (A-2)