Research on liquid sloshing model of partially-filled tank by nonlinear external excitation

1. Introduction

Truck crashes are particular importance because the size and weight of these vehicles results in a greater likelihood of fatalities when they are involved in a crash compared to passenger cars. Truck accidents in the state of New Jersey represent about 8 percent of all accidents. In 1998, 1999, and 2000, there were 20,025, 21,561, and 28,328 accidents respectively, involving trucks in the State. These accidents represent all fatal, injury and property damage accidents for all types of trucks, with the majority of these accidents property damage accidents [1]. And a US study has reported that the average annual number of cargo tank rollovers is about 1265, which takes up 36.2 % in the total number of heavy vehicle highway accidents [2]. Although there are many reasons that lead to tank vehicle rollover accidents, such as driver’s fatigue, overtaking, bad road and weather conditions, and so forth, liquid sloshing in the tank are the main factors [3, 4].

Liquid sloshing inside partially filled containers has been investigated extensively in the aerospace field in the 1950’s and 1960’s. The purpose was to study the effect of partially filled fuel tank of different shapes on the stability of rockets and airplanes. Several mathematical models as well as equivalent mechanical models have been developed to simulate liquid sloshing effects in partially filled containers. The use of simple and compound pendulums has been proved successful in simulating liquid-sloshing effects in partially filled containers of various shapes, especially spherical and cylindrical.

The roll stability of a partially filled tank is influenced by both the centre of gravity (C.G.) height and the magnitude of lateral load transfer. Differential cross sections have been analysed to optimize the roll stability of the tank. The circular cross section has a high centre of liquid location, but considerably less lateral load transfer with low filled volumes under a steady turning lateral acceleration field. On the other hand, modified elliptical tank compared to circular tank yield lower C.G. Height and relatively larger lateral load transfer with low filled volumes. Therefore, the modified elliptical tank represents higher roll stability limits than a circular tank [5].

Sloshing forces depend on the freedom that the liquid has to move into the tank, and on the amount of liquid participating in the motion. While such freedom of movement of the liquid is characterized by the free surface dimensions of the liquid, the amount of moving liquid depends on the filled level of the tank. Because the tank used for carrying liquids have in general a non-rectangular cross section, e.g. elliptic, modified elliptic and circular, the filled level will define the amount of moving liquid and the size of the free surface of the liquid [6].

Many studies have been carried out on liquid flow and sloshing characteristics that happened in the tank. The main methods are the quasi-static (QS) method and the experimental method. The cargo’s static moment at a specified point on a tank vehicle can be approximated by calculating the transient C.G. of the liquid in the tank. Then, the liquid sloshing effect on tank walls can be analysed. It is convenient and simple to obtain liquid sloshing force using the QS method. However, the analysis results have poor accuracy [7-9]. The experimental method, by building a test platform or using test tank vehicles, liquid sloshing phenomenon can be observed and relevant parameters can be monitored by reproducing liquid sloshing [10-11].

In a word, QS method is simple and intuitive and can be used for approximate calculation, but the accuracy is not high. The equivalent mechanical model in the impact of liquid in the tank was more accurate than the QS method. When the liquid impact is large amplitude, the study should be nonlinear coupled with liquid shock. Thus, the aim of this study is to 1) establish liquid free surface equation of motion, 2) set up model of liquid sloshing force and moment caused by nonlinear excitation, two models coupling to take simulate result, could be real reaction tank liquid sloshing situation.

2. The theoretical model

2.1. The free-surface equilibrium formulation of the liquid within a tank

The elliptical tanks overturning moment for liquid goods is caused by lateral and vertical displacement, those are the potential influencing factors of the vehicle roll. Steady state shifting of the liquid cargo C.G. has been modelled by different approaches, including geometrical modelling and experimental methods [11]. We would like to solve the equation of the elliptical cross section mass centre of the tank. Elliptical tank under a steady lateral acceleration field as shown in Fig. 1. The Angle is the $\phi$ between the horizontal direction and the free surface, changed the range is 0 to 90°, that is above 30° belongs to large amplitude.

The third and fifth order vibration modes are antisymmetric wave. Under this three order mode, liquid vibration can lead to C.G. of the lateral displacement. The fourth and sixth order vibration modes are symmetric wave. Under these three modes, liquid vibration will not be able to make C.G. lateral displacement. In addition, the vast majority of liquid C.G. happens larger displacement under the sloshing of first order modal. But the third order and fifth order modal, the size of the maximal displacement and displacement of the liquid C.G. number is greatly reduced. Therefore, under the sloshing of first order modal, liquid C.G. displacement is the largest. In a word, the first-order sloshing model is the lateral sloshing of the liquid in the tank is the most important modal. Furthermore, the sloshing effects can be described by liquid C.G. displacement. The first order model of the liquid free surface can be described by the inclined plane to approximate.

The $x$ and $y$ position of the liquid cargo C.G., as shown in Fig. 1.

$y_h$ is defined as:

1

$y_h=a_t\times \sqrt{1-{(2f-1)}^2}.$

The slope of the liquid along the corresponding ellipse, is:

2

$\frac{b_t}{a_t}=\frac{b_F}{a_F}=\frac{b_{CG}}{a_{CG}}=r.$

The C.G. position of the liquid depends on the free surface factors such as slope and cross section shape of tank, ignore the mass center of liquid along the longitudinal axis of the tank on the small changes that the cross section of the C.G. position of the same uniform distribution of liquid in the tank, there is no density difference.

The liquid C.G. position is given by:

1) $h\in \left(0,b_t\right]$, the filled level is 0 to 50 %, and $f=h/2b_t$, the slope of the C.G. ellipse:

3

$A_F=3a_t{y}_h\sqrt{1-{\left(\frac{y_h}{b_t}\right)}^2}+a_t{b}_t\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\sqrt{1-{\left(\frac{y_h}{b_t}\right)}^2}\right),$

4

$b_{cg0}\left(h\right)=-\frac{2b_t[1-(2f-1{)}^2{]}^{3/2}}{A_F}.$

2) $h\in \left(b_t,2b_t\right]$, the filled level is 50 % to 100 %, $f=h/2b_t$, the slope of the C.G. ellipse:

5

$A_F=a_t\left(y_h\sqrt{1-{\left(\frac{y_h}{b_t}\right)}^2}+b_t\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{y_h}{b_t}\right)+\frac{\pi }{2}\right)\right),$

6

$b_{cg0}\left(h\right)=-\frac{2a_t{b}_t[1-(2f-1{)}^2{]}^{3/2}}{3A_F}.$

The liquid C.G. position is given by:

7

$y\left(h\right)=\left(a_t-b_{cg0}\left(h\right)\frac{a_t}{b_t}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(-\frac{r}{g_y}\right),$

8

$z\left(h\right)=b_t-\left(b_t-b_{cg0}\left(h\right)\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(-\frac{r}{g_y}\right).$

Fig. 1Elliptical tank under a steady lateral acceleration field

Assumptions made include neglecting of sloshing damping. However, the most important simplification in the model is the uncoupled response of the liquid cargo with the truck. In that respect, it could be expected that the coupled sloshing frequency could be lower than the calculated, due to truck suspensions deformation. Hence, liquid sloshing frequency and is much less than the natural tanker vibration frequency. The resonance phenomenon can be ignored.

Due to $h\in \left(b_t,2b_t\right]$, the filled level f is in the 50 % to 100 %, and the roll angle 5° remains the same, the elliptic diameters is $a=$2. 2 m and $b=$1. 6 m, the circular diameter is $r=$0.94 m. The coordinates of the $z$ axis and $y$ axis change with the lateral acceleration ($a_y$) are shown in Fig. 2. Liquid C.G. height and lateral displacement in different filled level change with the lateral acceleration ($a_y$) and the liquid in the tank with ay all show a trend of rising, but higher filled level liquid mass centre height range is smaller, smoother. Compared with the cross section area of the circular tank, the elliptical cross section of tank centre of mass is lower. As Fig. 2 shown, circular cross section and the elliptical cross section in the equal area, the circular cross section of the $z$ axis coordinate values have been maintained a higher height of mass centre. As the condition of an equivalent filled level, an elliptical cross section of tank car has better roll stability than the circular.

Fig. 2The values of ay/z, y/z

a) The values of $a_y/z$

b) The values of $y/z$

2.2. The liquid sloshing model of partially-filled tank

When analysing the motion of the tank in 6-DOF, it is convenient to define two coordinate frames as indicated in Fig. 1. The moving coordinate frame $X_0{Y}_0{Z}_0$ is conveniently fixed to the vehicle and is called the body-fixed reference frame. The origin $O$ of the body-fixed frame is usually chosen to coincide with the C.G. When C.G. is in the principal plane of symmetry or at any other convenient point if this is not the case. For tank, the body axes $X_0$, $Y_0$ and $Z_0$ coincide with the principal axes of inertia, and are usually defined as:

$X_0$ – longitudinal axis (directed from aft to fore), $Y_0$ – transverse axis (directed to board), $Z_0$ –normal axis (directed from top to bottom).

Table 1Notation used for the tank

DOF		Forces and moments	Linear and angular velocity	Position and Euler angles
1	$X$-direction (surge)	$X$	$U$	$x$
2	$Y$-direction (sway)	$Y$	$V$	$y$
3	$Z$-direction(heave)	$Z$	$W$	$z$
4	Rotation about $X$-axis (roll)	$K$	$p$	$\varphi$
5	Rotation about $Y$-axis (pitch)	$M$	$q$	$\theta$
6	Rotation about $Z$-axis (yaw)	$N$	$r$	$\psi$

The motion of the body-fixed frame is described relative to an inertial reference frame. For tank, it is usually assumed that the accelerations of a point on the surface of the Earth can be neglected. Indeed, this is a good approximation since the motion of the Earth hardly affects slow speed tank. As a result of this, an earth-fixed reference frame $XYZ$ can be considered to be inertial. This suggests that the position and orientation of the vehicle should be described relative to the inertial reference frame while the linear and angular velocities of the vehicle should be expressed in the body-fixed coordinate system. The value range is from 0° to 90° about $\varphi$, $\theta$ and $\psi$. When their value is above 30, this situation is large amplitude:

9

$\eta =[{\eta }_1^T,{\eta }_2^T{]}^T,\mathrm{}\mathrm{}\mathrm{}{\eta }_1=[x,y,z{]}^T,\mathrm{}\mathrm{}\mathrm{}{\eta }_2=[\varphi ,\theta ,\psi {]}^T,$
$\nu =[{\nu }_1^T,{\nu }_2^T{]}^T,\mathrm{}\mathrm{}\mathrm{}{\nu }_1=[u,v,w{]}^T,\mathrm{}\mathrm{}\mathrm{}{\nu }_2=[p,q,r{]}^T,$
$\tau =[{\tau }_1^T,{\tau }_2^T{]}^T,\mathrm{}\mathrm{}\mathrm{}{\tau }_1=[X,Y,Z{]}^T,\mathrm{}\mathrm{}\mathrm{}{\tau }_2=[K,M,N{]}^T,$

where $\eta$ denotes the position and orientation vector with coordinates in the earth-fixed frame, $\nu$ denotes the linear and angular velocity vector with coordinates in the body-fixed frame and $\tau$ is used to describe the forces and moments acting on the vehicle in the body-fixed frame.

Fig. 3Body-fixed and Earth-fixed reference frames

2.2.1. Euler anglers and Euler parameter transformation

The vehicle’s flight path relative to the earth-fixed coordinate system is given by a velocity transformation:

10

${\dot{\eta }}_1=J_1\left({\eta }_2\right){\nu }_1.$

Let a be a vector fixed in $A$ and $b$ is a vector fixed in $B$. Hence, the vector $b$ can be expressed in terms of the vector $a$, the unit vector $\lambda ={\left[{\lambda }_1,{\lambda }_2,{\lambda }_3\right]}^T$ parallel to $L$ (the axis of rotation) which $B$ is rotated about and the angle frame $B$ is rotated. The rotation is described by (see Hughes (1986) or Kane, Likins and Levinson (1983)) [19, 20]:

11

$b=\mathrm{c}\mathrm{o}\mathrm{s}\beta \alpha +\left(1-\mathrm{c}\mathrm{o}\mathrm{s}\beta \right)\lambda {\lambda }^T\alpha -\mathrm{s}\mathrm{i}\mathrm{n}\beta \times \alpha .$

Consequently, the rotation sequence from $a$ to $b$ can be written as:

12

$b=Ca,$

where $C$ can be interpreted as a rotation matrix which simply is an operator taking a fixed vector and rotating it to a new vector $Ca$. From Section 2.3 we obtain the following expression for $C$:

13

$C=\mathrm{c}\mathrm{o}\mathrm{s}\beta I+\left(1-\mathrm{c}\mathrm{o}\mathrm{s}\beta \right)\lambda {\lambda }^T-\mathrm{s}\mathrm{i}\mathrm{n}\beta S\left(\lambda \right).$

The principal rotation matrices can be obtained by setting ${\lambda }_1={\left[\mathrm{1,0},0\right]}^T$, ${\lambda }_2={\left[\mathrm{0,1},0\right]}^T$, ${\lambda }_3={\left[\mathrm{0,0},1\right]}^T$, respectively, in the general formula for $C$. This yields the following transformation matrices:

14

$C_{x,\varphi }=\left(\begin{array}{lll}1& 0& 0\\ 0& \mathrm{c}\mathrm{o}\mathrm{s}\varphi & \mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ 0& -\mathrm{s}\mathrm{i}\mathrm{n}\varphi & \mathrm{c}\mathrm{o}\mathrm{s}\varphi \end{array}\right),\mathrm{}\mathrm{}\mathrm{}{C}_{y,\theta }=\left(\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\theta & 0& -\mathrm{s}\mathrm{i}\mathrm{n}\theta \\ 0& 1& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}\theta & 0& \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right),\mathrm{}\mathrm{}\mathrm{}{C}_{z,\psi }=\left(\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\psi & \mathrm{s}\mathrm{i}\mathrm{n}\psi & 0\\ -\mathrm{s}\mathrm{i}\mathrm{n}\psi & \mathrm{c}\mathrm{o}\mathrm{s}\psi & 0\\ 0& 0& 1\end{array}\right).$

The notation $C_{i,a}$ denotes a rotation angle a about the $i$-axis. Notice that all $C_{i,a}$ satisfy the following property: $C^{T}C=C^{T}C=I$, $\mathrm{d}\mathrm{e}\mathrm{t}C=$1.

2.2.2. Linear velocity and angler velocity transformation

Let $X_3{Y}_3{Z}_3$ be the coordinate system obtained by translating the earth-fixed coordinate system $XYZ$ parallel to itself until its origin coincides with the origin of the body-fixed coordinate system. Then, the coordinate system $X_3{Y}_3{Z}_3$ is rotated a yaw angle about the $Z_3$ axis. This yields the coordinate system $X_2{Y}_2{Z}_2$. The coordinate system $X_2{Y}_2{Z}_2$ is rotated a pitch angle $\theta$ about the $Y_2$ axis. This yields the coordinate system $X_1{Y}_1{Z}_1$. Finally, the coordinate system $X_1{Y}_1{Z}_1$ is rotated bank or roll angle $\varphi$ about the $X_1$ axis. This yields the body-fixed coordinate system $X_0{Y}_0{Z}_0$, shown on Fig. 2. Fig. 4 is shown the rotation sequence according to the $xyz$-convention showing both the linear ($u,v,w$) and angular ($p,q,r$) velocities. The rotation sequence is written as:

15

$J_1\left({\eta }_2\right)=C_{z,\psi }^T{C}_{y,\theta }^T{C}_{x,\varphi }^T.$

From liquid in $\Omega$ satisfy the continuity equation, meet in $\mathrm{\Sigma }$ impermeable conditions, meeting the conditions of kinematics and dynamics in the $S$, available control equations as follows:

16

$\left\{\begin{array}{ll}{\nabla }^2\phi =0,& \text{into}\mathrm{}\mathrm{\Omega },\\ \frac{d\phi }{dn}=0,& \text{above}\mathrm{\Sigma },\\ \frac{d\phi }{dn}=\frac{dH}{dt},\mathrm{}\mathrm{}\mathrm{}\frac{d\phi }{dt}+\frac{dU}{dt}\cdot r+gH,& \text{above}S,\end{array}\right.$

17

$\left\{\begin{array}{ll}{\nabla }^2\gamma =0,& \text{into}\mathrm{}\mathrm{\Omega },\\ \frac{d\gamma }{dn}=0,& \text{above}\mathrm{\Sigma },\\ \frac{d\gamma }{dn}=L,L=\frac{{\nu }^2}{g}\gamma ,& \text{above}S.\end{array}\right.$

In above type, $H\left(x,y,z,t\right)$ for free surface wave height function; $n$ for the outward normal direction of wet surface of tank; $r$ for liquid particle diameter bolts relative point $o$. When lateral excitation of tank is $U\left(t\right)=$ 0, liquid is freely sloshing. So on the type of boundary problem is transformed into a type Eq. (17). In above type, $\gamma$ for liquid sloshing state function; $L$ for the wave height state function; $\nu$ for the freedom liquid sloshing characteristic frequency, Hz.

In that case, the kinematic equations can be described by two Euler angle representations with different singularities. Another possibility is to use a quaternion representation. This is the topic of the next section. Summarizing the results from this section, the kinematic equations can be expressed in vector form as:

18

$\left(\begin{array}{l}{\dot{\eta }}_1\\ {\dot{\eta }}_2\end{array}\right)=\left(\begin{array}{cc}{J}_1\left({\eta }_2\right)& 0_{3\times 3}\\ 0_{3\times 3}& J_2\left({\eta }_2\right)\end{array}\right)\left(\begin{array}{l}{\nu }_1\\ {\nu }_2\end{array}\right)\iff \dot{\eta }=J\left(\eta \right)\nu .$

The following coordinate transformation matrix for the Euler parameters:

19

$C=\left({\eta }^2-{\epsilon }^T\epsilon \right)I+2{\epsilon }^T\epsilon -2\eta S\left(\epsilon \right).$

The transformation relating the linear velocity vector in the inertial reference frame to the velocity in the body-fixed reference frame can be expressed as:

20

${\dot{\eta }}_1=E_1\left(e\right){\nu }_1,$

where $E_1=C^T$ with $C$ defined is the rotation matrix. Hence:

21

$E_1\left(e\right)=\left(\begin{array}{ccc}1-2\left({\epsilon }_2^2+{\epsilon }_3^2\right)& 2\left({\epsilon }_1{\epsilon }_2-{\epsilon }_3\eta \right)& 2\left({\epsilon }_1{\epsilon }_2+{\epsilon }_2\eta \right)\\ 2\left({\epsilon }_1{\epsilon }_2+{\epsilon }_3\eta \right)& 1-2\left({\epsilon }_1^2+{\epsilon }_3^2\right)& 2\left({\epsilon }_2{\epsilon }_3-{\epsilon }_1\eta \right)\\ 2\left({\epsilon }_1{\epsilon }_3-{\epsilon }_2\eta \right)& 2\left({\epsilon }_2{\epsilon }_3+{\epsilon }_1\eta \right)& 1-2\left({\epsilon }_1^2+{\epsilon }_2^2\right)\end{array}\right).$

Fig. 4The rotation sequence according to the xyz-convention showing both the linear (u,v,w) and angular (p,q,r) velocities, boundary conditions. Ω is liquid domain, Σ is wet surface of tank, S is free surface

a) Rotation angle $\psi$ about $Z_3$ ($w_2=w_3$)

b) Rotation angle $\theta$ about $Y_2$ ($v_2=v_1$)

c) Rotation angle $\phi$ about $X_1$ ($u_2=u_1$)

d) The boundary conditions

The angular velocity transformation can be derived by differentiating:

22

$\dot{C}\left(t\right)C^T\left(t\right)+C\left(t\right){\dot{C}}^T\left(t\right)=0.$

Let us define the matrix $S\left(t\right)=C\left(t\right){\dot{C}}^T\left(t\right)$ as:

23

$S\left(\omega \left(t\right)\right)=\left(\begin{array}{ccc}0& -{\omega }_3\left(t\right)& {\omega }_2\left(t\right)\\ {\omega }_3\left(t\right)& 0& -{\omega }_1\left(t\right)\\ -{\omega }_2\left(t\right)& {\omega }_1\left(t\right)& 0\end{array}\right),$

where $\omega \left(t\right)$ is a unique vector defined as the angular velocity of the body-fixed rotating frame with respect to the earth-fixed frame at time $t$. Introducing the notation $\omega \left(t\right)={\left[p,q,r\right]}^T$:

24

$\left(\begin{array}{lll}{\dot{C}}_{11}& {\dot{C}}_{12}& {\dot{C}}_{13}\\ {\dot{C}}_{21}& {\dot{C}}_{22}& {\dot{C}}_{23}\\ {\dot{C}}_{31}& {\dot{C}}_{32}& {\dot{C}}_{33}\end{array}\right)=\left(\begin{array}{ccc}r{C}_{21}-q{C}_{31}& r{C}_{22}-q{C}_{32}& r{C}_{23}-q{C}_{33}\\ p{C}_{31}-r{C}_{11}& p{C}_{32}-q{C}_{12}& p{C}_{33}-q{C}_{13}\\ q{C}_{11}-p{C}_{21}& q{C}_{12}-p{C}_{22}& q{C}_{13}-p{C}_{23}\end{array}\right).$

2.2.3. Newtonian and Lagrangian mechanics

The 6-DOF nonlinear dynamic equations of motion can be conveniently expressed as:

25

$M\dot{\nu }+C\left(\nu \right)\nu +D\left(\nu \right)\nu +g\left(\eta \right)=\tau ,$

where: $M$ – inertia matrix (including added mass); $C\left(v\right)$ – matrix of Coriolis and centripetal terms (including added mass); $D\left(v\right)$ – damping matrix; $g\left(\eta \right)$ – vector of gravitational forces and moments; $\tau$ – vector of control inputs.

Consider a body-fixed coordinate system $X_0{Y}_0{Z}_0$ rotating with an angular velocity $\omega ={\left[{\omega }_1,{\omega }_2,{\omega }_3\right]}^T$ about an earth-fixed coordinate system $XYZ$. Arbitrary body-fixed coordinate system $X_0{Y}_0{Z}_0$ with origin $O$ in the body-fixed frame is defined as:

26

$I_0=\left[\begin{array}{ccc}{I}_x& -I_{xy}& -I_{xz}\\ -I_{yx}& -I_y& -I_{yz}\\ -I_{zx}& -I_{zy}& I_z\end{array}\right],I_0=I_0^T>0,$

where $I_x$, $I_y$ and $I_z$ are the moments of inertia about the $X_o$, $Y_o$ and $Z_o$-axes and $I_{xy}=I_{yx}$, $I_{xz}=I_{zx}$ and $I_{yz}=I_{zy}$.

${\rho }_A$ as the mass density of the body. Consequently, we can represent the inertia tensor $I_0$ in vectorial form as:

27

$I_0\omega ={\int }_{V}r\times (\omega \times r){\rho }_{A}dV.$

Another useful definition of $I_0$ is:

28

$I_0=-{\int }_{V}S\left(r\right)S\left(r\right){\rho }_{A}dV={\int }_V\left(r^{T}rI-r{r}^T\right){\rho }_{A}dV.$

As Parallel Axes Theorem, the inertia tensor about an arbitrary origin $O$ is defined as:

29

$I_0=I_C-mS\left(r_G\right)S\left(r_G\right)=I_C-m\left(r_G{r}_G^T-r_G^T{r}_G{I}_{3\times 3}\right).$

$I_0=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{I_x,I_y,I_z\right\}$, yields the following set of equations:

30

$I_x=\frac{1}{12}m{l}^2+\frac{1}{3}m{b}_t^2\left[\left(1+\frac{1}{k^2}\right)-\left(\frac{1}{k^2}-1\right)\frac{F\left(k\right)}{E\left(k\right)}\right],$
$I_y=\frac{1}{12}m{l}^2+\frac{1}{3}m{a}_t^2\left[\left(2-\frac{1}{k^2}\right)+\left(\frac{1}{k^2}-1\right)\frac{F\left(k\right)}{E\left(k\right)}\right],\mathrm{}\mathrm{}\mathrm{}{I}_z=\frac{1}{3}m\left[a_t^2+b_t^2+b_t^2\frac{F\left(k\right)}{E\left(k\right)}\right],$

where:

31

$k=\sqrt{1-\frac{b_t^2}{a_t^2}},\mathrm{}\mathrm{}\mathrm{}E\left(k\right)={\int }_0^{\frac{\pi }{2}}\sqrt{1-k^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\phi }d\phi ,\mathrm{}\mathrm{}\mathrm{}F\left(k\right)={\int }_0^{\frac{\pi }{2}}\frac{d\phi }{\sqrt{1-k^2\mathrm{s}\mathrm{i}{\mathrm{n}}^2\phi }}.$

Hence 6-DOF equation of motion:

32

$X=m\left[\dot{u}-vr+wq-x_G\left(q^2+r^2\right)+y_G\left(pq-\dot{r}\right)+z_G\left(pr+\dot{q}\right)\right],$
$Y=m\left[\dot{v}-wp+ur-y_G\left(r^2+p^2\right)+z_G\left(qr-\dot{p}\right)+x_G\left(qp+\dot{r}\right)\right],$
$Z=m\left[\dot{w}-uq+vp-z_G\left(p^2+q^2\right)+x_G\left(rp-\dot{q}\right)+y_G\left(rq+\dot{p}\right)\right],$
$K=I_x\dot{p}+\left(I_x-I_y\right)qr+m\left[y_G\left(\dot{v}-uq+vp\right)-z_G\left(\dot{v}-wp+ur\right)\right],$
$M=I_y\dot{q}+\left(I_x-I_z\right)rp+m\left[z_G\left(\dot{u}-vr+wq\right)-x_G\left(\dot{w}-uq+ur\right)\right],$
$N=I_z\dot{r}+\left(I_y-I_x\right)pq+m\left[x_G\left(\dot{v}-wp+ur\right)-y_G\left(\dot{u}-vr+wq\right)\right].$

The vectorial representation of the 6-DOF equations of motion:

33

$M_{RB}\dot{\nu }+C_{RB}\left(\nu \right)\nu ={\tau }_{RB},$

where $\upsilon ={\left[u,v,w,p,q,r\right]}^T$ is the body-fixed linear and angular velocity vector and ${\tau }_{RB}={\left[X,Y,Z,K,M,N\right]}^T$ is a generalized vector of external forces and moments.

3. Experiment test

In order to verify the accuracy of the simulation results, the common tanker on the market is contrasted with the 6-DOF equation of motion. Real vehicle experiment adopts the single line test method, test pegged diagram as shown in Fig. 5. $L=$ 6 m; $L_1=L_3=$30 m; $L_2=$25 m or $L_2=$40.6 m (the value can be adjusted according to the specific models); Distance: a lane change $D=$15 m. Benchmarking width: $B=$1.2$b+$0.25, $b$ for car width.

Fig. 5Pole layout diagram

The instantaneous maximum moments are contrasted by the real vehicle test and 6-DOF equation of motion (Fig. 6(a)). Two curves are the change tendency of the instantaneous maximum moments when the charging ratio from 40 % to 90 %. Both have the same change trend, which can verify the feasibility and applicability of the 6-DOF equation of motion.

Fig. 6Contrast test

a) Instantaneous maximum moment of the tanker and the model

b) Lateral force of Fluent and the model

In the field of computational liquid dynamics, the determination of solving domain and the control equation of discretization is the key point. According to the difference principle and way of discrete, the common discretization method can be divided into: finite difference method, finite element method (FEM) and finite volume method (FVM). The FVM is the most widely used and FEM is relatively mature. We used the Fluent software to verify the 6-DOF equation of motion, to simulate on the lateral force of tanks. It can be seen under the lateral acceleration of the time-varying, the use of Fluent software and dynamic model for the liquid sloshing values coincide well (Fig. 6(b)).

According to the different filled level tank liquid by simulating, the coordinates of the mass centre to change with the force and moment about $X$, $Y$, $Z$-axis. The related parameters in the simulation are listed as below, the friction coefficient is set as 0.5; the turning radius is 25 m; the truck mass is 5000 kg; the liquid mass is 4000 kg; suspension stiffness is 1000 kN/m; suspension damping is 15 kN/m; half track width is 2 m; half spring spacing is 3 m; the major axis and minor axis of ellipse are equal to 2.4 m and 1.8 m, respectively.

Fig. 7The moment under the filled level 20 %-50 %

a) The values of $Y/t$ (The force about $Y$-axis)

b) The values of $M/t$ (The moment about $Y$-axis)

c) The values of $N/t$ (The moment about $Z$-axis)

1) $h\in \left(0,b_t\right]$, the filled level is 20 to 50 %, and $f=h/2b_t$.

The numerical simulations are used for 20% to 50% filled level of the tank. Excitation parameters are that the left front and left rear wheel are the sine functions, the right front and right rear wheel are the cosine function. Fig. 7 shows the force about $Y$-axis. Fig. 7(b), (c) are shown the moments about $Y$-axis and $Z$-axis.

2) $h\in \left(b_t,2b_t\right]$, the filled level is 50 % to 100 %.

Due to the liquid mass center equation about liquid filled level 0-50 % and the equation of 50-100 % are different, that draw the following series of graphics (Fig. 8(a)-(c)). Fig. 9(a) shows the force about the $Y$-axis. Fig. 8(b), (c) is shown the moments about $Y$-axis and $Z$-axis.

4. Discussion

To solve the partially-filled tank truck liquid sloshing dynamics, this study established tank of liquid C.G. motion equation. Moreover, six degrees of freedom nonlinear liquid sloshing model were established by applying the Euler transform ideas. Finally, two sets of equations are coupled by using MATLAB ODE to simulate.

1) $h\in \left(0,b_t\right]$, the filled level is 0 to 50 %:

At the beginning, the lateral force of four different situations is basically the same they were slowly rising trend (Fig. 7(a)). After 2.5 s, the lateral force of $f=$ 0.4 (the filled level 40 %) is higher than other three situations. Then, $f=$ 0.3 for the liquid level is the smallest in the lateral force. At the same time it has been at a lower level, which accords with the actual situation.

The moment of different filled level related to the $Y$-axis is increased at the simulation phase (right-hand rule) (Fig. 7(b)). Four kinds of situations are presented slowly rising trend. But $f=$0.5 corresponding moment value is slightly larger than others. Four kinds of filled level of the yawing moment are equal that is the first time (at 1.8 s). In the process of the whole simulation, four equal yawing moments appear many times. In the process of the whole simulation, the corresponding moment value about $f=$ 0.5 is faster than others growth or decline, took two of local minimum or maximum on the four levels.

Fig. 8The moment under the filled level 60 %-90 %

a) The values of $Y/t$ (The force about $Y$-axis)

b) The values of $M/t$ (The moment about $Y$-axis)

c) The values of $N/t$ (The moment about $Z$-axis)

The four different liquid level corresponding moment’s absolute values were in a state of rise about the rotational moment of the $Z$ axis (Fig. 7(c)). After the start of about 1 s, four in moments of liquid level are equal to that of the latter continue to rise. Then, the $f=$0.5 corresponding moment value rises faster, always maintain a high state.

2) $h\in \left(b_t,2b_t\right]$, the filled level is 50 % to 100 %:

The tank of liquid sloshing numerical simulation results for the filled level 50 % to 90 % as shown Fig. 8(a)-(c). The shape of the graphic below is some similar, but there are differences in values. In Fig. 8(a), when simulation to 2.5 seconds, four different lateral force all showed a trend of rise, but the corresponding force about $f=$0.9 rose slowly. It demonstrates that, the mass centre of high liquid filled level relatively small variations, the lateral force is small in the process of liquid sloshing. As is Fig. 8(b), $f=$ 0.6 corresponding moment decreased rapidly and reach the maximum points. In addition, Roll moment curve about $f=$ 0.9 (the filled level 90 %) is more gentle, has good ability to resist roll. In Fig. 8(c), about the rotational moment of the $Z$ axis in the filled level from 50 % to 90 %, the four different liquid level corresponding moment’s absolute values were in a state of rise. After the start of about 1 s, moments of several liquid filled levels are equal to that of the latter continue to rise. In the meantime, the $f=$0.9 corresponding moment value rises faster and always maintain a high state.

5. Conclusion

The quality of the tank liquid will increase with the increase of liquid ratio. Its mass center height increases, the external incentives for the same, the tank lateral sloshing larger than the small ones. The more liquid tank is mass center height increase at the same time.

In the case of liquid ratio is less than 50 %, the roll, yawing and pitching moment smaller than others. But in the case of liquid ratio is more than 50 %, three axis force and moment obvious different. Because the free surface is small, high liquid filled levels in the center of mass of liquid tank is relatively small variations, received the lateral force is small, but quite a few requirements of speed. Also, $f=$ 0.9 (90 %) for the liquid level of the roll torque curve is more gentle, has good ability to resist lateral tilt.

This work in view of the friction coefficient is 0.5, the turning speed of is 80 km/h, turning radius of is 25 m, elliptical tank truck coupled simulation analysis, the model better reflects the external excitation caused by the sloshing of the tank of liquid sloshing in the tank. It provides a theoretical basis for further study on the different surface friction coefficient, turning speed and turning radius of the various cross section shape and size of the liquid tank design.