An Improvement of Port-Hamiltonian Model of Fluid Sloshing Coupled by Structure Motion

: The ﬂuid–solid interaction is an interesting topic in numerous engineering applications. In this paper, the ﬂuid–solid interaction is considered in a vessel attached to the free tip of a cantilever beam. Governing coupled equations of the system include the Euler–Bernoulli equation for bending of a beam, torsion of a beam, 2-D motion of the rigid vessel, and rotating shallow water equation of ﬂuid sloshing in the vessel. As an essential portion in the numerical simulation of the vibration control of this ﬂuid–plate system is the accurate modeling of sloshing; the partial differential equations of the system are modiﬁed by approximation of velocity proﬁle. The suggested method is validated by experimental results of a piezoelectric actuated clamped rectangular plate holding a cylindrical vessel. These sloshing interactions with elastic test cases illustrate the mass conservative characteristics of the method as well as its stability in a prompt change of the vessel situations.


Introduction
The problem of simulating the interaction between fluids (water and air) and structure (tank wall) in partially filled storage vessels is discussed in many kinds of research as it is of interest to numerous industrial applications, for instance, marine science, aircraft design, seismic engineering, and fuel vessels in vehicles. In the study of fluid-structure interaction (FSI) in sloshing tanks, the numerical investigation is reasonably efficient and cheaper than the full-scale experiments [1]. Nicolici and Bilegan [2] studied the earthquake excited sloshing in flexible tanks by one-way coupling of finite element method (FEM) and computational fluid dynamics (CFD) techniques. They found that the wall elasticity augments the impulsive pressure.
As the excitation frequency of the system approaches the sloshing fundamental frequency, it can hurt the other structures in the system. Since the study of sloshing is more important in this field [3], time-independent smoothed particle hydrodynamics [4] and finite difference analysis [5] of viscous and inviscid fluid sloshing in a rectangular domain is performed by researchers. The impact of sloshing water on walls [3,6,7] is not negligible and can affect the coupled system response. The smoothed-particle hydrodynamics (SPH) computational method is a mesh-free Lagrangian method which considers the fluid as the collection of separate particles [8]. These particles affect each other over a smoothing length with a kernel function (such as cubic spline and Gaussian) [9].
The advantages of SPH over traditional grid-based methods are higher accuracy, guaranteed mass conservation, higher performance in large distortion conditions, pressure evaluation by neighboring particles, and easy free surface tracking [10]. On the other hand, the disadvantages of SPH are in the evaluation of free surface by polygonization method, lower performance in small distortion conditions, and numerous particles requirements [11,12]. Although the SPH method is used for solid mechanics computations [13][14][15], due to its instability in tensile case (stretching cause the tear) [15], first-type boundary conditions (Dirichlet) [16], and high computational cost [17], the solid part in this paper is done by finite element method. The common method in the literature for solid-liquid interaction problems is the use of SPH for the fluid part and FE for solid part [17][18][19][20]. In the works of Johnson [17,18] and Attaway [19], the SPH-FEM interface is controlled by FE motion while at each FE point the force is calculated from the nearest SPH particles. Since this method cannot be considered a conservative interface method, to have an interface preserving method with a higher computational cost method, Sonia Fernández-Méndez et al. [20] defined a buffer zone in FSI interface that is affected by both zones but does not allow the SPH and FEM to interrupt each other. Recently, the study of sloshing-plate interaction is performed in various studies [21][22][23]. In the study of Yu et al. [21], the plate is located within the sloshing system, while, in the studies of Kotyczka and Maschke [22] and Ribeiro et al. [23], the solid system is confined the sloshing. Ribeiro et al. proposed the FEM-based port-Hamiltonian formulation with the Saint-Venant approximation for Shallow Water Equations (SWE). Their port-Hamiltonian model of liquid sloshing in moving containers is applied to a fluid-structure system [23]. To control the liquid sloshing dynamics [24], Lasiecka and Tuffaha [25] used the boundary feedback control in fluid-structure interactions. This method is the solution to some automatic control problems for sloshing phenomena and water-tank systems [26,27]. Coron [28] investigated the local controllability of a 1-D tank containing a fluid modeled by applying the calculus of variations on the shallow water equations. In addition, the stabilization problem of a 1-D tank containing a fluid modeled by the shallow water equations was solved by Prieur and Halleux [29]. This method was used by Robu [30] to active vibration control of a fluid-plate system.
As it is important in the port-Hamiltonian systems theory, the modeling of each part plays a significant role in the numerical calculation of the whole system [31]. In the various research for a coupled fluid-plate system excited by piezoelectric actuators, the authors use the linearized shallow water equation for the fluid part in the water-tank system [30,[32][33][34]. Such simple models are also used for the control problems for water-tank systems [26], dynamic coupling between shallow-water sloshing and horizontal vehicle motion [35], sloshing in vessels undergoing prescribed rigid-body motion in three dimensions [36], and stabilization of a one-dimensional tank containing a fluid modeled by the shallow water equations [29]. The shallow water equation is also recently studied in the port-Hamiltonian framework but with a different goal: modeling and controlling the flow on open channel irrigation systems [37,38]. In this case, there is neither rotation nor translation of the channels. The effect of the nonlinear term has not been considered in the literature.
The aim of the current manuscript is to present a new dynamic simulating of the coupled system by driving the equation of motion from the port-Hamilton method and consider the nonlinear effects of fluid on the system by a comprehensive model. By the new modifications, the accuracy of the port-Hamiltonian formalism is improved. First, the state equation of the moving container system is derived and benchmarked by the frequency response of the coupled system. Second, the code is validated by experimental results and shows the improvement of the prediction of liquid sloshing response. New improvements offer a better model of sloshing which augment the accuracy of simulation in association with arbitrarily complex vibrational systems. Finally, the importance of each mode of energy is discussed in the coupled system.

Mathematical Modeling
The fluid, structure and rigid body considered in the current problem are plotted in Figure 1. The system shown in Figure 1 includes a cantilevered aluminum plate with a plastic vessel mounted by fixation rings near the free end. The actuator/sensor patches are used to control the vibration of the moving container system. Table 1 presents the geometrical parameters and physical properties of the system. Here, the solid plate is assumed as a bending and torsion element with Euler-Bernoulli beam model where the bending and torsional modes are operating individually. Subsystems are considered to affect each other by kinetic and kinematic interactions. The experimental system from Ribeiro et al. [23] was chosen to compare the numerical simulation with experimental results. Their experimental rig is designed for the application of airplane fuel tank and wing design where the light materials adopted for the wing caused low fundamental frequency as the order of sloshing frequency. Navier-Stokes fluid flow is modeled through the following mass and momentum conservation equations: where u, ..
x b , and g are vectors and presented in bold. The equation inside the boundary layer is in and for y-direction  Outside the bottom boundary layer, the potential flow is dominated with the zero velocities on the walls u(x = 0, y) = 0 (7) constant pressure at the free surface p(x, y = η) = 0 (10) And kinematic velocity at the free surface The Cartesian components of acceleration of the rotating system are: By the use of velocity potential: and integration of continuity equation versus y η y=−h 0 ∂u ∂x which leads to ∂η ∂t To derive a linearized form of momentum equation, first the pressure term should be evaluated. By y integration of the momentum equation in y-direction and differentiate with respect to x the first term in Equation (13) is simplified by the aid of Leibniz rule where the use of velocity potential the second term in the above equation is simplified as In addition, the second term in Equation (18) is simplified by the aid of Leibniz rule Finally, Equations (17) where the last term in Equation (22) could be rewritten by the use of 1 inside the boundary layer and outside the bottom boundary layer is Integrating the above equations from the bottom to free with surface respect to z-direction divided by the total length, ∂u(x,y=η) ∂t where δ is the length of the boundary layer. Replacing the vertical velocity in inertial term with (26) causes to the following simplification and the time variance of vertical velocity in Equation (24) with (inviscid limit The viscous term in Equation (21) is simplified by Stokes first problem. The viscous term in Equation (25) is rewritten as If the bottom wall velocity is u b e i(ωt+ϕ) and considering the viscous effect and neglecting the inertial effect, the governing boundary layer equation by If the semi-infinite fluid zone is considered, which has the solution of error function, there is no geometric length scale If we consider a bounded solution, the bottom shear is sinhy cos y + cosh y sin y (sinhy cos y) 2 + (cosh y sin y) 2 + i sinhy cos y − cosh y sin y (sinhy cos y) 2 + (cosh y sin y) 2 (33) where y = gω 2µ h. Equation (33) can be simplified to Finally, Equation (23) is rewritten as The friction is found with a linear relation of top velocity by the coefficient of where the shear at the boundary layer is ignored (τ −h 0 +δ = 0) and the shear at the bottom is replaced 2 factor because of π 4 phase difference between velocity and shear, and the factor (2 + 2h 0 b ) because of side effect. The governing equations presented in Equations (15) and (25) with the simplifications of Equations (26) where σ = tanh(kh 0 ) (16) and (37), the 1-D equation of motion is completed while the coupling the frequencies are a significant problem in design.
To model the dynamics of the system, it could be considered in two subsystems: a beam with torsion and bending mode of freedom ( Figure 2) and a fluid in a moving vessel (Figure 3).  The equations presented here could be derived from the classical Lagrangian method. The Lagrangian of vessel moving by external force and solid connected in base is where the system is described by the following quantities: a is the tank length; b is the tank width; ρ is the fluid density; g is the gravity; m v is the mass of the empty vessel; X(t) represents the horizontal position of the vessel (with respect to an inertial frame); x is the position along the vessel (measured in local coordinates from the middle of the vessel); h(x,t) is the height profile of fluid (−a/2 < x < a/2); and u(x,t) is the fluid speed (measured relative to the vessel). The assumptions are: • The fluid and the tank are viewed as an open system. • The fluid governing equations are obtained from Euler-Lagrange formulation.

•
There is no flow through the tank walls.

•
The effect of rigid-body inertias of the tank and the mass of the tank is considered.

•
The structure is considered a beam with two independent bending and torsion motions.

•
Because of symmetrical installation of the piezoelectric patches along the torsion axis, their effect on torsional motion is neglected.

•
The piezoelectric patches actuation affects only the bending motion of the beam.

•
Here, the linear Euler-Bernoulli beam with distributed piezoelectric actuators is used for the bending.

•
The electric field is constant through the piezoelectric patch thickness.

•
The voltage is uniform along the patch.

•
The beam is considered as a long beam under small displacements.

•
The fixed end for the torsion and for the bending is considered at the clamped point.

•
A weak formulation is used to overwhelm the derivation of the non-smoothness function.
By introducing the constraint of which is obtained from mass conservation, the following variational problem of finding the minimum Lagrangian (where constraints are adjoined by Lagrange multiplier) is obtained: To find the governing equations of the system, for any small variation of δh of h such that and then, thanks to an integration by parts, In addition, for any small variation of δu of u such that δu(x, 0) = δu(x, τ), and the boundary condition of δu −a 2 , t = δu a 2 , t = 0, Equation (36) yields and then, thanks to an integration by parts, Replacing the value of ∂λ ∂x by Equation (45) into Equation (43), A differentiation with respect to x and use of Equation (45) yields ∂u ∂t which is the momentum conservation equation. where m e f f = m v + z=L z=0 ρ s t s w s dz and k e f f = z=L z=0 E s I s X 2 ∂w ∂z 2 dz. As the system is controlled by a piezoelectric sheet, the following function is added to the bending function: The Hamiltonian of the parts of the system, equation of motion of the parts of the system, and boundary condition of the system are presented in Tables 2-4, respectively. Table 2. Hamiltonian of the parts of the system.

Sub-system. Type Formula
Beam bending Potential  Table 3. Governing equations of the parts of the system.

Sub-Formula
Beam bending Table 4. Boundary condition of the parts of the system.

Sub-Formula
Beam bending w For the node i = 1 to N b , the discretized bending equations are ∂ ∂t For the node i = 1 to N t , the discretized torsion equations are ∂ ∂t For the node i = 1 to N f , the discretized shallow water equations are ∂ ∂t The equations of motion of the rigid body of vessel are linear momentum in x-direction, angular momentum in z-direction, and angular momentum in y-direction: The total mass of the rigid body contains mass of the empty vessel (m e = ρπ·L out ·(D out 2 − D in 2 )/4), the tip mass (m tip ), the mass of two fixation rings (m ring = 2ρ ring ·π·t ring ·((D out ring ) 2 − (D in ring ) 2 )/4), and reduction of plate hole mass (m h = ρ b π·t b ·D out 2 /4). The tip mass (m tip ), located at the vessel tip, consists of plastic disc, glue, and a plug. Although it is not a simple geometry, here, the effects on the system mass and inertia are approximated as a point mass (m R = m e + m tip + m ring − m h ). The I r is total vessel mass moment of inertia around the z-direction, since vessel and fixation rings are symmetric around the z-axis. The total vessel mass moment of inertia is the sum of the empty vessel inertia (I e = m e (3/4(D out 2 − D in 2 ) + L out 2 )/12), the tip mass inertia (I tip = m tip L tip 2 ), the mass of two fixation rings (I ring = m ring /12(3/4 ((D out ring ) 2 + (D in ring ) 2 ) + (t ring ) 2 )), and reduction of plate hole mass (I h = m h ·D out 2 /16). The J R is total mass moment of inertia around the y-direction (J R = I r + I f ) Total vessel mass moment of inertia around the y-direction is the sum of vessel moment of inertia around the y-direction, which is similar to I r ·(J R = I e + I tip + I ring − I h ) and the frozen fluid mass moment of inertia I f = m f L eq 2 /12 Finally, the fluid-solid dynamic coupling equations are If the cylindrical vessel is filled to the height h (0 < h < D in ), then the total volume of the fluid is where the filling angle, α, is found from h = D in 2 1 − cos α 2 . The effective depth of the equivalent rectangular vessel is found from the depth of the fluid in the cylindrical vessel at rest and the effective length of the equivalent rectangular vessel is found from depth and equivalent volume (L eq = V f /W v ).

Numerical Results
Having implemented the FE in a MATLAB-based code, the resulting transient problem is then solved to determine the coupled response of illustrative beam motion, subjected to the various boundary and load conditions. In this section, first, the code is validated by comparison of the frequency response of the coupled system with experimental results of reference [23]; and, second, the state equation of the system is solved and the dynamic of the variables of the system is shown versus time. The most important parameters are the sloshing phenomenon (regarding evolutions of the free surface) in a moving tank and displacement and rotation of the structure and rigid body which is numerically investigated. By varying the time, interesting characteristics, such as sub-system energies, and dynamic responses of the structures in both time and frequency domains are presented.
To choose basis functions ϕ i for finite element method with the property of being one at the point i and zero at the other points, the simple linear functions are used. As shown in Figure 4, the simplest choice of course linear functions (blue lines) caused the values of basis functions equal to the maximum value (here is one) at the desired point and zero other grid nodes (plus symbols). During the semi-discretization of the equations, a weak formulation is used to overcome the problem of discontinuous derivate which made some difficulties in existence and uniqueness results for such systems. Finite-element method discretization of the second order derivate in Equation (50) using spline approximation leads to the matrix equation as where The FE method makes it possible to effectively simulate three-dimensional FSI problems in which a vessel is under arbitrary rotation, translation, and deformation in the fluid domain. In Table 5, the comparison of the first natural frequency of the coupled system obtained from generalized eigenvalues with those of experimental results of Reference [23] is presented. The first natural frequency of fluid dynamics introduces modes that are symmetric with respect to the center of the tank. Calculated natural frequencies obtained from the semi-discretization model computed for different values of 40 basis functions and comparison with the experimental results. As shown for the same filling ratios, the modified port-Hamiltonian method presents a lower error in comparison with the previous numerical method. Notice that a quite good agreement with experimental results was obtained. For the 25% filling ratio, first modes agree with an error less than 6%. For the 50% filling ratio, first modes agree with an error less than 6%. Larger errors appear for larger filling ratios in [23]. One of the reasons is that the fluid equations used in this paper assume the shallow water hypothesis, which is more accurate for small filling ratios of the tank. The sets of governing equations are rearranged in state space formulation and solved by Runge-Kutta time integration in MATLAB with ∆t = 10 −4 s. Figure 5 illustrates the tip displacement versus time. Initial condition considered here is the stationary case, with four Newtons tip force. The figure shows the maximum downward deflection of the beam is 3 mm. Figure 5 also shows that in steady conditions beam tip deflection started from 3 mm with two small oscillations in each period. The classical laminated plate theory gives the relationships among the plate bending moments, torsional moment and curvatures. The dimensionless coupling term ψ = K/ √ EIGJ here is small, where EI, GJ and K are bending stiffness about the y-axis, torsion stiffness and bending-torsion coupling stiffness, respectively. Different modes of operation lead to a wide range of frequencies, which are partially visible for example in the fluid forces acting on the structure. Angular deflection versus time graph for a vibrating cantilever beam subjected to the initial conditions can show the effect of bending motion on rotational mode. The tip torsion versus time is shown in Figure 6. As shown, the amplitude of orientation oscillation orientation of the sloshing vessel becomes larger with time. This is because the cylindrical shape of the vessel reduces resistance from the fluid flow so that the increasing deformation due to the time evolution leads to the reduction of the resistance in the vessel motion. The figure shows steady behavior started from equilibrium with two small oscillations in each period. In steady conditions, beam tip torsion started from zero degrees with two small oscillations in each period where the maximum reaches to 0.0003 radians.   Figure 7 demonstrates the outcome after 10 s of recreation. As shown, the direct and nonlinear reenactments correspond and the state of the waves is similar to the primary modular bending of the liquid. Figure 8 presents the nonlinear horizontal momentum of the fluid versus time. The fluid starts in still condition and is stimulated using symphonious voltages for the piezoelectric patches, with a repeat close to the principle trademark repeat of sloshing. The maximum value of momentum fluctuations is 0.006 duration of motions on the sloshing phenomena in water tanks. Figure 9 reveals time dependency of the fluid angular momentum. The role of the rotational mode is stabilizing the coupling between fluid and structure. The smoothness of the liquid nonlinear conduct influences the basic mode reaction, lessening the amplitude of the vibrations. The vessel is moved by the piezoelectric wafers and the water reacts to it. The torsion mode was made a degree of freedom to release the energy of the system in the rotational degree of freedom and help maintain a stable coupling between fluid and structure.     Figure 5, it demonstrates the time reaction of the tip speed of the bar for these same reenactments. From the momentum equation in the inertial reference system, this phase lag evolves with slow time dynamics. Increasing torsional stiffness increases lateral stability. On the other hand, similar to the horizontal displacement, the time histories for rotation were obtained. Figure 11 displays rigid-body torsion versus time. In the sloshing-structure simulation shown in Figure 11, same as Figure 6, the level of torsion is low in comparison with the bending mode. The vessel is moved by the piezoelectric wafers and the water reacts to it. The torsion mode was made a degree of freedom to release the energy of the system in the rotational degree of freedom. The values calculated by Euler-Bernoulli theory completely prove that for the current study the bending can be considered to act independently in this case, well higher than that computed by the coupled Euler-Saint Venant bending-torsion theory. Although the Euler-Bernoulli theory would predict correctly the dynamics of the beam in Figure 11, the bending-torsion coupling effects cannot be neglected when applying random loads.   Simulating such a complex phenomenon, related to bending rotations of the beam, reveals the higher level of the bending in comparison with the torsion and smoothness of this motion because of the sloshing coupling. This quasi-periodicity can also be confirmed by looking at the predicted motion depicted in Figure 12. This effort sought to measure the dynamic nonlinear bending response of a cantilever beam. That shows the Euler theory for bending and Saint Venant theory for torsion are sufficient for current analysis. Note that, in most of the three-dimensional numerical datasets (see , the periodic motion is predicted. Figure 13 portrays percent change of the Hamiltonian of the system versus time. A small decrease of the total Hamiltonian of the system shows the general conservative for of the current system with low damping ratio. Notice that a very decent concurrence with test results was obtained. For the 75% filling proportion, the Hamiltonian of the system decrease by time. Figure 14 exposes the variation of the Hamiltonian of the sub-systems versus time. As shown, most energy is distributed between the fluid and bending modes with a 90-degree phase difference. The next order is the rigid-body motion. The lowest energies relate to the torsion of the system and rigid-fluid coupling.   A result of the fully coupled system exemplifies in Figure 15 essential response frequencies of the system. The presented eigenvalues of the coupled system were obtained for filling ratio of 75 percent. The initial five modes are basically because of the coupling between the sloshing elements and the principal bending method of the structure. The sixth mode is commanded by the torsion elements. The seventh and eighth modes are bending modes. The liquid element additionally presents modes that are symmetric regarding the focal point of the tank. As can be seen in Figure 15, the first five normal modes are pure sloshing modes and the sixth normal mode is fundamental torsional mode when the effect of bending-torsion coupling is ignored. If the effect of bending-torsion coupling is considered, there is no strong coupling between bending and torsional displacements for the first six modes, since the sixth mode remains the fundamental torsional mode. Figure 15 shows that, although the bending-torsion coupling does not significantly affect the natural frequencies of sloshing, there can be appreciable differences in the mode shapes.

Conclusion
An improvement of the port-Hamiltonian model of fluid sloshing coupled by structure motion was performed in this study. A new method to simulate the nonlinear effects of sloshing in confined vessels is presented. A set of governing equations including the shallow water equations in the onedimensional case and bending and torsion of a beam are solved by Runge-Kutta time integration. This paper offers a modified shallow water equation which leads to a proper distribution in each computational cell, which prevents numerical errors in frequency response analysis of the system and maintaining continuity condition. The essential enthusiasm of utilizing this formalism, with regards to this work, is that it gives an efficient, secluded approach for demonstrating complex multimaterial science marvels. Especially in contrast to past work utilizing shallow water conditions for the reenactment of sloshing in moving compartments, this paper expressly characterizes information and yields association ports that can be utilized for coupling components of an unpredictable framework in a simple and methodical way. The main contributions of this paper are the following:

Conclusions
An improvement of the port-Hamiltonian model of fluid sloshing coupled by structure motion was performed in this study. A new method to simulate the nonlinear effects of sloshing in confined vessels is presented. A set of governing equations including the shallow water equations in the one-dimensional case and bending and torsion of a beam are solved by Runge-Kutta time integration. This paper offers a modified shallow water equation which leads to a proper distribution in each computational cell, which prevents numerical errors in frequency response analysis of the system and maintaining continuity condition. The essential enthusiasm of utilizing this formalism, with regards to this work, is that it gives an efficient, secluded approach for demonstrating complex multi-material science marvels. Especially in contrast to past work utilizing shallow water conditions for the reenactment of sloshing in moving compartments, this paper expressly characterizes information and yields association ports that can be utilized for coupling components of an unpredictable framework in a simple and methodical way. The main contributions of this paper are the following:

1.
A new model for rotating-translating shallow water equations in moving containers is offered within the port-Hamiltonian equations. 2.
The suggested model was used to solve a fluid-structure interaction problem, which includes a piezoelectric actuator on a flexible beam (bending element and torsional element), and a rigid vessel. Each subsystem was considered within the port-Hamiltonian equations and coupled using the interaction ports.