An integrated coupling framework for highly nonlinear fluid-structure problems

Computers & Fluids 60 (2012) 36–48

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

An integrated coupling framework for highly nonlinear fluid-structure problems Qun Zhang a, Baoshan Zhu b,⇑ a b

INTESIM (Dalian) Co., Ltd., Dalian 11602, China Department of Thermal Engineering, State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 9 February 2011 Received in revised form 11 January 2012 Accepted 21 February 2012 Available online 7 March 2012 Keywords: Fluid–solid interaction Morphing Remeshing Strong coupling INTESIM

a b s t r a c t An integrated coupling framework is provided for the transient simulations of large scale highly nonlinear fluid-structure problems with extreme large domain changes. The ALE finite element method is employed for viscous fluid flow, both geometrical and material nonlinearity is considered for dynamic structure analysis. The strong coupling behavior between fluid and solid that causes the high linearity of the coupled system is treated by the strong coupling methods. Advanced morphing and automatic remeshing technology are used for dealing with the extreme fluid domain changes. Time and spatial stabilization, auto-time step and bisection schemes are employed in this research to improve the stability and efficiency of the coupling simulation. Hybrid parallel technologies are adopted for large scale simulation. The fully coupled fluid-structure simulation of a Hydraulic Engine Mounts (HEM) has been carried over by this framework. The mechanical characteristics of the HEM are clarified, and the comparisons of the numerical results with experimental results are demonstrated to show the reliability of this platform. The advanced morphing and automatic partial domain remeshing schemes are employed for the fluidstructure coupling simulation of a flapping wing structure in a water channel. These features have been fully implemented in a general purposed multi-physics simulation and design optimization software named INTESIM. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction The fluid-structure interaction (FSI) analysis has made remarkable progress in recent years [1–6]. The key features for simulating a fluid-structure interaction analysis consist of the coupling method for handling the coupled problems at different level of difficulties; mesh moving or mesh adapting technique for non-structural (fluid) field with large domain change; the stabilization method, auto-time stepping and bisection scheme to improve the stability and efficiency of the multi-physics simulation; and parallel computing technology to deal with the large scale problems. As mentioned above, in FSI analysis, the method used to satisfy the geometrical compatibility and the equilibrium conditions on the interface between the fluid and the structure is one of the key factors. The strong coupling method (direct method) and the weak coupling method (iterative method) have been mainly used to solve the same fundamental continuum and discretized equations. For the strong coupling method, also referred to by some other authors as direct coupling or simultaneous solution, the coupled equations are solved directly and the variable vectors of the coupled system are updated simultaneously. Zhang and Hisada [4,5] have investigated the strong coupling and weak coupling ⇑ Corresponding author. Tel.: +86 10 62796797; fax: +86 10 62770209. E-mail address: [email protected] (B. Zhu). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.02.019

methods in FSI analysis and used the strong coupling method with stabilization schemes for the numerical analysis of an artificial heart. As indicated by Zhang and Hisada [4,5], the convergence rate of the strong coupling method is of the second order and for some FSI problems with a strong dependence between the fluid and structure, a strong coupling method with appropriate stabilizations is required. However, the use of strong coupling is limited because there are some drawbacks in the strong coupling method: (1) the need for modification of the existing fluid and structure solver; (2) less flexibility between the fluid dynamic (FD) solve and structure dynamic (SD) solver in the time integration and mesh discretization, and (3) requirements of more storage and computational time for one iteration. There are a few limitations in Zhang’s technology [4,5] to simulating the large scale highly nonlinear fluid-structure problem with large domain changes: (1) The mesh morphing method without automatic remeshing make it difficult to handle the FSI problems with extremely large boundary movement, (2) the strong coupling method is developed for highly coupled problems but the interface meshes need to be identical across fluid-structure interface, this will affect the modeling convenience and the computational efficiency, (3) no parallel capability to solve large scale problems in an efficient way. There is some general purpose commercial CAE software featuring fluid-structure interaction capability, e.g., ANSYS multi-physics

37

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

(MFX ANSYS-CFX), ABAQUS with StarCD or Fluent through MpCCI, Algor, ADINA-FSI, etc. Most of these FSI products limit to weak coupling method, this make it impossible to solve strongly coupled highly nonlinear fluid-structure interaction problems [4,5]. ADINA software has both strong coupling and weak coupling capacity and has been used for solving some FSI problems [1,6,7], although it is not suitable for solving large scale FSI problem due to the lack of high performance parallel computing algorithm. In this research we will provide a systemic coupled physics framework to simulate large-scaled strongly coupled fluid-structure interaction problem with high nonlinearity and extreme large fluid domain changes in a robust, efficient, accurate, and convenient way. This framework consists of strong interface coupling methods that can deal with not only the compatible but also incompatible mesh discretization across physics interfaces. It also provides an advanced morphing technology to move the interior nodes in fluid domain according to the interface movement. Laplace equation or elasticity equation with adjustable mesh stiffness can be employed for solving this moving boundary problem, the mesh stiffness of the fluid element can be controlled to achieve the best mesh quantity. The morphing scheme for improving the mesh quality without changing topology of the mesh is limited to certain level of mesh deformation. For problems with extreme large domain change automatic remeshing technology in which a new mesh is recreated may be needed. In this research, we combine advanced morphing and automatic remeshing to deal with the domain changes from small to extreme large. Numerical stabilization method in both time and spatial domain, and auto-time step and bisection algorithms are employed to improve the stability and efficiency of simulation for the strongly coupled FSI problems. To handle the large scale problem as well as the ill-conditioned global matrix properties of the coupled system, a MPI based parallel direct sparse solver is used. For element matrices and vectors calculations, OpenMP based shared memory parallel scheme is employed. All of those features listed above have been implemented in a general purposed multi-physics simulation and optimization software INTESIM [8,9]; it has a very convenient and consistent graphical user interface for all physics models. To demonstrate the effectiveness of the proposed framework, FSI simulation of a hydraulic engine mount has been carried over. The hydraulic engine mount has been developed and used for isolating vibration and noise from the engine to the body of the vehicle for more than 30 years. The mechanical behavior of HEM under the excitation of high frequency oscillating forces includes complex fluid–solid interaction [10]. There are a lot of challenges in numerical simulations: the large domain changes, strong coupling between the hyper-elastic rubber and the viscous fluid flow, the high frequency vibration, and large scale problem. The numerical simulations will be very useful for clarifying the mechanical characteristics and doing optimum design of HEM, the simulation capabilities provided here are able to solve these numerical difficulties and successfully simulate the dynamic responses of HEM in three dimensions and in a fully coupling way. The difficulties for the fluid-structure interaction simulations of a flapping wing structure not only comes from the coupling between the flow field and vibrations of the wing structure but also from the extreme large fluid domain changes caused by the wing’s motions. The fully coupled simulation of the falling wing structure in a water channel has been carried out. It shows the effectiveness of the proposed morphing and remeshing to deal with the coupling problem with arbitrary and extreme domain changes. In Section 2 we will briefly introduce the ALE formulation for fluid flow and the FE equation for nonlinear hyper-elasticity materials that will be used to model the rubber material of Engine Mount. As one of the major topics, the technical details of strong coupling methods will be given in Section 3. Section 4 will discuss

advanced mesh morphing and automatic remeshing method for handling extreme large CFD domain changes. The numerical schemes to improve the robustness and efficiency in the simulation of highly nonlinear coupled problems are introduced in Section 5. The hybrid parallel computing technique for coupled system will be elaborated in Section 6. In Section 7 two examples of FSI problems with simulation modeling and comparisons with experimental results will be provided. Finally the concluding and remarks will be given in Section 8. 2. Fundamental equations 2.1. Arbitrary Lagrangian-Eulerian formulation for barotropic fluid In the ALE description, the mesh is allowed to move arbitrarily, i.e., it does not necessarily coincide with particle motion. Therefore, the application of the ALE formulation to fluid flows makes it possible to directly couple a fluid with a solid, the interface of which moves according to the dynamic equilibrium of the whole system. The material time derivative of any variable is expressed as:

  @ðÞ @ðÞ ¼ þ ðv  v m Þ  rðÞ @t X @t v

ð1Þ

where vm is the velocity of the mesh, o()/ot|X the material time derivative, and o()/ot|v the time derivative in the ALE coordinates, the velocity of which is vm. We assume the viscous fluid to be isothermal and barotropic (i.e., F(P, q) = 0) and that oP/oq = B/q with B, P and q being the fluid bulk modulus, the pressure, and the fluid density, respectively. By applying Eq. (1) to the Eulerian form of Navier–Stokes equations one obtains the corresponding ALE form [11]:

 1 @p 1 @p @ v i þ ci þ ¼0 B @t v B @xi @xi

q

 @ v i  @ v i @ rij þ qc j ¼ þ qg i @t v @xj @xj

in Rft

ð2Þ

in Rft

ð3Þ

where r is the Cauchy stress tensor, g the acceleration of gravity, and c = v  vm. Rft denotes the spatial domain bounded by the boundary @Rft of interest at any instant t. Here the superscript f stands for the fluid component. The fluid is assumed to be Newtonian, and the constitutive equation is:

1 2

rij ¼ pdij þ lðv i;j þ v j;i Þ

ð4Þ

where d is the identity tensor and l is the dynamic viscosity of the fluid. The boundary @Rft is composed of @Rgt and @Rht corresponding to Dirichlet- and Neumann-type boundary conditions, respectively.

v i ¼ gi

on @Rgt

ð5Þ

rji nj ¼ hi

on @Rht

ð6Þ

2.2. Finite element equations By using the Galerkin method and finite element discretization for the ALE form Navier–Stokes Eqs. (2) and (3), we can obtain the FE Eq. [4]. 

Mf uf þCf uf ¼ Ff

ð7Þ

where, " Mf ¼

" # # GT KP MP 0 ; ; Cf ¼ G Kl þ K 0 M

(

uf ¼

Pf Vf

) ; Ff ¼



0 Ffv



ð8Þ

38

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

MP and M are the generalized mass matrices for pressure and velocity, respectively, and KP and K are the generalized matrices of convective terms for pressure and velocity, respectively. Kl is the fluid viscosity matrix, G is the divergence operator matrix, Pf and Vf are the pressure and velocity vectors, respectively, and Ff is the external force vector.  indicates the time derivative in the ALE coordinates. In the present work the secant matrices are used for nonlinear iterations, so the incremental forms of Eq. (7) can be simply obtained as shown later.

3. Coupling methods 3.1. Interface conditions For a fully coupled fluid solid system the geometrical compatibility conditions and the equilibrium conditions on the interface must be satisfied. In this research the no slip conditions are assumed on the fluid and solid interface, i.e., the compatibility conditions:

v fi

¼ v si

ði ¼ 1; 2; 3Þ on @Rct

2.3. Incompressible hyper-elastic material for the rubber

and the equilibrium conditions:

2.3.1. High Mooney–Rivlin model The hyper-elastic material satisfies the following equation [12]:

rfji nfj þ rsji nsj ¼ 0 ði ¼ 1; 2; 3Þ on @Rct

Sij ¼

@W @Eij

ð9Þ

where E is the Green–Lagrange tensor, S the 2nd Piola–Kirchhoff stress tensor, and W the elastic potential function. The high Mooney–Rivlin model in Reduced Invariants form can be expressed as:

~ C  3Þ þ c3 ð~IC  3Þ2 þ c4 ð~IC  3ÞðII ~ C  3Þ W HR ¼ c1 ð~IC  3Þ þ c2 ðII 2 3 2 ~ ~ ~ ~ þ c5 ðIIC  3Þ þ c6 ðIC  3Þ þ c7 ðIC  3Þ ðIIC  3Þ ~ C  3Þ2 þ c9 ðII ~ C  3Þ3 þ c8 ð~IC  3ÞðII ~IC ¼ IC ; 1=3 IIIC

ð10Þ

~ C ¼ IC II 2=3 IIIC

ð11Þ

where IC IIC and IIIC are the main constants of the right Cauchy– Green strain tensor C(E = (C  I/2)). The parameters of c1  c9 can be determined by experiments. We use the mixed finite element method for the spatial discretization of the hyper-elastic solid.

€ s þ Q s ðU s ; P s Þ ¼ F s Ms d

ð12Þ

3.2. Strong coupling methods We can solve the coupled system by a strong coupling method, in which the coupled equation is solved implicitly, it means the variable vector of the coupled system is obtained simultaneously by solving the linearized coupled equation, and also the variable vectors of the coupled system are updated simultaneously. 3.2.1. Direct interface coupling method If mesh discretizations of the fluid and solid on the interfaces are compatible i.e., the fluid and solid mesh share the same interface nodes and the underline elements (of element faces) on the interface (Fig. 1), then the direct interface coupling can be used to satisfy the interface coupling conditions. Using Lagrangian description for both of the fluid and the solid on the interface, then the variable vector of the coupled system can be divided into three parts,



8 f 9 > > < ui = fs

uc > : s > ;

;

ui

8 9 > < > = fs s d ¼ dc > : s > ; di

ð17Þ

where ufi ðufi ¼ fPf ; Vfi gT Þ is the pressure and the internal velocity vector of fluid, ufsc ðufsc ¼ Vfsc Þ is the coupled velocity vector, 

usi usi ¼ ðfVsi ; Ps gT Þ, the internal velocity vector and the time derivafs

where,

Ms ¼

ð16Þ

are given, where n is the outward normal. Superscripts f and s indicate the value corresponding to the fluid and the solid, respectively, and the @Rft interface of the fluid and the solid.

ufs ¼ 2.3.2. Finite element equation The equilibrium equations of solid in matrix form can be expressed as follows:

ð15Þ

tives of the pressure vector of the solid domain. dc ðUfsc Þ is the coupled



Msv 0

;

s

d ¼



Us Ps

 ;

( Qs ¼

Q sv Q sp

) ;

Fs ¼



Fsv

s s di ðdi



0 ð13Þ

fUsi ; Ps gT Þ,

¼ the internal displacement displacement vector, vector and pressure vector of solid. Based on the definition of the variable vector of the coupled system, and combining Eqs. (7) and (12), the coupled equation can be expressed as:

Ms is the mass matrix, Us and Ps are the displacement and pressure vectors, respectively, and Qs and Fs are the internal and external force vectors, respectively. Subscripts v and p represent the value corresponding to velocity and pressure, respectively. .. indicates the second order time derivative in the Lagrangian coordinates. In the present work the tangent matrices are used for nonlinear iterations [12].

€ s þ t Ks Dds ¼ tþDt Fs  t Q s M s Dd

ð14Þ

where tKs is the tangent stiffness matrix of solid at time t, and D stands for the increments of variable vector of solid. Hyper-elasticity material is used for the rubber modeling in the hydraulic engine mount analysis in Section 7.1. The linear elasticity material with geometrical nonlinearity is employed for the large deflection analysis of the wing structure in the Section 7.2.

Fig. 1. Direct interface DOF coupling.

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

Q fs ¼ F fs

39

ð18Þ

where the internal force vector Qfs includes the inertial force of the coupled system, fs

s

Q fs ¼ M fs u þC f ufs þ Q s ðd Þ

ð19Þ

Ffs is the external force vector. In Eq. (19),  denotes the time derivative in ALE coordinate for fluid and time derivative in Lagrangian coordinate for interface and solid, respectively. Mfs is the mass matrix composed of those for the fluid and the solid, Cf consists of the divergence, viscous and convective terms of fluid. The details of these coupled matrices and vectors can be expressed as follows.

2

MP 6 6 0 6 M fs ¼ 6 6 0 6 4 0 2

0 M fii

0 M fic

0 0

M fci

M fcc þ M scc

M sci

0

M sic

M sii

0

0

0

0

KP 6 6 Gi 6 Cf ¼ 6 6 Gc 6 4 0

3 0 7 07 7 07 7; 7 05

Fig. 2. Schematic of a mapping technology.

0

3

GTi

GTc

Kii þ K lii

Kic þ Klic

Kci þ K lci

Kcc þ K lcc

0

0

7 0 07 7 0 07 7 7 0 05

0

0

0 0

0 8 9  > > > > > > > > > > > < > = s s Q ¼ Qc ; > > > > Qs > > > i > > > > : s> ; QP

0 0

8 9 0 > > > > > > Ff > > > > > < i > = fs fs F ¼ Fc > > > > > > > F si > > > > : > ; 0

ð20Þ

ð21Þ

In the above equations, the subscripts c and i denote the coupling parts and the internal parts of variables, respectively. The incremental form of the coupled Eq. (18) can be obtained by taking linearization. t

fs

s

M fs D u þt C f Dufs þ t K s Dd ¼ tþDt F fs  t Q fs

ð22Þ

where tKs is the tangent stiffness matrix and the coefficient matrices t fs M and tCf are secant matrices. Although the linearization of tMfs and tCf has been discussed by Zhang [13], in practice, if an effective time step control algorithm [13] is used to control the time step Dt, the secant matrices of tMfs and tCf frequently work quite well. Using implicit time integration method (Newmark-b) for the linearized coupled Eq. (22), one can obtains the linear Eq. (23): t

fs

t

f

2t

s

fs

½ M þ cDt C þ bDt K D u ¼

tþDt

t

F Q

fs

ð23Þ

where c and b are the parameters thatfs can be chosen so as to obtain  integration accuracy and stability, D u is the increment of variable vector of the coupled system. Furthermore Eq. (23) can be expressed as the simplest form Eq. (24). fs

Kfs D u ¼ Rfs

ð24Þ

fs

where K is the effective mass matrix of the coupled system, Rfs is the residual vector of the coupled system. 3.2.2. Multipoint constraint equation based strong coupling When different physics models occupy same geometry location (e.g., vertex, edge, face, solid) but not sharing the same interface geometry, and nodes and the underline elements has co-related DOFs, a constraints equation based coupling can be employed (Fig. 2). By using mapping technology, un-matching mesh discretization is allowed across the coupling interface.

In the case the fluid mesh and structure mesh on the coupling interface are incompatible, appropriate searching tools need to be used in order to find the correlations between the fluid nodes and the structure nodes. An efficient bucket search based mapping algorithm is employed for mapping the fluid nodes to the solid mesh and vice versa [14]. In this research we treat the fluid nodes on the interface as slave nodes, on the contrast the structure nodes on the interface as master nodes and forced the displacement, velocity and acceleration of the fluid nodes on the interface equal to the corresponding value of structure. The constraint equations and the force distributions of the fluid forces to structure are based on the mapping results. The constraints for velocity vector are given in Eq. (25).

vf

¼

n X

W i v si

ð25Þ

i¼1

where vf is the velocity vector of a fluid node on the interface, n the number of related master/structure node, Wi the weight on structure node i and v fi is the velocity vector of structure node i. Using multipoint constraint method one can get similar coupled equations as in direct interface coupling method. 4. Advanced mesh control technique 4.1. Advanced morphing method 4.1.1. Boundary conditions for mesh motion For the fluid–solid interaction analysis the boundaries for controlling the mesh motions of fluid consist of rigid walls, open boundary, and the interfaces with deformable solid. Under these boundary conditions, the following methods are employed for mesh control in the fluid domain. The mesh is fixed in space for rigid walls and open boundaries, i.e.,

um ¼ 0

ð26Þ

is used. On the interface with a deformable solid, fluid nodes are attached to the solid, i.e., the Lagrangian motion is assumed as:

um ¼ us

ð27Þ

where us is the displacement vector of the structure. For the incompatible mesh interface in the present research, an appropriate mapping and interpolation described is used to get appropriate displacement boundaries for fluid domain.

40

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

4.1.2. Morphing equations for solving the mesh motion of interior nodes For controlling the motions of interior nodes, mesh distortion must be effectively avoided. Because the geometrical configuration of the fluid domain can be complex and the deformation of the fluid domain can be large, a Dirichlet-type boundary value problem is usually solved for determination of nodal displacements in each direction (Eq. (28))

r  ðjrum Þ ¼ 0

in Rft

ð28Þ

where um is the mesh displacement, j the element mesh stiffness, generally as a function of element size, level of element distortion, or the distance from the wall. We apply Eq. (28) to each direction of the simulation domain; it means the mesh deformations in each direction are uncoupled. 4.2. Automatic remeshing technology Although the automatic mesh updating algorithm (also called morphing) does not change the element topology, however it makes difficult to handle the FSI problems with extremely large boundary movement. Therefore in this case remeshing (generating a new mesh) is helpful to maintain the element quality and to avoid morphing failure. 4.2.1. Mesh quality measurement In this research we use three simple ways to measure element qualities to decide when remeshing is necessary [15]. 1. Generalized element aspect ratio AR. For 3D tetrahedron element:

AR ¼ ððlav g Þ3 =VÞ=8:48

ð29Þ

2. Changes of element size. For 3D tetrahedron element:

    VðtÞ  VOCH ¼ exp log Vð0Þ 

ð30Þ

3. Changes of element aspect ratio:

    ARðtÞ  ARCH ¼ exp log ARð0Þ 

ð31Þ

Fig. 4. Inertial track.

where lavg is the average length of the element edges, V the element volume, V(t) the volume at time t, and V(0) is the volume at time 0. The number of 8.48 in Eq. (29) is the aspect ratio of a regular tetrahedron element. This framework provides a convenient way to specify criteria for these quantities, if any element quality reaches the specified value, remeshing will be considered for that element. These mesh quality measurements are easy to implement and suitable for evaluating the qualities of triangle mesh and tetrahedral mesh. 4.2.2. Remeshing tool Based on the exterior surface mesh and the total displacement at the exterior nodes of a selected element group, our remeshing tool uses the Delaunay triangular method [16] to create a new mesh and then do swapping, coarsening, refining, smoothing etc. to improve the mesh qualities. This can guarantee a new mesh with more desirable element qualities than the one which uses mesh improvement operations only. 4.2.3. Result mapping and interpolation Interpolating the nodal or element based values (include boundary conditions) from the old mesh (before remeshing) to the new mesh (after remeshing) is the last step in remeshing process. We use the bucket search method [17] to search the mapping between the old mesh and the new mesh. It loops over all elements on the old mesh for each new node or new element. Element local coordinates based linear interpolation will be used to get the nodal solution for the interior nodes. Because there is no topology change on the boundary nodes of the remeshing domain, the nodal values on the boundary of the new mesh can be obtained directly from the coinciding node of the old mesh. We also update the element based FSI boundary conditions on the new fluid mesh. 5. Numerical schemes for the highly nonlinear multi-physics problems For a highly nonlinear coupled system, numerical stabilization schemes are needed to improve the stability and computational efficiency during time integration and nonlinear iteration process. The typical stabilization methods used in this FSI solver are given in this section. 5.1. Stabilization method in time and spatial domain For fluid flow with high convection term the Streamline Upwind Petrov/Galerkin (SUPG) method [18] is used to stabilize the velocity field. Since we used the same shape function and interpolation Table 1 Material properties of the rubbers. Locations

Fig. 3. Structure of HEM.

Upper spring Un-coupler Lower diaphragm

c1 (Pa)

q (kg/m3)

c2 (Pa) 6

0.409016  10 0.241209  10 0.354556  106

4

1.991666  10 0.689460  106 0.0

1.11  103 1.18  103 1.098  103

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

(a) t=0.0

(b) t=T/4

(c) t=T/2

(d) t=3T/4

41

Fig. 5. Deformations and pressure distributions of the solid parts for f = 10 Hz.

(a) t=0.0

(b) t=T/4

(c) t=T/2

(d) t=3T/4

Fig. 6. Velocity distributions of the fluid domain for f = 10 Hz.

for velocity and pressure degree of freedoms in our element, the pressure-Stabilizing/Petrov–Galerkin (PSPG) stabilization method [19] is employed to avoid the pressure oscillation. To avoid high

frequency oscillation in time, numerical damping can be introduced by adjusting the time integration parameter, for details please refer to Bathe [20].

42

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

Fig. 8. Dimensions of the model, front view and top view.

(a) f=5.0Hz

Fig. 9. Motion of the wing structure in one flapping cycle.

(b) f=20.0Hz

Fig. 10. The mesh discretization and specified motion location on the wing structure.

(c) f=50.0Hz Fig. 7. Reflection force on the lower steel boundary (computations vs. experiments).

5.2. Auto-time stepping and bisection scheme For a time transient and nonlinear problem, it is common that the level of nonlinearity changes during the analysis, it means that the appropriate time step size may vary during the simulation. It is

Fig. 11. The near wing CFD meshes with higher element morphing stiffness (1000 times higher) and excluded from remeshing.

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

43

Besides the auto-time stepping, we also provide bisection capability to allow redo a diverged solution with half time step size of the original one. When the time step size reached the minimum value and the number of iteration also reached the maximum value, then the solver will automatic back to the previous time step, reduce the time size by half and start a new time step simulation. This bisection operation will also be performed when negative volume of fluid element occurred, or any other divergence that is caused due to high nonlinearity. 5.3. Relaxation of the solution vector Relaxation method is another well known approach to stabilize a nonlinear simulation procedure. Relax the solution may slow down the convergence rate but it will make the nonlinear iteration more stable. It is helpful for highly nonlinear coupled physics problems. Eq. (33) is used to relax the changes of solution vector in this research.

Fig. 12. CFD region excluding the near wall elements for remeshing.

optimal and sometime necessary to adjust the time step size to achieve a stable and fast solution. The auto-time stepping scheme we used here is shown in the following equations.

Dt nþ1 ¼ Dt n  Nt arg et =Nactual if Dtnþ1 > Dt max then Dtnþ1 ¼ Dt max

ð32Þ

if Dtnþ1 < Dt min then Dtnþ1 ¼ Dt min where Dt is the time step size, Ntarget the user specified target number of iteration, Nactual the actual number of iteration in time step n. If the actual number of iteration is larger than the user specified target number then time step size will be reduced for the next time step, and vice verse, but the time step size will not exceed the user specified range.

uiþ1 ¼ ui þ aDu; 0 < a 6 1:0

ð33Þ

where ui+1 is the solution vector at iteration i + 1, ui the solution vector at iteration i, Du is the incremental of the solution vector from the coupled incremental equation solver. a is the relaxation factor, with a = 1.0 without relaxation, by reducing the a value from one to zero it will make the nonlinear iteration stable but will also slow down the convergence speed. 6. Parallel computing It is necessary to use parallel computing for large scale fluid-structure problem. In the study a hybrid parallel coupling technique, shared memory OpenMP for element calculation and

Fig. 13. Time histories of the fluid forces on the fluid-wing interface.

44

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

MPI based distributed memory for linear equation solver, are used.

scalability, and very easy to implement into a computer program, it applicable to most computer platforms especially on multi-core machines.

6.1. OpenMP based element matrices and vectors calculations 6.2. MPI based distributed memory direct sparse solver For element matrices and vectors calculation of fluid field and structure field we use OpenMP based shared memory parallel scheme to improve the performance. It has almost linear

Since the coupled equations of the fluid-structure system is usually ill-conditioned, e.g., the thin shell structure coupling with

(a) time=0.375s

(b) time=1.375s

(c) time=1.5s Fig. 14. Velocity field of the flow channel in one flapping cycle.

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

45

(d) time=1.625s

(e) time=2.625s

(f) time=3s Fig. 14 (continued)

water, in this case it is hard to get a converged solution by iterative solver or need a considerable amount of iterations to get a converged solution, instead of using parallel iterative solver the distributed parallel direct sparse solver will be more robust and efficient. In this research we use the Massage Passing Interface (MPI) based distributed memory multi-frontal parallel direct sparse solver to solve the coupled fluid-structure equations [21].

7. Numerical analyses 7.1. Strong coupling simulation of a hydraulic engine mount 7.1.1. Physical model of HEM A typical HEM consists of a rubber spring, two fluid chambers, an inertial track and a un-coupler. Fig. 3 shows the cross section of a

46

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

(a) the morphed mesh

(b) the remeshed mesh Fig. 15. The morphed meshes and the meshes after remeshing at time = 1.375 s.

HEM, and the inertial track is shown in Fig. 4. The rubber spring acts as the main support spring. The inertial track plays a role as tuned isolator damper to provide higher damping in the low frequency range. The un-coupler connecting the upper and lower chamber can provide a lower stiffness to isolate the higher frequency vibrations when the inertial track is clogged. The material properties of fluid are as follows: mass density q = 1.08156  103 kg/m3, dynamic viscosity l = 19.23925 Pa s, bulk modulus B = 2.0  107 Pa. The rubbers are assumed to be hyper-elastic materials, the coefficients of the high Mooney–Rivlin model are listed in Table 1. The material properties of the steel (engine connector) are Young’s modulus E = 2.0  1011 Pa, Poisson’s ratio v = 0.3, and the mass density q = 7.8  103 kg/m3. 7.1.2. 3D numerical model and analysis conditions of HEM Because both the deflection and the deformation of the lower steel boundary are small, to minimize the computational costs we subtract the steel parts and use fixed boundaries for the interfaces between the fluid and the steel and the interfaces between the rubber and the steel. We use the Q1Q0 mixed elements for both the solid and fluid domains of the simplified model. 6510 fluid elements and 3354 solid elements are used for the mesh discretization; the total number of nodes is 21,837. To make sure the mesh discretization listed above is enough for the spatial resolution and able to produce accurate reflection forces on the lower steel boundary, a simulation with a finer mesh (the total number of nodes is 53,720) was carried out. The differences of the amplitude of the reflections forces between two models is less than one percent and the phase difference also is negligible, therefore, the coarse mesh has been used for case studies to reduce the computing time. The analysis procedure consists of the static load phase and high frequency forced vibration phase. First, the initial load F0 is

applied slowly (in 1000 s) to the upper steel boundary, then the vibration displacements are applied to the upper steel boundary: d = Asin(2pft), here A is the amplitude of the vibration, f is the frequency and t indicates the time. In this paper, we select the interesting case of F0 = 1200 N, and f = 5, 10, 20, and 50 Hz for the calculations, and the comparisons of the numerical results with the experimental results are presented. As the temporal resolution concerns, in this study 40 time steps for each time cycle at different vibration frequencies are adopted to produce a converged nonlinear solution as well as accurate and smooth solution curves for reflection forces in the time domain. 7.1.3. Results and comparisons The deformations and pressure (p ¼  13 ðr11 þ r22 þ r33 Þ) distributions of the rubber domain in one vibration circle are shown in Fig. 5a–d for f = 10 Hz, and the corresponding velocity distributions of the fluid domain are shown in Fig. 6a–d. Although there is a limitation of demonstrating the behavior of HEM by those figures, we know from the animations at different frequencies that for the high frequency forced vibration the un-coupler works well, some animations can be downloaded from the website [22]. On the other hand the fluid flow in the inertial track provides more damping for low frequency vibration. The comparisons of the time histories of the reflection force on the lower steel boundary from the numerical calculations and experiment measurements are given in Fig. 7a–c for f = 5 Hz, f = 20 Hz, and f = 50 Hz. Fig. 7 demonstrates that the numerical results have a good agreement with those of experimental ones, and the differences are between 5% and 15% on average. The differences are getting larger with higher frequency, which is probably due to the material damping of the rubber materials having more effect at higher frequency. From Fig. 7c, we also find an oscillation in the solution in the first loading cycle, higher amplitude at the first half cycle and lower in the

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

second half. This oscillation maybe induced by the inertial effect, which is important at high frequency. After three time cycles the initial oscillation effect is damped out and the solution reaches a stable periodic state. 7.2. Fluid structure simulation of flapping wing structure in a water channel 7.2.1. Model of flapping wing in a water channel This is an example of a 3D wing structure that flaps in a water channel.ThedimensionsofthemodelareshowninFig.8(frontviewandthe topview).Forthefluiddomain,theinletvelocityof1.0 m/sisspecified, an opening boundary with zero pressure is applied on the outlet, and fixedwallboundariesareappliedonallofthewallboundaries.Themotions of the wingstructure in one flapping cycle are dramaticallydemonstrated in Fig. 9. In this simulation a forced motion is applied on the surrounding nodes of the rotating axis (red1 colored nodes through rotating axis in Fig. 10), and 4276 tetrahedral elements are used to discretize the wing. The motions of the wing in one flapping cycle consists of clockwise rotation from initial position with zero pitching angle to a position with pitching angle of 30° at 0.375 s, in the second phase from 0.375 s to 1.0 s the wing translates vertically in the channel, then from 1.0 s to 1.625 s the wing structure does counter clockwise rotation that makes it back to zero pitching angle and at the top position of the channel. By contrast, in the second half of one flapping cycle the wing will do clockwise rotation, vertical translational motion, and counter clockwise rotation back to original position, sequentially. Due to the large motion of the wing structure to avoid the mesh distortion, advanced morphing scheme with automatic partial domain remeshing technique is needed. In the morphing controls, sliding mesh moving conditions are applied on the inlet, outlet and all wall boundaries. In order to maintain the mesh size and to minimize the mesh deformation in the near wall region, we use varying mesh morphing stiffness for different analysis domain, i.e., a 1000 times larger mesh morphing stiffness is specified in the near wall region (Fig. 11) than the far away regions from the wing (Fig. 12, the fluid elements excluding the elements in the near wall box); we also exclude those near wall elements from remeshing process. The numbers of the fluid tetrahedral elements in the near wall region (Fig. 11) and the far away regions from the wing (Fig. 12) are 23,774 and 36,516, respectively. The material properties of fluid are as follows: mass density q = 1000 kg/m3, dynamic viscosity l = 0.001 Pa s, bulk modulus B = 1.0  106 Pa. The material properties of the flexible wing structure in the fluid-structure interaction simulation are Young’s modulus E = 7.0  1010 Pa, Poisson’s ratio v = 0.3, and the mass density q = 2.7  103 kg/m3. Due to the high Reynolds number flow, the LES Smagrinsky model [4,23] is employed in this study case to model the turbulence behavior of the fluid flow.

7.2.2. Simulation results The time histories of the integrated fluid force on the coupled fluid-wing interface in FSI analysis are plotted on Fig. 13a, by contrast, Fig. 13b provides the time history of the fluid forces on the moving boundary from a pure CFD moving boundary simulation, and it assumes a rigid body motion of the wing structure. It shows obviously from Fig. 13 that due to the flexibility of the wing structure in the coupled fluid-structure interaction, the maximum value of the fluid forces in the X and Y direction are different and slightly larger than those of the pure CFD moving boundary problem, also 1 For interpretation of color in Fig. 10, the reader is referred to the web version of this article.

47

due to the high stiffness of the wing structure the pattern of the time history curves from the two simulations are quite similar. The motion and velocity distribution of the flow field at different time points in one flapping cycle are demonstrated in Fig. 14. The deformed meshes after the morphing process and the meshes after remeshing at time 1.375 s are given in Fig. 15. This example shows that the advanced morphing and partial domain remeshing scheme proposed in this research can handle the periodic extremely large domain changes well. More investigations about the mesh resolution on the boundary layer, detailed discussions about the accuracy of the simulation results and comparisons with experimental results exceed the scope of this paper, and we will provide more information on this subject in coming studies. 8. Concluding remarks In order to solve the complex fluid–solid interactions, a numerical framework is proposed, in which ALE finite element formulation is employed for the fluid, and total Lagrangian formulation for the structure. The full interactions between fluid and structure are considered by the strong coupling method that allows incompatible mesh discretization across physics. An advanced Laplace equation based mesh morphing with automatic remeshing is used for treating the fluid domain changes from small to extremely large. OpenMP based shared memory parallel for element calculation and distributed sparse solver for solving linearized equation make this framework capable to simulate large scale coupled problems. In addition, stabilization schemes in time and spatial domain, auto-time stepping, and bisection schemes guarantee this framework more reliable and efficient for highly nonlinear coupled fluid-structure problems. The developed numerical schemes were successfully used for the numerical analysis of a hydraulic engine mount. The mechanical characteristics of HEM were clarified and the numerical results have a good agreement with those of the experimental results. The flapping wing structure in a fluid channel with complex motions were also studied, it tells us that the proposed advanced mesh morphing scheme with automatic partial domain remeshing can deal with the fluid-structure interaction problems with extreme large domain changes. The present results show great application opportunities of the advanced coupled physics simulation framework in Automotive Engineering and Aerospace Engineering. All of these features have been implemented in a general purposed multi-physics simulation and optimization software INTESIM. Acknowledgements The project is supported by the National Science Foundation of China Grant (Nos. 10872113 and 51179090) and research fund of State Key Laboratory of Hydroscience and Engineering, Tsinghua University (No. 2009T3). The example of HEM used in the paper was from the State Key Laboratory of Automotive Safety and Energy in Tsinghua University, China, thanks to Professor Zhenhua Lu and Dr. Lirong Wang for providing the model and experimental results. The model of fluid structure simulation of flapping wing structure in a water channel is prepared by Mr. Yepeng Han, project manager of INTESIM (Dalian) Co., Ldt., and he also carried out intensive simulations and results analysis for this study case. References [1] Bathe KJ, Zhang H, Ji S. Finite element analysis of fluid flows fully coupled with structural interactions. Comput Struct 1999;72(1-3):1–16. [2] Stein K, Benney R, Kalro V, Tezduyar T, Leonard J, Accorsi M. Parachute fluidstructure interactions: 3-D computation. Comput Methods Appl Mech Eng 2000;190:373–86.

48

Q. Zhang, B. Zhu / Computers & Fluids 60 (2012) 36–48

[3] Taylor CA, Hughes TJR, Zarins CK. Finite element modeling of blood flow in arteries. Comput Methods Appl Mech Eng 1998;158:158–96. [4] Zhang Q, Hisada T. Analysis of fluid-structure interaction problems with structural buckling and large domain changes by ALE finite element method. Comput Methods Appl Mech Eng 2001;190:6341–57. [5] Zhang Q, Hisada T. Studies of the strong coupling and weak coupling methods in FSI analysis. Int J Numer Methods Eng 2004;60:2013–29. [6] Rugongyi S, Bathe KJ. On the finite element analysis of fluid flows fully coupled with structural interactions. Comput Models Methods Appl Sci 2001;2:195–212. [7] Shangguan WB, Lu ZH. Experimental study and simulation of a hydraulic engine mount with fully coupled fluid structure interaction finite element analysis model. Comput Struct 2004;82(22):1751–71. [8] INTESIM2010 Users’ Manual, Dalian, China; 2010 . [9] INTESIM2010 Theory Manual, Dalian, China; 2010 . [10] Geisberger A, Hajepour AK, Colnaraghi F. Nonlinear modeling of hydraulic mounts: theory and experiment. J Sound Vib 2002;249(2):371–97. [11] Huerta A, Liu WK. Viscous flow with large free surface motion. Comput Methods Appl Mech Eng 1988;69:277–324. [12] Watanabe H. Study of mixed finite element analysis for incompressible hyperelastic materials. PhD thesis, University of Tokyo; 1995. [13] Zhang Q. Analysis of structure-fluid interaction problems with structural buckling and large domain changes by ALE finite element method. PhD thesis, University of Tokyo; 1999.

[14] Jansen K, Shakib F, Hughes TJR. Fast projection algorithm for unstructured meshes, computational nonlinear mechanics in aerospace engineering; 1992 [ISBN: 1563470446 chapter 5]. [15] Johnson AA, Tezduyar TE. Simulation of multiple spheres falling in a liquidfilled tube. Comput Methods Appl Mech Eng 1996;134:351–73. [16] Yagawa Y, Yoshimura S, Nakao K. Automatic mesh generation of complex geometries based on fuzzy knowledge processing and computational geometry. Intr Comput – Aid Eng 1995;2:265–82. [17] Jansen K, Shakib F, Hughes TJR. Fast projection algorithms for unstructured meshes. American Institute of Aeronautics and Astronautics, Inc.; 92 [chapter 6] [18] Brooks AN, Hughes TJR. Stream upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equation. Comput Methods Appl Mech Eng 1982;32:199–259. [19] Tezduyar TE, Sathe S. Stabilization parameters in SUPG and PSPG formulations. J Comput Appl Mech 2003;4(1):71–88. [20] Bathe KJ. Finite element procedures. Prentice Hall; 1996. [21] Amestoy PR, Duff IS, L’Excellent L-Y. Multifrontal parallel distributed symmetric and unsymmetric solvers. Comput Methods Appl Mech Eng 2000;184:501–20. [22] http://www.intesim.cn/business.php?id=743. [23] Ghosal S. An analysis of numerical errors in large-eddy simulations of turbulence. J Comput Phys 1996;125:187–206.