The field of fluid dynamics has been extensively studied for decades using, for the most part, Eulerian approaches. Another area of application that has recently seen significant advances is vehicle dynamics, which is often examined using MBS algorithms based on a total Lagrangian approach. Nonetheless, fluid-vehicle interaction impacts many areas of science and technology including rail, highway, aerospace, and marine transportation. Although materials, including crude oil and other hazardous materials (HAZMAT), are transported using a variety of methods, including rail, shipping vessels, and pipelines, transportation by highway vehicle dominates the industry, generating more revenue and creating more jobs than the other modes of transportation combined, as shown by the data presented in Table 1. Due to the extent of public roads in the US and the sheer volume of freight vehicles, the tonnage of materials transported using highway vehicles far outweighs all other methods. This is true for both non-hazardous and hazardous materials, as shown in Table 2 [1].
Rollovers are more common in tanker trucks than passenger vehicles because trucks have a higher center of gravity. Rollovers can occur due to a variety of reasons, including vehicle and road conditions, load size, and the most common, driver error, which accounts for up to 78% of tanker truck rollovers [2]. Hazardous materials are regularly transported by tanker trucks, and accidents in which the tank is compromised and the contents are released can lead to damage to the environment and the surrounding infrastructure, fires and explosions, and civilian injuries and casualties [3,4]. In the last decade alone, highway transportation accidents comprised the majority of all HAZMAT incidents, with 144,296 out of a total of 166,494 incidents; other incidents include air, railway, and water transportation accidents. Highway accidents have also proven to be the most deadly and costly, accounting for 100 out of 105 documented fatalities and 1,520 out of 2,129 injuries, at a cost of $6.1 billion out of $8.2 billion in damages [5]. Therefore, thorough testing and virtual prototyping are necessary to ensure better vehicle design and stability. However, because physical prototyping is expensive, inefficient, and time-consuming, it is necessary to develop accurate predictive models to investigate the effect of liquid sloshing on the dynamics of highway vehicles subject to different loading conditions and motion scenarios.
Although recent advances allow for modeling more accurate fluid behavior, most commonly used models are insufficient in adequately capturing the dynamics of the fluid in complex motion scenarios, particularly in the cases of three-dimensional finite rigid body rotations. Early sloshing models represented the fluid as a series of planar pendulums or massspring systems [6 -9]; spherical and compound pendulums were later used to capture nonlinearities in the motion and damping was added to include the effect of energy dissipation [10]. Discrete inertia models have been used extensively in studying sloshing dynamics in the aerospace industry since 1960s [8,11,12]. However, while these discrete inertia models have been improved over time, such models cannot be used to accurately capture the change in inertia due to a change in fluid shape and the complex dynamics that results from the vehicle motion [13]. Furthermore, the discrete rigid body models do not allow for modeling the continuous free surface of the fluid, and it has been found that such models significantly under-predict the maximum amplitude of oscillation and the sloshing frequency [14]. Nonetheless, discrete pendulum models are still being used to study the effect of fluid sloshing on vehicle dynamics and stability because of the difficulties encountered in integrating computational fluid dynamics (CFD) and vehicle modeling algorithms [15]. This is attributed to the fact that CFD investigations are mostly focused on the use of the Eulerian approach and do not consider the effect of three-dimensional rotations, which cannot be accurately captured using existing incremental FE formulations. Consequently, most FE fluid models are simplistic due to the inherent difficulties in accurately capturing the fluid behavior, particularly in applications related to vehicle dynamics. For example, a sloshing fluid will experience large deformation and finite rotation, and many commercial FE formulations are not capable of handling such behaviors accurately or efficiently. Furthermore, because of the interaction between fluid and the tank walls, liquid sloshing results in variety of dynamic behaviors, including symmetric and asymmetric motion, planar and non-planar motion, and rotational and irrotational motion [16 -18]. Furthermore, in flexible body dynamics, the centrifugal forces which result from curve negotiation are not simply measured by the rigid body dynamics equation 2 s mV r , where m is the mass of the body, s V is the forward velocity, and r is the radius of curvature of the curve [19]. In flexible body dynamics the centrifugal forces take a more complex form that depends on the body deformation. Additional difficulties arise when more accurate or complex fluid models are integrated into full vehicle models. Many FE and CFD formulations, for example, do not easily lend themselves to integration with MBS algorithms as discussed in the literature [18].
In order to address the challenges encountered in the analysis of the important economic, environmental, and safety problem of liquid sloshing, this investigation makes the following specific contributions:
1. A formulation that correctly captures the geometry of the fluid and tank is used in order to accurately represent the distributed inertia and elasticity of the fluid. In order to develop these new and unique sloshing geometry models, ANCF elements that produce accurate geometry are used, eliminating the need for using B-spline and NURBS (Non-Uniform Rational B-Spline) representations for developing the complex fluid geometry.
The effect of the initially curved fluid geometry, which cannot be captured accurately using existing FE formulations, is properly accounted for in this investigation, leading to a systematic integration of the geometry and analysis by adopting one fluid mesh from the outset. Such an important goal cannot be achieved using other MBS formulations that employ modal representation for the fluid displacements, as in the case of the floating frame of reference (FFR) formulation [18].
2. The ANCF geometry/analysis mesh developed is used to formulate the inertia forces using a non-modal continuum-based approach. Proper definition of the inertia forces is necessary in order to be able to predict the effect of the sloshing on the vehicle dynamics and stability. In particular, a continuum-based and general definition of the centrifugal forces in terms of the fluid displacement is developed and used to shed light on the approximation made using the simple rigid body dynamics formula 3. An ANCF fluid/tank car walls penalty contact formulation is developed and used to determine the generalized contact forces associated with the ANCF nodal coordinates which include absolute position and gradient vectors. The penalty contact formulation developed in this investigation takes into account the fluid large displacement and complex geometry that result from the sloshing effect.
4. It is shown in this investigation how general constitutive fluid models can be developed and integrated with ANCF complex fluid geometry models, thereby opening the door for future investigations that focus on adopting new and highly nonlinear constitutive fluid models as well as experimenting with different tank designs that have different, complex, and unconventional geometries. In so doing, the field of liquid sloshing can be significantly advanced to a new level.
5. The analysis presented in the paper demonstrates for the first time how an ANCF liquid sloshing model can be integrated with an MBS system computational algorithm that ensures that the kinematic algebraic constraint equations are satisfied at the position, velocity, and acceleration levels. Such new ANCF fluid/MBS algorithms will allow for investigating a large class of liquid sloshing problems that cannot be solved using existing approaches. The purpose of this analysis is to create a high fidelity model which is capable of capturing more details than can be described by existing modeling methods.
It is important to note that simple models can still be valuable if real-time simulations are required. In these cases, both simple vehicle and fluid models can be used to significantly reduce the computer simulation time. High fidelity continuum-based models, on the other hand, are necessary in order to account for the distributed inertia and viscoelasticity of the fluid. The interaction between the rigid tank walls and the ANCF fluid is formulated using the penalty approach. The MBS model includes a suspension system and Pacejka's brush tire model is introduced to represent the ground-tire contact [20]. Specified motion trajectories are used to examine three different working conditions -deceleration under straight-line motion, rapid lane changing, and negotiating a curve. Reduced integration is used to increase computational efficiency when the fluid viscosity forces are calculated.
7. The results of this investigation show that sloshing has the effect of increasing contact forces on some wheels and decreasing contact forces on other wheels. Severe sloshing behavior can cause vehicle instability; in extreme cases, wheel lift and vehicle rollover may occur.
In liquid sloshing problems, the main focus is on the nominal motion of the fluid that affects the vehicle dynamics. Therefore, turbulence is not an issue and is not normally considered in liquid sloshing investigations since it has a negligible effect on the vehicle dynamics and stability. For
In this section, a simplified planar vehicle model subjected to discrete forces is analyzed in order to have an understanding of how the contact forces on the tires change as a tanker truck enters a curve. A force diagram for this model during straight-line motion is presented in Fig. 1, where , w t F is the tank gravity force at the tank center of mass located at a vertical distance t z from the ground, , F and the lateral displacement of the fluid f x will be equal to 0.
Using this simple analysis, one can examine how the contact forces on the tires change when a vehicle enters a curve. In Fig. 3, the steady-state normal force equations are used to calculate the contact forces for the first 0.7s, then the constant-radius curve contact force equations are used for the following 9.3s. This represents a vehicle driving in a straight line initially before entering a constant-radius curve at 0.7s, where it remains for the rest of the simulation. While the results of this figure, obtained using the simple analysis and the simple force equations previously presented in this section, do not capture the oscillations of the fluid because the lateral shift of the fluid f y is assumed to remain constant for simplicity, it is evident from these results that the contact force on the outer tire increases and that on the inner tire decreases. This change is due to the lateral shift of the center of mass of the fluid, which is a result of the outward inertia force acting on the fluid. The lateral shift of the fluid and thus the outward inertia forces act to increase the roll moment and thus increase the contact force on the outer tire.
These simplified results can be made more realistic by using simulation results in the equations instead of constant theoretical values. By replacing the position of the center of mass of the fluid f y and f z and the centrifugal force on the fluid , C f F with the simulation results that will be presented in detail in Section 9.3, the resulting contact forces calculated by the previously derived equations will capture the sloshing behavior. This effect is evident from the results presented in Fig. 4, where the contact forces oscillate with time due to the oscillatory motion of the liquid. The discontinuity in the plot is due to the fact that the theoretical calculations assume a sudden change from straight-line to constant-radius curve trajectories. In more realistic scenarios, a spiral segment is used to connect the straight and curved sections in order to ensure a smooth transition.
Inertial forces play an important role in the dynamics and stability of a vehicle negotiating a curve. The centrifugal force is exerted on the vehicle in the outward normal direction of the curve. If the bank angle of the curve is zero, the only opposing force is the inward lateral friction force due to the contact between the tires and the ground. When the bank angle is different from zero, the centrifugal force of a vehicle with mass m is also opposed by the component of the gravity force along the normal to the curve. If the rigid body assumptions are used and additionally the vehicle forward velocity s V along the tangent to a curve of radius of curvature r is assumed constant, one must have an upper limit on the velocity s V , called the balance speed, such that When ANCF finite elements are used, the expression of the outward inertia force differs significantly from the expression used in rigid body dynamics. For ANCF finite elements, the vector of nodal coordinates can be written as the sum of two vectors as
e e e , where o e is the vector of nodal coordinates before displacement and d e is the vector of displacements that include large liquid reference displacements including finite rotations as well as the liquid deformations. Therefore, the outward inertia force, as will be demonstrated in this section, becomes function of the liquid motion and the simple expression 2 s mV r is no longer applicable for the calculation of the balance speed or for accurate force analysis during curve negotiations.
Furthermore, the vector o e can be used to systematically account for the initial curved geometry of the liquid. As described in the literature, this can be accomplished by using the matrix of position vector gradients o J , where In order for the vehicle to safely remain on the road, the outward inertia force must not exceed the sum of the inward friction and gravity forces. Although the centrifugal force on a rigid body negotiating a curve takes a simple form, as previously mentioned, the same expression does not apply to curve negotiation of a flexible body, because such a force expression is function of the deformation [17]. In general, the outward inertia force in F of a flexible body or an ANCF finite element negotiating a curve is defined as d
that the component of the acceleration along the tangent to the curve will not contribute to the outward force vector. When ANCF finite elements are used, the absolute acceleration vector of an arbitrary point can be written as r Se . If the flexible body is discretized using e n ANCF elements, the outward inertia force vector that must be used to define the vehicle balance speed can be written as
where superscript j refers to the element number. One can also write integration can be systematically used to evaluate the outward inertia force
analytical integration of the element shape functions is to be avoided. In this case, one can create a mesh of p n points on the flexible body and if an assumption is made that the mesh consists of only one type of ANCF elements, then approximation of in F can be written as
where k m is the lumped mass associated with the mesh point k , k S is the assembled matrix of the element k j j k j S S x B matrices, and jk
x is the vector of the element spatial coordinates
evaluated at the mesh point k that corresponds to element j .
Alternatively, one can use the moment of mass to write d
, where m is the total mass of the liquid, C r is the global position vector of the liquid center of mass, and r Se
when ANCF finite elements are used. It follows that
which can be simply written as
, and
. Therefore, the outward inertia force vector can be written in an alternate form as
Because of the liquid oscillations, C r will not remain constant relative to the curve, and as a consequence, the outward inertia force is not in general constant as in the case of a rigid body negotiating a curve.
In this section, the development of the initially curved ANCF geometry of the fluid that assumes the shape of a rigid cylindrical tank is discussed. The tank used in this investigation has a cylindrical geometry as shown in Fig. 5, and therefore, it is required for the ANCF fluid mesh to have the same the shape of the container it fills and at the same time represent different levels of the free surface. The use of the ANCF absolute positions and gradients as nodal coordinates allows for efficient shape manipulation and for obtaining the accurate complex geometry without the need for using the CAD B-spline and NURBS representations that have rigid recurrence structure [21,22]. As previously mentioned, in the ANCF description, the assumed displacement field can be written as ( , ) ( ) ( )
S x e , where r is the global position vector, [ ] T x y z x is the vector of the element spatial coordinates, t is time, S is the time-independent element shape function matrix, and e is the vector of the element nodal coordinates that include absolute position and gradient coordinates [23]. The superscript j that refers to the element number is omitted here for notational simplicity. The vector of element nodal coordinates e can be written as
e e e , where o e is the vector of nodal coordinates in the reference configuration and d e is the vector of nodal displacements. The assumed displacement field can then be written as
S x e e . Using the general continuum mechanics description By choosing the elements in the vector o e appropriately, initially curved structures can be defined in a straightforward manner using ANCF elements [21].
The fully-parameterized ANCF solid element [24,25], based on an incomplete polynomial representation, is used in this investigation to represent the fluid by applying the proper fluid constitutive model which will be discussed in Section 5. In this case, the arbitrary fluid material point on element j can be written as
where I is the 3 3 identity matrix; detailed shape function and nodal coordinate expressions can be found in Appendix A of this paper. For example, consider the element j which has the initially curved structure shown in Fig. 6. The matrix of position vector gradients at node k can be written as
. For the specific element geometry shown in Fig. 6, by adjusting the magnitude of the gradient vector
without changing the gradient vector orientation, the position vector gradients at node 5 will be
where is the stretch factor used to represent the stretch of the edge; the value of can be obtained by taking the ratio between the arc lengths of curves 5-8 and 1-4. Following this procedure, the complex geometry of the fluid structure can be created, as shown in Fig. 7. The mesh used in this investigation consists of 48 ANCF solid elements and the mesh has a total number of degrees of freedom of 1260.
A general ANCF fluid constitutive model that can account for the initially curved configuration is developed in this section. The proposed fluid model ensures the continuity of the displacement gradients at the nodal points and allows for imposing a higher degree of continuity across the element interface by applying algebraic constraint equations that can be used to eliminate dependent variables and reduce the model dimensionality at the prepossessing stage [25]. In order to describe the fluid-structure interaction, the penalty approach, described in Section 6, is used to evaluate the contact and friction forces between the fluid and the rigid tank. By using the non-modal ANCF approach, the fluid elastic forces can be formulated without imposing restrictions on the amount of deformation and rotation within the elements. Figure 8 shows the three configurations of the fluid; the straight, curved reference, and current configurations. As previously mentioned, the volume of the fluid in the curved reference configuration o V is related to the volume in the straight configuration V using the relationship d d
is the determinant of the matrix of position vector gradients
Therefore, integration with respect to the reference domain can be converted to integration with respect to the straight element domain. This allows for using the original element dimensions to carry out the integrations associated with the initially curved configuration. The matrix
x is constant and can be evaluated at the integration points using the ANCF element shape function and the vector of nodal coordinates in the reference configuration [21]. The matrix of position
, where
e J r x Se x . The relationship between the volume defined in the current configuration v and the volume in the curved reference configuration o V can be written as
The linear fluid constitutive equations can be defined using the Cauchy stress tensor and can be assumed as
where the temperature effect is neglected and the fluid is assumed to be incompressible, σ is the symmetric Cauchy stress tensor, p is the hydrostatic pressure, and are Lame's material constants, I is a 3 3 identity matrix, tr refers to the trace of a matrix, and D is the rate of deformation tensor [23,26]. If the incompressibility condition is imposed using a penalty method, the first two terms will vanish and the constitutive model reduces to 2 σ D . It is convenient to use the second Piola-Kirchoff stress tensor since it is associated with the Green-Lagrangian strain tensor defined in the reference configuration. One has
. For an arbitrary element j in the fluid body, the virtual work of the fluid stress forces can be written as
where
The virtual work of the fluid viscous forces can then be written as
where T j j j r C J J is the right Cauchy-Green deformation tensor, and upon using the identity 1 1
, the vector of generalized viscosity forces j v Q associated with the ANCF nodal coordinates can be written as
In this case, the integration over the current configuration domain is converted to integration over the straight configuration domain.
The incompressibility condition is imposed using the penalty method. Figure 8 shows that the volume relation between the reference and current configuration is d d
therefore, 1 with the ANCF nodal coordinates that result from imposing the two penalty conditions can be defined as
where
. Knowing that
r r r r r r , j j J e can be written more explicitly, by differentiating any of the three expressions for j J with respect to j e , as
By defining the generalized forces associated with the fluid element coordinates j e , the generalized forces associated with the fluid body coordinates e can be obtained using a standard FE assembly procedure.
The fluid should remain inside the tank regardless of the severity of the sloshing and these boundary conditions of the mesh can be defined in multiple ways. One can impose constraints on the boundary nodes, using either Lagrange multipliers or elimination of dependent variables.
Because this method is often more computationally expensive, in this investigation, the penalty method is used to formulate the interaction between the fluid body and the rigid tank walls. The tank deformation is not considered in this analysis because the main focus of this investigation is on studying the sloshing. Figure 9 shows the contact geometry in the radial direction; the radius of the tank is t r , superscript t refers to the tank, f r is an arbitrary point on the fluid body, and t R is the position vector of the tank reference point defined in the global coordinate system. The position of an arbitrary point of the fluid defined in the tank local coordinate system can be written as
, where t A is the 3 3 transformation matrix which defines the tank orientation, and the bar notation means the vector or matrix is defined in the body local coordinate system. The inequality
implies that penalty forces must be applied in order to prevent the fluid from penetrating the tank walls. In this case, the penetration can be evaluated as
The unit normal n at the fluid/tank contact point can be defined as
. The location of the contact point on the tank wall with respect to the tank coordinate system can be defined as 1 0 0
. This equation can be used to define the global position vector of the contact point on the tank as
r R u . The relative velocity vector can be defined as ft Knowing the ANCF element j on which the fluid/tank contact point lies, one can develop an expression for the generalized penalty contact forces associated with the fluid element nodal coordinates. The virtual work of the penalty force acting on the fluid and tank can be written as
, which can be written as
, where j S is the element shape function matrix evaluated at the contact point, t u is the skew symmetric matrix associated with the vector t u , and t G is the matrix that relates the absolute angular velocity vector t ω of the tank to the time derivatives of the tank orientation parameters t θ , that is, t t t ω G θ . It follows that the generalized reaction forces exerted on the element j of the ANCF fluid body can be written as
, while the generalized penalty forces exerted on the tank and associated with the tank reference coordinates t R and t θ are given, respectively, as
A similar procedure can be used to evaluate the interaction forces between the fluid and the rigid tank in the longitudinal direction. This contact geometry is shown in Fig. 10.
The MBS model used in this numerical investigation consists of 21 bodies which have 147 absolute coordinates because Euler parameters are used to describe the body orientations. These bodies are subjected to 115 constraint equations, leading to a model with 32 degrees of freedom.
The 10 tires are modeled using Pacejka's brush tire model, which is discussed in Section 7.1 [20].
The four-bar Ackermann steer axle, which allows the truck to turn, is described in Section 7.2.
Other bodies include two rear axles, the cab, the tank, the frame rails, and a ground body. Inertial properties for the individual bodies are listed in Table 3; the products of inertia are assumed to be zero. The wheelbase of this vehicle is 4.064 m and the track width is 1.939 m. The front wheels are attached to the steer axle using the steering knuckles with revolute joints; this allows the wheels to rotate about the lateral axis to produce the desired forward motion, and about the vertical axis to allow the vehicle to turn. For simplicity, the steering knuckles and tie rod arms are modeled as a single body. The rear wheels are connected to the drive axles with revolute joints allowing rotation about the lateral axis. The cab and tank are rigidly attached to the frame rails and these three bodies are assumed in this model to represent the chassis (sprung mass). The capacity of the tank is roughly 3,000 gallons, which is typical of a medium-duty commercial vehicle that services residential areas. In order to induce the most extreme sloshing scenarios, the tank is assumed to be half-filled with water (with viscosity of 0.001 kg/m.s [28]). The penalty coefficients used in this investigation to enforce the incompressibility conditions 1 4.
In this investigation, the forces exerted on the tires by the ground are calculated using the brush tire model [20]. Four coordinate systems which define the orientations of the tire body i and the ground body j are introduced. The ground coordinate system which describes the orientation of the ground body is given by the matrix j j j j A i j k , where the columns j i , j j , and j k are unit vectors along the ground coordinate axes j j j X Y Z . In this investigation, j A is assumed to be the identity matrix except in the case of uneven terrains, such as a hill, bumpy road, or inclined ramp. The three other coordinate systems are used in the tire formulation and are depicted in Fig. 11. The tire coordinate system defined by the matrix
attached to the center of the tire and rotates with the tire. The matrix
describes the intermediate tire coordinate system (ITCS) which is also rigidly attached to the center of the tire, but does not share the pitch rotation with the tire. The axes of i o A are defined as
Finally, the contact point coordinate system located at the contact point between the ground and the tire is defined by the matrix
A , the contact point between the tire and the ground can be defined as
, where i R is the global position of the ITCS origin and i oc u is the position of the contact point with respect to the ITCS origin. The velocity vector of the contact point can then be obtained by differentiating the position vector with respect to time, and is defined as
, where i R is the velocity of the ITCS origin, i ω is the absolute angular velocity vector of the tire defined in the global coordinate system, and
The brush tire model accounts for the normal, lateral, longitudinal, and rotational friction forces, as well as an aligning torque, also referred to as the spin moment in rail vehicle dynamics.
In this investigation, it is assumed that the material properties are the same in the lateral and longitudinal directions, which is a simplifying assumption often made in the literature [20,22,29,30,31]. The normal force in the contact point coordinate system is calculated simply as
where K is the radial stiffness coefficient of the tire, z is the vertical penetration of the tire with the ground, C is the radial damping coefficient, and z is the rate of change of the penetration. In order to determine the lateral and longitudinal friction forces, the slip angle and slip ratio vector ξ are needed. First, the slip velocity s v is defined as
, where rc v is the velocity of the tire with respect to the ground at the contact point. The slip ratio vector ξ is 0 0
, where r v is the tire forward velocity. The slip angle is determined as
. The model parameter is defined as The lateral and longitudinal friction forces on tire body i in the contact point coordinate system are determined by
is a simplifying model parameter. The aligning torque i z M is also calculated in the contact point coordinate system depending on the slip angle , as
in the case of pure sliding.
To determine the moment due to the rotational friction force The forces and moments obtained in this section for the tire can be defined in the appropriate coordinate system for the inclusion in the Newton-Euler equations that govern the motion of the tire which is treated in the brush model as a rigid body.
In order for a vehicle to be able to negotiate a turn with minimal tire scrub, the Ackerman steering mechanism is often used. The Ackermann steering system is a four-bar mechanism which allows the wheels to be oriented at different angles with respect to the forward direction.
The four-bar mechanism consists of the front axle as the ground link, two tie rod arms, and a tie rod, as shown in Fig. 12. The tie rod is connected to the tie rod arms with spherical joints to avoid over-constraining the mechanism. By specifying the geometry of the Ackerman mechanism, each of the front wheel forward velocity vectors remains tangent to a circular arc whose origin is located at the instantaneous center of rotation of a line element connecting the centers of the two wheels, thus reducing tire scrub. The linkage geometry can be defined by two equations, cos h r and 2 sin
, where h is the distance between the axle and the tie rod, r is the length of the tie rod arm, is the angle between the tie rod arm and the normal to the axle, s is the length of the tie rod, and t is the length of the axle, as shown in Fig. 13.
Because the length of the axle and the wheelbase are known for the vehicle model used in this investigation to be 1.939 m and 4.064 m, respectively, can be calculated as 13.42 o . In the optimization study by De-Juan et al. [32], it was found that for an axle length of 1.5 m, the optimum tie rod length and tie rod arm length are, respectively, 1.27 m and 0.3 m. To determine the remaining parameters of the steering mechanism, it was assumed that the steering mechanism dimensions are proportional to the dimensions of the optimal mechanism geometry obtained in the study by De-Juan et al. [32]. Using this assumption, the geometry parameters obtained were r 0.3878 m, s 1.7590 m, and h 0.3772 m.
The suspension system of the truck is modeled using linear spring-damper elements. The springdampers are oriented to provide restoring forces in the longitudinal, lateral, and vertical directions -the vertical springs are used to support the weight of the chassis, cab, tank, and fluid, while the longitudinal and lateral springs are used to prevent relative motion in the longitudinal and lateral directions, respectively. The spring-dampers are located at each end of the three axles, resulting in 18 elements total. The vertical springs on the rear axles are initially compressed to aid in supporting the weight of the tank and fluid and to minimize oscillations at the beginning of the simulation when the vehicle reaches equilibrium.
In numerical simulations, two methods can be used to produce the MBS motion; the first is to apply forces on the system components, while the second is to specify motion trajectories using algebraic constraint equations. The latter approach of using constraint equations is more appropriate when it is required to precisely follow certain trajectories and correctly capture their geometry, as is the case in this investigation. In order for the tanker truck model to follow different specified paths necessary to create the motion scenarios to be investigated in this study, trajectory coordinate constraints must be imposed. Three coordinate systems are used to define the type of the trajectory coordinate constraint used in this investigation: the global coordinate system XYZ , the trajectory coordinate system ti ti ti X Y Z , and the body coordinate system ir ir ir X Y Z , as shown in Fig. 14 [33,34]. Six trajectory coordinates i p are used to specify the motion of a body, where
, i s is the arc length along the user-specified trajectory, ir y and ir z are the lateral and vertical displacements of the body with respect to the specified trajectory, and ir , ir , and ir are the three Euler angles describing the relative rotations of the body coordinate system with respect to the trajectory coordinate system.
Because a curve can be completely defined using one parameter, Frenet frame geometry is employed to write the matrix ti A that defines the orientation of the trajectory coordinate system in terms of three Euler angles ti , ti , and ti which can be expressed in terms of the arc length parameter i s as
, where k C is the constraint function corresponding to the trajectory coordinate R,
, and ( ) f t is the time-varying function defining the values of the trajectory coordinate , i R k p . For example, to constrain the vehicle to move along a specified path with a constant forward velocity s V , the constraint applied to the body coordinate system of the front axle can be written as 0
, where ir R is the global position of the body coordinate system with respect to the trajectory coordinate system and ir u is the position vector of the center of mass of the body in the trajectory coordinate system. For a given set of absolute Cartesian coordinates, this set of nonlinear equations can be solved iteratively to determine the arc length parameter i s as well as the coordinates ir y and ir z . The constraint Jacobian matrix associated with the absolute Cartesian coordinates can also be systematically developed and used to enforce the constraints at the position, velocity, and acceleration levels. The driving constraint forces that produce the desired motion can be determined using Lagrange multipliers and the trajectory coordinate constraint Jacobian matrix [33].
The vectors and matrices defined in the previous sections enter into the general formulation of the equations of motion for the MBS vehicle model, which may include rigid and flexible bodies.
The equations of motion can be written as [33]:
where r M and e M are the mass matrices associated with the rigid and flexible ANCF body coordinates, respectively, r q C , e C , and s C are the constraint Jacobian matrices associated with the rigid, elastic, and non-generalized trajectory coordinates, respectively, r q , e , and s are the accelerations of the rigid reference, elastic, and non-generalized trajectory coordinates, respectively, λ is the vector of Lagrange multipliers associated with the constraints, r Q and e Q are the vectors of generalized forces associated with the rigid and elastic coordinates, respectively, and d Q is a vector resulting from the second time derivative of the vector of constraint equations. The equations of motion are solved numerically using the Adams-Bashforth integration technique, and a solution algorithm that ensures that the constraint equations are satisfied at the position, velocity, and acceleration levels. Because the equations of motion are second order differential equations, two sets of initial conditions (coordinates and velocities) are required to obtain a unique solution. The initial coordinates and velocities of the rigid bodies, 0 r q and 0 r q , the initial coordinates and velocities of the nodes of the flexible ANCF bodies, 0 e and 0 e , and the initial non-generalized coordinates and velocities 0 s and 0 s are user-specified and known at the beginning of the simulation. The initial velocities of all bodies and nodes vary depending on the maneuver being considered in the numerical simulation. These maneuver scenarios are presented in Table 5.
The dynamic behavior of the tanker truck model, created in the MBS dynamics software SIGMA/SAMS and shown in Fig. 15, was examined using three different motion scenarios, each of which produces different fluid motion. In the first scenario, the truck decelerates under straight-line motion, such that the fluid primarily exerts longitudinal forces on the tank. In the second scenario, the truck performs a lane change, which causes the fluid to exert alternating lateral forces on either side of the tank. In contrast, in the third scenario the truck is assumed to negotiate a wide curve, such that the lateral motion of the fluid is continuous and exerted on one side of the tank only. Steady state is achieved before the simulation results are reported in order to eliminate the transient effects. Another model in which the fluid is represented by a rigid body with equivalent inertial properties was also created; the rigid fluid body is rigidly attached to the tank so that the sloshing motion is prevented but the fluid inertia is correctly accounted for. This model was also examined using the same three scenarios so that the effect of the fluid sloshing on vehicle dynamics can be isolated and evaluated.
For this scenario, the truck begins at a highway speed of 55 mph and brakes to 20 mph as seen in the velocity -position plot in Fig. 16. This scenario could occur if braking is suddenly applied in an attempt to avoid a rear-end collision. As seen in Fig. 17, the sloshing phenomenon is clearly evident as the fluid moves in the longitudinal direction towards the front end of the tank as a result of the sudden braking. It is important to note that the section on the top of the tank is an exterior section that the fluid cannot enter (refer to Fig. 5); therefore, the fluid fills the front portion of the tank during braking.
The normal forces increase on the front tires and decrease on the rear tires, as evident in Fig. 18, because the center of mass of the fluid moves towards the front of the tank, as seen in the plot of the fluid center of mass presented in Fig. 19. A second peak occurs in the contact forces on the front tires near 7s due to the rebounding motion of the fluid after it impacts the rear of the tank and again sloshes longitudinally towards the front of the tank. Although the normal forces also increase on the front tires and decrease on the rear tires in the equivalent rigid body model due to the shift in inertia of the chassis on the suspension system, also seen in Fig. 18, the magnitudes are much less significant because the center of mass of the rigid fluid body remains constant with respect to the tank and the relative displacement of the chassis on the suspension system is small compared to the longitudinal displacement of the flexible fluid. The fluid free surface of the ANCF model returns to a flat shape once the truck reaches the lower speed, as seen in Fig. 20, and the normal forces return to approximately equilibrium in both models, as evident by the results presented in Fig. 18.
In the second motion scenario, the truck performs a lane change over a standard-width US highway lane of 3.7m, as seen in the plot of lateral -longitudinal position in Fig. 21. The lane change is completed in a relatively short time of 4s so that the fluid sloshing readily occurs, as seen in the series of images in Fig. 22 depicting the change in the free surface of the fluid mesh.
The shift in the center of mass of the flexible fluid mesh as the truck negotiates the turns of the lane change causes the normal forces exerted on the tires on the outer edge of the curve to be greater than those exerted on the inner tires, as seen in Fig. 23. For the first half of the lane change, the left tire (Fig. 23a) is the inner tire, and for the second half of the lane change, the right tire (Fig 23b ) is the inner tire. This effect is also evident in the rigid body model due to the outward inertia forces, however the peak forces exerted are greater for the ANCF model than for the rigid body model due to the sloshing behavior, which is evident in the plot of the fluid center of mass presented in Fig. 24. Furthermore, it can be seen that the peaks in the normal forces on respect to the ground is greater for the ANCF model tire, which negates the effect of the decreased normal force. This is clear from the results presented in Fig. 26, where the lateral slip velocity of a left tire is greater for the ANCF model during the lane change.
For the third scenario, the truck negotiates a wide curve as seen in Fig. 27, similar to an onramp or exit ramp of a highway, except that the bank angle is assumed zero for simplicity. While driving along a road of constant curvature and in the case of zero bank angle, the outward centrifugal force on a rigid vehicle due to the curve of the road is counteracted by the lateral friction force exerted on the tires; that is,
where m is the mass of the vehicle, s V is the forward velocity, r is the radius of curvature of the road, is the coefficient of friction between the tires and the road, and g is the gravitational constant. As previously discussed in this paper, the maximum speed at which the vehicle can traverse a curve without sliding can then be calculated as s V gr . For example, for a radius of curvature of 115 ft and a coefficient of friction of 0.7, the maximum calculated speed is 33.9 mph. However, due to the high center of gravity of the truck and increased chance of rollover, in practice vehicles cannot traverse a curve at the theoretical maximum speed. According to the National Highway Traffic Safety Administration (NHTSA), the maximum speed at which an average fully-loaded tractor-trailer can negotiate a curve of 150 ft is 30mph; at greater speeds, the chance of rollover is greatly increased [34]. Therefore, for this analysis, the constant forward speed is chosen to be 25 mph and the radius of curvature of the track is set to 150 ft, as shown in Fig. 27. As previously mentioned, the road is assumed to be flat with no super-elevation; this assumption is consistent with methods used in the literature [15].
The normal and lateral forces exerted on the tires of the rigid body model are larger on the outer tires and smaller on the inner tires, as seen in Figs. 28 and29, respectively. This is due to the roll moment that is exerted from the centrifugal force on the vehicle. Because the radius of curvature is constant and the center of mass of the fluid in the rigid body model is not able to move, the centrifugal force is constant once the truck enters the curve, and thus the normal and lateral forces are constant as well. The contact forces on the tires of the ANCF model, however, overshoot the constant value exerted on the rigid body model tires, and oscillate due to the sloshing motion and the change in the center of mass of the fluid, as seen in Figs. 28 and29.
Furthermore, while the centrifugal force on the rigid fluid is 2 s mV r as discussed in Section 3, it has been demonstrated that this is not the case for flexible bodies [19]. This is evident in Fig. 30,
where the outward inertia force on the flexible fluid mesh oscillates with a maximum amplitude mV r used to calculate the centrifugal force in rigid body dynamics is no longer applicable in the case of liquid sloshing. Figure 32 shows the components of the normalized velocity of the liquid center of mass obtained by dividing by the vehicle forward velocity. Figure 32a shows the dimensionless velocity component tangent to the curve, while Fig. 32b shows the other two components.
The magnitudes of the contact forces are not identical to those predicted by the analytical model in Section 2 because only two wheels were included in that analysis, and the weight of the tank and fluid is actually distributed over 10 wheels. However, the orders of magnitude of the contact forces are the same and the relative changes in the forces were well predicted, and thus the analytical model verifies the simulation results.
A total Lagrangian ANCF fluid formulation that can be systematically integrated with fully nonlinear MBS vehicle algorithms is proposed in this investigation. The new approach can capture the fluid distributed inertia and viscosity, can accurately predict the change in inertia due to the change in shape of the fluid, and can visualize the change in the fluid free surface, unlike other discrete inertia models which do not capture these significant details. The outward forces on the fluid during curve negotiation are derived and it is shown that these forces do not take the same simple form as the case of a rigid body negotiating a curve. As discussed in the paper, accurate modeling of the fluid geometry using ANCF elements can be achieved without the need for using computational geometry methods such as B-spline and NURBS representations which have a rigid recurrence structure unsuitable for MBS analysis. By using the approach proposed in this investigation, one geometry/analysis mesh is used from the outset. The fluid constitutive law and the fluid/tank interaction forces are developed. The penalty method is used to ensure that the fluid remains within the boundaries defined by the tank geometry. Both normal and tangential penalty contact forces are considered in this study. The MBS vehicle model components are described and the dynamics of the vehicle is examined using three contrasting motion scenarios in order to study the effect of sloshing on vehicle dynamics. The braking scenario examines the case of longitudinal sloshing of the fluid, a rapid lane change produces alternating lateral fluid forces on the tank, and curve negotiation sheds light on the case of steady-state outward forces due to the centrifugal effect.
The results presented in this investigation demonstrate that depending on the scenario, the sloshing phenomenon can increase the contact forces on some wheels while decreasing contact forces on other wheels, and this can lead to vehicle instability. In the case of brake applications, the fluid in the partially-filled tank surges forward, causing uneven wheel loading; in cases of severe braking, wheel lift may occur on the rear wheels of the vehicle. This can lead to difficulty controlling the vehicle and increased stopping distances due to lessened road contact, and the decrease in stability may result in jack-knifing for tractor-trailer vehicles. When entering a curve or performing a lane change, the outward centrifugal forces cause lateral displacement of the fluid, which also causes uneven wheel loading. In extreme cases of high speed or small radius of curvature, these changes could be significant enough to induce wheel lift on the inner wheels and increase the possibility of rollover compared to an equivalent truck carrying rigid materials.
Furthermore, in cases where tire friction forces are decreased, such as on wet or icy roads, the possibility for vehicle instability increases even further and drivers must exercise extra caution.
Future analysis can result in defining general rules for speed reduction of a tanker truck entering a curve in order to maximize vehicle stability and driver safety, proposing modifications to the tank geometry to reduce sloshing amplitudes and forces, and studying the effect of viscosity and incompressibility on the fluid dynamic behavior.
The three-dimensional ANCF solid element, with an incomplete polynomial representation, used in this investigation is an 8-node element. The nodal coordinates jk e at the node k of the finite element j can be defined as Where I is the 3 3 identity matrix, j S and j e are, respectively, the element shape function matrix and the vector of nodal coordinates which can be written as