Mixed finite element–computational fluid dynamics method for nonlinear fluid–solid interactions

Mixed finite element–computational fluid dynamics method for nonlinear fluid–solid interactions

Chapter 11 Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter Outline 11.1 Updated Lagrangian fo...

6MB Sizes 1 Downloads 59 Views

Chapter 11

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter Outline 11.1 Updated Lagrangian formulation in finite element methods 11.1.1 Principle of virtual work and dynamic equilibrium equations 11.1.2 Expression of stress virtual work 11.1.3 Total Lagrangian formulation in solid dynamics 11.1.4 Updated Lagrangian formulation in solid dynamics 11.1.5 Solution of nonlinear equations 11.1.6 A one-dimensional example 11.2 Updated arbitrary LagrangianEulerian formulations in computational fluid dynamics 11.2.1 History and development of arbitrary LagrangianEulerian description 11.2.2 Basic discretization techniques of computational fluid dynamics 11.2.3 Consistency, stability, and convergence of numerical schemes 11.2.4 Updated arbitrary LagrangianEulerian formulations in finite difference method and finite volume method

410 410 411 412 414 416 418 423 423 424 434

11.3 Mixed finite elementcomputational fluid dynamics solutions for nonlinear fluidsolid interaction problems 11.3.1 Direct or simultaneous integration 11.3.2 Partitioned iteration 11.3.3 A simple example 11.4 Numerical examples by the mixed finite elementfinite difference method 11.4.1 Description of nonlinear rigid bodywater fluidsolid interaction systems 11.4.2 Arbitrary LagrangianEulerian description of the water motion 11.4.3 Numerical formulations 11.4.4 A fluidmassspring interaction system 11.4.5 A spring-supported fluidrigid-dam interaction system 11.4.6 Two-dimensional rigid body floating in a free-surface fluid 11.4.7 Prescribed motion of a rigid cylinder floating in the water 11.4.8 Flows around a bluff body

449 449 450 451 457 458 462 463 467 469 471 474 475

443

The application of finite element (FE) methods (FEM) in solid structural analysis (see, e.g., Bathe, 1982; Zienkiewicz and Taylor, 1989, 1991) and computational fluid dynamics (CFD) to describe fluid flows (see, e.g., Hirsch, 1988, 1990; Anderson, 1995; Chung, 2002) have achieved increased sophistication and success in recent years. How to combine this well-developed FEM for structures and the powerful CFD approaches for fluids in order to deal with nonlinear fluidsolid interaction (FSI) problems is discussed in this chapter. We call this numerical approach simulating nonlinear fluidstructure interaction problems as a mixed FECFD method. The benefit of this mixed method is effectively using the available powerful commercial software for FEM and CFD to simulate complex nonlinear FSI problems in engineering. Generally, we assume that the structure undergoes large rigid motions as well as large elastic deformation, while the fluid flow is governed by nonlinear, viscous, or nonviscous field equations with nonlinear boundary conditions applied to the free surface and FSI interfaces. As we have mentioned, for numerical simulations of nonlinear fluidstructure interaction problems, when powerful FEM methods for solids and CFD methods for fluids are used to solve fluidstructure interaction problems, a difficulty arises. In solid mechanics, a Lagrangian description usually adopts a material point formulation of the solid structure, whereas in fluid dynamics, a Eulerian description describes the behavior of the fluid at each fixed point in space. Thus variables describing a structural motion are functions of the material coordinates fixed to each particle in the solid and time, whereas variables describing a fluid flow are functions of the coordinates fixed in space and time. When the structure moves, the material coordinates also move from their original FluidSolid Interaction Dynamics. DOI: https://doi.org/10.1016/B978-0-12-819352-5.00011-2 © 2019 Higher Education Press. Published by Elsevier Inc. All rights reserved.

409

410

FluidSolid Interaction Dynamics

positions to new positions in space, but the Eulerian coordinates remain unchanged, so that a separation of the mesh points on the solid from those in the fluid, which initially coincided before the system was disturbed, is produced. For linear small disturbance problems, as discussed in Chapters 47, this separation is neglected and the original statically equilibrium configuration of the FSI system forms the basis of the mathematical model on which a numerical analysis is developed. However, for nonlinear problems, the difference between Lagrangian and Eulerian descriptions of a motion has to be fully described to keep the validity of compatibility conditions on the coupling interfaces. The main way to overcome this difficulty is to replace the original configuration of the system by its configuration at time t as the reference position in order to investigate the motions of the system from time t to t 1 Δt. Based on this reference configuration, for the solid structure, its position at time t is chosen as its material coordinates to describe its motion, while the fluid spatial mesh, with a suitable moving velocity, is structured to derive fluid formulations, which are, respectively, called the updated Lagrangian (UL) description in FEM analysis of solids and updated arbitrary LagrangianEulerian (ALE) mesh description in CFD, as discussed in Chapter 3, Fundamentals of continuum mechanics. For readers effectively to understand this mixed FECFD method, we will discuss the UL formulation in FEM and the updated ALE formulation in CFD before tackling the combined FECFD approach for nonlinear FSI problems. The detailed numerical implication process and application examples will be presented after the fundamental theory is described.

11.1

Updated Lagrangian formulation in finite element methods

The various theories, formulations, and applications of FEM in solid structural analysis have been well reported (see, e.g., Bathe and Wilson, 1976; Bathe, 1982, 1996; Zienkiewicz and Taylor, 1989, 1991). We assume that readers have known the fundamental concepts of FEM, as explained by a simple example presented in Section 1.4.4, so that we will not give more details on this process but only discuss the elementary knowledge on the UL formulation in FEM, which closely involves the mixed FECFD method investigated in this chapter. Readers who are not familiar with FEM may read one of the references previously listed.

11.1.1

Principle of virtual work and dynamic equilibrium equations

The dynamic equilibrium equation of a nonlinear dynamic system at the current time t 1 Δt can be derived by using the d’Alembert’s principle and the principle of virtual work as follows. The d’Alembert principle assumes that an imaged virtual inertial body force 2ρS u€i (a negative multiplication of mass density and acceleration) is applied to the dynamic system, so that the system can be investigated for its dynamic equilibrium, like a static system. The principle of virtual work states that a system is in its static equilibrium state if and only if its total virtual work, done by the all internal and external forces applied in/on the system through virtual displacements, vanishes. A virtual displacement is one of admissible displacements discussed in Section 1.4.3.1, which is a continuous derivative function satisfying the displacement boundary condition. This implies that a virtual displacement vanishes on the displacement boundary where the displacement is prescribed. Based on these two principles, the equilibrium of the body at the current time t 1 Δt requires that ð   t 1 Δt σ ij δt1Δt eij 2 t 1 Δt ρS u€i δui dt 1 Δt Ω 5 t 1 Δt W ; (11.1) t1Δt

Ω

where, as shown in Fig. 11.1, we define that ui is the displacement vector on the current configuration at time t 1 Δt, a function of t 1 Δt xi and time t 1 Δt; u€i is the acceleration vector on the current configuration at time t 1 Δt, a function of t 1 Δt xi and time t 1 Δt; δui is the virtual displacement vector imposed on the current configuration at time t 1 Δt, a function of t 1 Δt xi and time t 1 Δt; t 1 Δt xi is the Cartesian coordinates of a material point at time t 1 Δt; t 1 Δt Ω is the volume of the body at time t 1 Δt; t 1 Δt σi is the Cauchy stress tensor at time t 1 Δt, that is, the forces per unit areas in   the deformed geometry; δt 1 Δt e ij 5 ð1=2Þ ð@δui =@t 1 Δt xj Þ 1 ð@δuj =@t 1 Δt xi Þ , strain tensor corresponding to virtual displacement at time t 1 Δt; t 1 Δt

ð W5

t 1 Δt t1Δt

Ω

F^ i δui dt 1 Δt Ω 1

ð

t 1 Δt t1Δt

ST

T^ i δui dt 1 Δt S;

(11.2)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

411

t+Δt t t+Δt

ST

t=0

x3 0

t

ST

t+Δt

t

ST

0

Ω

Ω

Ω δui

o x2

Su

p ( 0 xi )

x1

p( t+Δt xi )

FIGURE 11.1 Motion and its configuration at different time instants: Su is the displacement boundary on which u3 5 0 is imposed, and t xi 5 0 xi 1 t ui , t 1 Δt xi 5 0 xi 1 t 1 Δt ui 5 t xi 1 ui , the increment displacement ui 5 t 1 Δt xi 2 t xi .

F^ i is the external body force vector per unit volume at time t 1 Δt; t 1 Δt T^ i is the external traction force vector per unit surface area at time t 1 Δt; t 1 Δt ST is the service of the body at time t 1 Δt on which external traction forces are applied. Compared with linear analysis in which it is assumed that the displacements are infinitesimally small so that in Eqs. (11.1)(11.3) the original configuration is used, now for nonlinear analysis, a fundamental difficulty in the general equations given by Eqs. (11.1)(11.3) is that the configuration of the body at time t 1 Δt is unknown. With large deformations and motions, the configuration of the body undergoes a continuous change, which results in some important consequences for the development of nonlinear FE analysis. For example, the Cauchy stress at time t 1 Δt cannot be obtained simply by adding to the one at time t a stress increment caused only by the strain of the material. It must also take into account the rigid body rotation of the material because the Cauchy stress tensor also changes due to that rotation. In Eq. (11.1), the first term on the left-hand side is the internal virtual work, and the second term is the virtual variation of the kinetic energy of the body; the right-hand side is the external virtual work. Using the Green theorem and considering that the virtual displacement δui is arbitrary, from Eqs. (11.1) and (11.2), we obtain the following dynamic equilibrium equation describing the dynamics of the body at time t 1 Δt: t 1 Δt

t 1 Δt

σij 1 t 1 Δt F^ i 5 t 1 Δt ρS u€i ;

in

where t 1 Δt ηj is the unit normal vector to the surface

11.1.2

t 1 Δt

t 1 Δt

Ω; t 1 Δt σij t 1 Δt ηj 5 t 1 Δt T^ i ;

on

t 1 Δt

ST ;

(11.3)

ST at time t 1 Δt, pointing from inside to outside the body.

Expression of stress virtual work

The principle of virtual work in Eq. (11.1) is given in the form based on the spatial coordinate system, which must transform into its corresponding form in the Lagrange coordinate system for FE analysis in solid mechanics. Based on the fundamental knowledge presented in Chapter 3, Fundamentals of continuum mechanics, we can obtain the following formulations. From Eq. (3.26), it follows that the variation of the Green strain tensor is given by "    # 1 @xm @xm @xm @xm δEij 5 1 δ δ 2 @Xi @Xj @Xi @Xj " # 1 @δxm @xm @xm @δxm 5 1 2 @Xi @Xj @Xi @Xj " # 1 @δum @xm @xm @δum (11.4) 5 1 2 @Xi @Xj @Xi @Xj (   @xr 1 @δum @δur @xm @xr @xm 5 1 5 δerm ; @Xi 2 @xr @xm @Xj @Xi @Xj δerm 5

@Xi @Xj δEij ; @xr @xm

412

FluidSolid Interaction Dynamics

where the definition of the variation δerm given in Eq. (11.1) and the relationship δxm 5 δum ;

@δum @δum @xr 5 ; @Xi @xr @Xi

(11.5)

as well as the rules of tensor index presented in Chapter 2, Cartesian tensor and matrix calculus, have been used. Furthermore, the Cauchy stress tensor given by Eq. (3.50) is expressed by the second PiolaKirchhoff stress tensor, that is, σij 5 J 21

@xi @xj ρ @xi @xj Σrs 5 Σrs ; ρ0 @Xr @Xr @Xs @Xs

(11.6)

where ρ0 5 ρJ in Eq. (3.72) is introduced. Substituting Eqs. (11.4) and (11.6) into the first term on the left-hand side of Eq. (11.1) and considering the configuration of the body at time t, we obtain   ð ð t ρ @xi t @xj @Xn t @Xl t t σij δt eij dt Ω 5 Σ δ E dΩ 0 rs t t @Xs @xi 0 nl @xj Ω Ω ρ0 @Xr (11.7)  ð t ð t ð ρ t ρt t t 0 t t t t 5 0 Σnl δ 0 Enl d Ω; ρ δnr 0 Σrs δ ls δ0 Enl d Ω 5 ρ 0 Σnl δ 0 Enl d Ω 5 t

Ω

0

t

Ω

0

0

Ω

where ρd Ω 5 ρd Ω is used. Eq. (11.7) indicates that the Cauchy stress and the Almansi strain are a pair of variables to calculate the work for the current reference configuration, while the second PiolaKirchhoff stress and the Green strain are another pair of variables for the work calculation based on the Lagrangian configuration. The work is scalar so that it is not changed during the coordinate transformation process. Eq. (11.7) is the fundamental expression of the total Lagrangian (TL) formulation and of the UL formulation used in the FE analysis of solid/structure mechanics. In Eq. (11.7), the final integration is performed over the initial volume of the body; any other previously calculated configuration with the coordinate τ xi at time τ; ðτ , tÞ; could be a reference configuration to define the second PiolaKirchhoff stress and the Green strain as tτ Σ nl and tτ E nl , respectively. Eq. (11.7) can be written in a more general form: ð ð τ t t t σ ij δt eij d t Ω 5 (11.8) τ Σ nl δτ E nl d Ω: t

t

0

0

t

Ω

τ

Ω

Based on Eq. (11.8), we can recall that the solutions from times 0; Δt; 2Δt; . . .; t have already been calculated, so we can employ Eq. (11.8) and refer the stress and strain to one of these known equilibrium configurations. In principle, the options lies essentially between two formulations that have been termed TL and UL formulations (see Bathe, 1996). In TL formulation, all static and kinematic variables are referred to the initial configuration at time 0. However, in UL formulation all static and kinematic variables are referred to the last calculated configuration. The only advantage of choosing one formulation rather than the other mainly depends on its greater numerical efficiency. For nonlinear FSI analysis, UL formulation may be more effective to match the CFD model based on ALE formulation. For comparison purposes, we list both TL and UL formulation in the following subsections.

11.1.3

Total Lagrangian formulation in solid dynamics

11.1.3.1 Dynamic equation From Eqs. (11.1), (11.2), and (11.8), in the TL formulation, we develop the numerical model based on the basic equation ð  ð ð  0 0 t 1 Δt ^ t 1 Δt ^ t 1 Δt t 1 Δt 0 t 1 Δt t 1 Δt 0 € F Σ δ E 2 ρ δu Ω 5 W ; W 5 δu d Ω 1 u d (11.9) 0 T i δui d S; 0 i i S i i 0 ij 0 ij 0 0 0 Ω

where

Ω

ST

3 2 t 1 Δt ui;j 1 t 1 Δt0 uj;i 1 t 1 Δt0 uk;i t 1 Δt0 uk;j 0 @X @X i t 1 Δt j t 1 Δt 5: σmn ; δt 1 Δt0 Eij 5 δ4 0 Σij 5 t 1 Δt 2 ρ @xm @xn 0

ρ

(11.10)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

413

11.1.3.2 Incremental decompositions To arrive at the linearized equation of motion, incremental decompositions are needed: t 1 Δt 0 Σij

0 eij

5

5 t0 Σij 1 0 Σij ;

t 1 Δt 0 ui

t 1 Δt 0 Eij

5 t0 ui 1 0 ui ;

5 t0 Eij 1 0 Eij ;

0 Eij

5 0 eij 1 0 εij ;

0 εij

5

1 0 uk;i 0 uk;j ; 2

 1 t t 0 ui;j 1 0 uj;i 1 0 uk;i 0 uk;j 1 0 uk;i 0 uk;j : 2 (11.11) t 0 ðÞ

at time t and Here, the stress and strain at time t 1 Δt are, respectively, represented by a summation of the one its increment 0() from time t to time t 1 Δt. Furthermore, the incremental strain 0 Eij is represented by a summation of the linear part 0 eij and the nonlinear part 0 εij . In the linear part 0 eij , the last two terms on the right-hand side denote the effect of the initial displacement t0 ðÞ.

11.1.3.3 Dynamic equation with incremental decompositions Calculating δt 1 Δt0 Σ ij 5 δ0 Σij due to the stress t0 E ij at time t and substituting Eq. (11.11) into Eq. (11.9), we obtain the dynamic equation ð  ð  0 0 t t 0 t 1 Δt € Σ δ E 1 Σ δ ε 2 ρ δu Ω 5 W 2 u d (11.12) 0 ij 0 ij 0 Σij δ0 eij d Ω: S i i 0 ij 0 ij 0 0 Ω

Ω

11.1.3.4 Linearization of dynamic equation Using the approximations 0 Σij

5 0 Cijkl 0 ekl ;

δ0 Eij 5 δ0 eij ;

we obtain an approximate dynamic equation in the form ð   t 0 €i δui d0 Ω 5 0 Cijkl 0 ekl δ0 eij 1 0 Σij δ0 εij 2 ρS u 0

Ω

t 1 Δt

(11.13) ð

W2

0

Ω

0 t 0 Σij δ0 eij d Ω:

(11.14)

Here, 0 Cijkl denotes the incremental stressstrain tensor at time t referring to the configuration at time 0.

11.1.3.5 Matrix equation of displacement-based element Dynamic equation with implicit time integration   € 1 t K 1 t K U 5 t 1 Δt R 2 t F; Mt 1 Δt U 0 L 0 NL 0

(11.15a)

and its explicit time integration form € 5 t R 2 t F; Mt U 0 can be derived by considering the dynamic equilibrium of the system subject to the force increment, that is,   € 1 t K 1 t K U 5 t 1 Δt R 2 t R; MU 0 L 0 NL

(11.15b)

(11.15c)

where M is the time-independent mass matrix; t0 KL is the linear strain incremental stiffness matrix; t0 K NL is the nonlinear strain incremental stiffness matrix; t 1 Δt R is the vector of externally applied nodal point loads at time t 1 Δt; t0 F is the vector of nodal force equivalent to the element stress at time t; U is the vector of increments in the nodal point dis€ is the vector of nodal point acceleration at time t 1 Δt. placements; t 1 Δt U Finite element matrices Integral and evaluation for mass matrix:

ð 0 0

Ω

ρS u€i δui d 0 Ω;

(11.16a)

414

FluidSolid Interaction Dynamics

Mt1Δt u€ 5



t1Δt 0

Ω

Integral and evaluation for external load matrix: ð ð 0 t 1 Δt ^ t 1 Δt W5 δu d Ω 1 F i 0 i 0

t 1 Δt

ð R5

Ω

0

HT 0

Ω

€ u:

(11.16b)

t 1 Δt ^ 0 0 T i δui d S;

(11.17a)

t 1 Δt ^ 0 0 Td S

(11.17b)

HT 0 ρS Hd 0 Ω

t 1 Δt ^ 0 0 Fd Ω

ð 1

ST

HT 0

ST

Integral and evaluation for stiffness matrices and nodal point stress vector: ! ð ð 0 0 t t t T 0 Cijkl 0 ekl δ 0 eij d Ω; 0 KL u 5 0 B L 0 C 0 B L d Ω u; 0

Ω

0

Ω

!

ð

ð 0

Ω

0 t 0 Σ ij δ 0 εij d Ω;

t 0 K NL u 5

0

ð

Ω

t T t t 0 B NL 0 Σ 0 B NL

ð

0

Ω

0 t 0 Σ ij δ 0 eij d Ω:;

t 0F

5

0

Ω

(11.18a)

d 0 Ω u;

(11.18b)

0 t T t ^ 0 B L 0 Σd Ω:

(11.18c)

In these equations, the following notations are used for calculation of the element matrices: u is the element displacement vector, which is a time function for dynamics; H is the element displacement interpolation matrix, from which the element displacement vector is interpolated by uðx; tÞ 5 HðxÞuðtÞ; t 1 Δt0 F^ is the element body force vector at ^ is the element traction force vector at time t 1 Δt; t B ; t B time t 1 Δt; t 1 Δt0 T 0 L 0 NL is the linear, nonlinear ^ is the matrix and straindisplacement transformation matrices, from which it follows 0 e 5 t0 BL u; 0 ε 5 t0 B NL u; t0 Σ; t0 Σ vector of second PiolaKirchhoff stresses; and 0 C is the increment stressstrain material property matrix, 0 S 5 0 C 0 e:

11.1.4

Updated Lagrangian formulation in solid dynamics

11.1.4.1 Dynamic equation From Eqs. (11.1), (11.2), and (11.8), in the UL formulation we develop the numerical model based on the basic equation referring to the configuration at time t: ð ð ð ð t 0 0 t 1 Δt ^ t 1 Δt ^ t 1 Δt t 1 Δt 0 0 t 1 Δt t 1 Δt € F Σ δ E d Ω 2 ρ δu d Ω 5 W ; W 5 δu d Ω 1 u i i i 0 T i δui d S; (11.19a) t ij 0 i S t ij t

Ω

0

Ω

0

where

2 t 1 Δt t Σij

ρ @Xi t 1 Δt @Xj t 1 Δt 4 5 t 1 Δt σmn ; δ t Eij 5 δ ρ @xm @xn t

t 1 Δt t u i;j

Ω

0

ST

1 t 1 Δtt u j;i 1 t 1 Δtt u k;i t 1 Δtt u k;j 2

3 5:

(11.19b)

11.1.4.2 Incremental decompositions For the linearized equation of motion, the corresponding incremental decompositions are as follows: t 1 Δt t Σij t 1 Δt t ui

5 tt Σij 1 t Σij 5 t σij 1 t Σij ;

5 t ui ;

t 1 Δt t Eij

5 t Eij ;

t Eij

t t Σij

5 t σij ;

5 t eij 1 t εij

t eij

5

 1 t ui;j 1 t uj;i ; 2

t εij

5

1 t uk;i t uk;j : 2

(11.20)

Here, the second PiolaKirchhoff stress at time t is just the Cauchy’s stress, and the Green strain increment is just the summation of linear and nonlinear parts in the straindisplacement Eq. (3.30).

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

415

11.1.4.3 Dynamic equation with incremental decompositions Substituting Eq. (11.20) into Eqs. (11.19a) and (11.19b), we obtain the dynamic equation ð ð ð t 0 t t 1 Δt 0 ðt Σij δt Eij 1 σij δt εij Þd Ω 2 ρS u€i δui d Ω 5 W 2 t σij δt eij dt Ω: t

Ω

0

Ω

t

Ω

(11.21)

11.1.4.4 Linearization of dynamic equation Using the approximations t Σij

5 t Cijkl t ekl ;

we obtain an approximate dynamic equation in the form ð  ð  t t t Cijkl t ekl δ t eij 1 σ ij δ t εij d Ω 2 t

Ω

0 0

Ω

δt Eij 5 δ t eij ;

ρS u€i δui d Ω 5 0

(11.22)

t 1 Δt

ð W2

t t

Ω

σij δt eij d t Ω:

(11.23)

Here, t Cijkl denotes the incremental stressstrain tensor at time t.

11.1.4.5 Matrix equation of displacement-based element Dynamic equation with implicit time integration   € 1 t K 1 t K U 5 t 1 Δt R 2 t F Mt 1 Δt U t L t NL t

(11.24a)

and its explicit time integration form € 5 tR 2 tF Mt U t

(11.24b)

can be obtained by noting the increment equilibrium condition when the system is subjected to a force increment, that is,   € 1 t K 1 t K U 5 t 1 Δt R 2 t R; MU (11.24c) t L t NL where M is the time-independent mass matrix; tt KL is the linear strain incremental stiffness matrix; tt K NL is the nonlinear strain incremental stiffness matrix; t 1 Δt R is the vector of externally applied nodal point loads at time t 1 Δt; tt F is the vector of nodal force equivalent to the element stress at time t; U is the vector of increments in the nodal point dis€ is the vector of nodal point acceleration at time t 1 Δt. placements; t 1 Δt U Finite element matrices Integral and evaluation for mass matrix:

ð 0 0

Mt1Δt u€ 5

Ω

ρS u€i δui d 0 Ω;



 0

Ω

HT 0 ρS Hd 0 Ω

Integral and evaluation for external load matrix: ð ð t 1 Δt ^ 0 t 1 Δt W5 0 F i δui d Ω 1 0

t 1 Δt

ð R5

0

Ω

Ω

(11.25a)

0

^ 0Ω HT t 1 Δt0 Fd

ð 1

0

ST

ST

t1Δt

€ u:

t 1 Δt ^ 0 0 T i δui d S;

(11.26a)

^ 0 S: HT t 1 Δt0 Td

(11.26b)

Integral and evaluation for stiffness matrices and nodal point stress vector: ! ð ð t t t t t T t C ijkl t e kl δ t e ij d Ω; t KL u 5 t B L t C t BL d Ω u; t

Ω

t

Ω

(11.25b)

(11.27a)

416

FluidSolid Interaction Dynamics

ð

!

ð t t

Ω

t t K NL u 5

σij δt ε ij d Ω; t

t

ð

u;

(11.27b)

ð t

t

Ω

t t T t t t B NL σ t B NL d Ω

Ω

σij δt e ij dt Ω;

t tF

5

t

Ω

t T t ^ t Ω: t B L σd

(11.27c)

In these equations, the following notation is used for calculation of the element matrices: u is the element displacement vector, which is a time function for dynamics; H is the element displacement interpolation matrix, from which the element displacement vector is interpolated by uðx; tÞ 5 HðxÞuðtÞ; t 1 Δt0 F^ is the element body force vector at time t 1 Δt; t 1 Δt ^ t t 0 T is the element traction force vector at time t 1 Δt; t BL ; t B NL is the linear, nonlinear straindisplacement t t transformation matrices, from which it follows t e 5 t BL u; t ε 5 t B NL u; t σ; t σ^ are the matrix and vector of Cauchy stresses; t C is the increment stressstrain material property matrix, t S 5 t C t e:

11.1.5

Solution of nonlinear equations

The dynamic FE equations derived in TL and UL formulation are nonlinear, in which the matrices involve the unknown variables. To solve these nonlinear equations, an iteration scheme is required. The most often adopted scheme is the NewtonRaphson iteration as follows.

11.1.5.1 NewtonRaphson iteration A nonlinear equation can be represented in a generalized vector form: ^ 5 R; fðUÞ

(11.28a)

where   ^ T 5 U^ 1 f T 5 f1 f2 ? fN ; U  RT 5 R1 R2 ? RN :

U^ 2

? U^ N ;

Assume that in the iterative solution, we have evaluated UðI21Þ ; then a Taylor series expansion gives

  ^ 2 UðI21Þ 1 higher-order terms: ^ 5 fðUðI21Þ Þ 1 @f

U fðUÞ @U ðI21Þ

(11.28b)

(11.28c)

Substituting Eq. (11.28c) into Eq. (11.28a) and neglecting the higher-order terms, we can calculate an increment in the variable vector:

  @f

; KðI21Þ ΔUðIÞ 5 R 2 f UðI21Þ ; KðI21Þ 5 @U

ðI21Þ (11.28d)     ^ 2 UðI21Þ 5 R 2 f UðI21Þ ; ΔUðIÞ 5 U where KðI21Þ is the current tangent stiffness matrix, of which an element KijðI21Þ 5

@fi ; @Ui

(11.28e)

and the improved solution is UðIÞ 5 UðI21Þ 1 ΔUðIÞ : ð0Þ

(11.28f)

ð0Þ

The initial conditions in this iteration are K and U previously determined according to the physical problems. The iteration is continued until convergence criteria are satisfied, such as the relative error less than a convergence tolerance εD : :ΔUðIÞ : # εD : :UðIÞ :

(11.28g)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

417

11.1.5.2 Solution of explicit nonlinear dynamics equation The explicit integration Eq. (11.15b) in TL formulation and the one in Eq. (11.24b) for UL formulation can be rewritten as € 5 t R 2 t F; Mt U

(11.29a)

The solution for the node displacement at time t 1 Δt is obtained using the central difference approximation for € in preceding equation: accelerations given by Eq. (A.20.1) to substitute for t U t

€ 5 U

t 2 Δt

U 2 2 t U 1 t 1 Δt U ; Δt2

(11.29b)

from which, when substituted into Eq. (11.29a), it follows t 2 Δt

M

U 2 2 t U 1 t 1 Δt U t 5 R 2 t F; Δt2

(11.29c)

that gives the solution t 1 Δt U. As discussed in Appendix A, this explicit integration solution simply corresponds to a forward progression in time, and especially when a diagonal mass matrix is used, the solution of t 1 Δt U does not involve a triangular factorization of a coefficient matrix. The central difference method is conditionally stable, which requires that the time increment step satisfies Eq. (A.22), with the smallest time period Tm in the FE system involved. For nonlinear analysis, this time period is not constant during response calculation, the time step Δt needs to be decreased if the system stiffens, and its adjustment must be performed conservatively so that the stable condition in Eq. (A.22) is certain to be satisfied at all times.

11.1.5.3 Solution of implicit nonlinear dynamics equation The implicit integration Eq. (11.15a) in TL formulation and the one in Eq. (11.24a) for UL formulation can be rewritten as € 1 t KΔU 5 t 1 Δt R 2 t F; M t 1 Δt U

(11.30a)

from which, the corresponding NewtonRaphson iteration can be written as € ðIÞ 1 t KΔUðIÞ 5 t 1 Δt R 2 M t 1 Δt U

t1Δt ðI21Þ

F

;

t 1 Δt

UðIÞ 5 t 1 Δt UðI 2 1Þ 1 ΔUðIÞ :

(11.30b)

The implicit time integration schemes discussed in Appendix A can be used in nonlinear dynamic response calculations. For example, based on Newmark’s method, we have the following assumptions: t 1 Δt

Δt  t 1 Δt _ t _  U1 U ; 2

U 5 tU 1

t 1 Δt

  _ 5 tU _ 1 Δt t 1 Δt U € 1 tU € ; U 2

(11.30c)

which, combined with Eq. (11.30b), gives t 1 Δt

  € ðIÞ 5 4 t 1 Δt UðI 2 1Þ 2 t U 1 ΔUðIÞ 2 4 t U _ 2 t U: € U Δt2 Δt

Substituting Eq. (11.30d) into Eq. (11.30b), we obtain 4  t 1 Δt ðI 2 1Þ t  4 t_ t€ t ^ ^ 5 t K 1 4 M: KΔUðIÞ 5 t 1 Δt R 2 t1Δt FðI21Þ 2 M U 2 U ; tK U 2 U 2 Δt2 Δt Δt2

(11.30d)

(11.30e)

The convenience is achieved when the force tolerance RTOL and the energy tolerance ETOL, which can be defined by computational accuracy requirement, such as 1026, are satisfied, that is, : t 1 Δt R 2

t1Δt ðI21Þ

F

€ ðI 2 1Þ : 2 M t 1 Δt U

€ : t 1 Δt R 2 t F 2 M t U: ΔUðIÞT



t 1 Δt

 € ðI 2 1Þ 2 M t 1 Δt U  # ETOL: t 1 Δt € R 2 tF 2 M tU

R2 

ΔUð1ÞT

# RTOL;

(11.30f)

t1Δt ðI21Þ

F

(11.30g)

418

FluidSolid Interaction Dynamics

tu 1

o

= ûeiωt

h tu

2

x

FIGURE 11.2 One-dimensional elastic rod inserted into the water and subject to a harmonic excitation at its top end.

11.1.6

A one-dimensional example

Here, we discuss a simple 1D example to illustrate the formulations presented in Sections 11.4 and 11.5. As shown in Fig. 11.2, a uniform elastic rod of length L, cross section area A 5 1, mass density ρs , and elastic modulus C is excited ^ iωt at its top end while the lower part of length h , L is inserted into the water. We by a harmonic motion u1 5 ue assume that the rod undergoes a large motion with nonlinear geometrical deformation but neglects its change of cross section area during the motion. Also, we do not consider the radiation effect on the water pressure to find the dynamic water pressure applied at the lower end of the rod. We solve this simple dynamic problem referring to its static equilibrium position in Fig. 11.2 as follows. For simplicity, we consider the rod as a two-node element; its nodal point coordinates are used to determine the configuration of the rod at a time based on the following interpolations.

11.1.6.1 Coordinate interpolations 0

  x r 5 H0 x;

t

  x r 5 H t x;

h H 5 h1

i h2 ;

x 5

0 T

h

0

i x1 0 x2 ;

x 5

t T

h

t

i x1 t x2 ;

(11.31a)

where ð1 2 rÞ ; 2

h1 5 r5

0

L 5 0 x2 2 0 x1 5 L;

t

L 5 t x2 2 t x1 ;

  20 x 2 0 L 0

L

h2 5 5

ð1 1 rÞ ; 2

  2 tx 2 tL t

dr 2 dr 2 50 ; t 5 t ; 0 d x L L dx

L

t

(11.31b) ;

J5

t dt x L ð0 L 1 ΔLÞ 5 5 : 0 0 d0 x L L

(11.31c)

11.1.6.2 Displacement interpolations Based on an isoperimetric displacement element formulation, the displacement at a given time and its increment are, respectively, given by the following formulation: t

uðrÞ 5 Ht u;

uðrÞ 5 Hu; t uT 5

t

u1 t u2 ;

uT 5 ½ u1

u2 ;

du dr 5 H;r u; d0 x d0 x

dt u dr 5 0 H;r t u: d0 x d x

(11.32a)

Here, we can also denote the displacement at time t in the form t

u 5 t x 2 0 x;

dt u dt x 5 0 2 1 5 t J 2 1; 0 d x d x

which gives the same result as obtained by Eq. (11.32a).

(11.32b)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

419

11.1.6.3 Total Lagrangian formulations Longitudinal strain t 0 E 11

5

dt u 1 dt u dt u 1 : d0 x 2 d0 x d0 x

(11.33)

Increment strain 0 E 11

5 0 e 11 1 0 ε 11 ;

0 e 11

5

du dt u du 1 0 0 ; 0 d x d xd x

0 ε 11

1 du du : 2 d0 x d0 x

5

(11.34)

Straindisplacement matrices From Eqs. (11.32a), (11.32b), and Eq. (11.34), it follows 0 1 t t t t @1 1 d u A du 5 t J du 5 2 J H;r u 5 0J ½ 21 1 u 5 t B u; t B 5 0J ½ 21 1 ; e 5 0 11 0 L 0 L 0 0 0 0 L L d x d x d x L 0 12 0 12 1 du du 1@ dr A T T dr 5 u H;r H;r u; δ0 ε11 5 @ 0 A δuT HT;r H;r u 5 δuT t0 B TNL t0 BNL u; 0 ε11 5 2 d0 x d0 x 2 d0 x d x t 0 BNL

5

dr 1 H;r 5 0 ½ 21 0 d x L

(11.35)

1 :

Stiffness matrices and nodal point stress vector From Eqs. (11.18a)(11.18c), it follows ðL t 2 21  t T t t 0 5 J C 21 K 5 B C B d x 0 L 0 L 0 L 0 1 L 0 t 0 KNL

5

ðL 0

t 0 t T t 0 B NL 0 Σ11 0 BNL d x;

5

0

2

L t σ 11 0

2

L tL



21  1

2 tL C 1 1 5 3 21 L 21

21 ; 1

σ 11 1 1 5 t L 21

t

21 1

(11.36a) ;

(11.36b)

where t σ 11 is the current Cauchy stress equaling the current force divided by the current cross section area (assuming its change is neglected here). Therefore, from Eq. (11.10), we have the second PiolaKirchhoff stress  2 0 2 0 0 ρ 0L t L t L t t Σ 5 σ 5 J σ 11 5 t t σ 11 : (11.36c) 11 0 11 t t tρ L L L ðL t 0 J L 21 0 L t 21 t t T t ^ t 0 F 5 B d x 5 σ 5 σ 11 : (11.36d) Σ 0 L 0 0 t 0 1 1 L 11 L 0 11 Mass matrix and external load vector From Eq. (11.16b), we obtain the mass matrix ð ð0x 0 ð1 0 0 2 L L ρS 1 1 2 r  12r HT 0 ρS Hd0 x 5 HT 0 ρS Hdr 5 M5 2 21 8 0x 21 1 1 r 1

0 0 L ρS 1 1 1 r dr 5 1=2 3

1=2 : (11.37) 1

Noting that we choose the static equilibrium of the system as our reference state, the body force is neglected. The dynamic water force t 1 Δt ^ 0T

5 2 t 1 Δt0 pðtÞ 5 2 ρf g t 1 Δt0 u2

(11.38a)

420

FluidSolid Interaction Dynamics

is applied at the rod bottom and is a function of the displacement at the nodal point 2. Also, to keep the prescribed motion at the top node 1, a reaction force t 1 Δt0 T 1 is generated at node 1 during the motion. Therefore we obtain the external load vector # " t 1 Δt # " " # t 1 Δt 0 0 t  0T 1 0T 1 t 1 Δt 2 ρf g (11.38b) 5 R5 0u 1 u : t 1 Δt 0 1 2 0 0p Dynamic Eq. (11.15) L ρS 1 1=2 3

0 0

1=2 1



t 1 Δt



u€1 u€2



0

2

t

LC 1@ 3 L



t σ 1 21 1 t 11 1 L 21

1 21

1 t 1 Δt 21 A u1 21 t T1 σ11 : 2 5 1 u2 1 2 t 1 Δt pðtÞ (11.39a)

Here, t 1 Δt pðtÞ denotes the dynamic water pressure applied at node 2 at time t 1 Δt, which generally is determined by considering the governing equation of the fluid dynamics for FSI problems. However, for this simple example, we are neglecting the effect of radiation and free surface wave on the water pressure and are involving only the static water pressure change by the rod motion, so that the water pressure is given by Eqs. (11.38a), (11.38b), and (11.39a) can be further written as ! t 0 t 2 0 0 u1 21 1=2 t 1 Δt u€1 21 L 0 ρS 1 LC 1 σ 11 1 1 t 1 1 ρf g L3 3 L 21 0 1 21 1 u2 1=2 1 1 u€2 (11.39b) t 1 Δt 0 0 t u1 21 t T1 σ11 1 5 2 ρf g : 2 0 1 u2 1 0  In the preceding dynamic equations, the uT 5 u1 u2 represents the increment vector of the nodal point displacements; the nodal displacement, velocity, and acceleration at time t 1 Δt are calculated by t 1 Δt

t 1 Δt

u 5 t u 1 u;

_ u_ 5 t u_ 1 u;

t 1 Δt

€ u€ 5 t u€ 1 u;

(11.40a)

from which the length and the rod’s section stress at time t 1 Δt can be obtained by t 1 Δt

L 5 t 1 Δt x2 2 t 1 Δt x1 ;

Δ t 1 Δt L 5 t 1 Δt L 2 L;

t 1 Δt

σ 11 5 C11 Δ t 1 Δt L:

(11.40b)

Introducing the displacement boundary condition ^ iωt ; therefore the node displacement, acceleration, and its displacement At node 1, the motion is given by t u1 ðtÞ 5 ue increment at node 1, respectively, take the following values: t

^ iωt ; u1 5 ue

t

^ iωt ; u€1 5 2 ω2 ue

0

u1  0;

(11.41a)

from which, when substituted into Eqs. (11.39a) and (11.39b), it follows that 0 0 L ρS t 1 Δt L ρS ω2 u^ iωðt1ΔtÞ e ; u€2 1 t K u2 5 2 t 1 Δt pðtÞ 2 t σ 11 1 6 3

(11.41b)

0 0 L 0 ρS t 1 Δt L ρS ω2 u^ iωðt1ΔtÞ e ; u€2 1 ðt K 1 ρf gÞu2 5 2 ρf g t u2 2 t σ11 1 6 3

(11.41c)

0 0

0

where the tangent stiffness t

K5

t

2

t LC σ 1 t 11 : 3 L L

(11.41d)

These two equations derived from the second line of the matrix Eqs. (11.39a) and (11.39b) govern the motion of node 2 from which its dynamic variables can be solved. The first line of matrix Eqs. (11.39a) and (11.39b) provides another equation to determine the dynamic node force at node 1 to keep its motion as prescribed, that is,

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

0

0 0 L 0 ρS t 1 Δt L ρS ω2 u^ iωðt1ΔtÞ e 5 2 t 1 Δt T 1 : u€2 2 t K u2 1 t σ 11 2 3 6

421

(11.41e)

11.1.6.4 Updated Lagrangian formulations Longitudinal strain t 1 Δt t E 11

5 t E11 :

(11.42)

Increment strain t E 11 5 t e 11 1 t ε 11 ;

t e 11 5

dt u ; dt x

t ε 11 5

1 dt u dt u : 2 dt x dt x

(11.43)

Straindisplacement matrices From Eqs. (11.31a)(11.31c) and (11.43), it follows d t u dr 2 1 5 t H;r u 5 t 21 1 u 5 tt BL u; t dr d x L L  t  t  1 d u dr 2 T T d u dr ε 5 t 11 dr d t x 5 t 2 u H;r H;r u; 2 dr dt x L t e11

5

δt ε11 5

4 t

L

2

δuT HT;r H;r u 5 δuT tt B TNL tt BNL u;

t t BNL

t t BL

5

1 21 L

t

1 ;

(11.44)

2 1 5 t H;r 5 t 21 L L

1 :

Stiffness matrices and nodal point stress vector From Eqs. (11.27a)(11.27c), it follows ð1 t Lt T t C 1 C 21  21 t K 5 21 1 5 t ; t L t B L C t BL dr 5 t 1 L 1 L 21 21 2 ð1 t t t σ 11 1 Lt T t σ 11 21  21 t t 21 1 5 ; K 5 B Σ B dr; 5 t NL t NL t 11 t NL t t 1 1 L L 21 21 2

(11.45a)

(11.45b)

where t σ 11 5 tt Σ 11 is the current Cauchy stress equaling the current force divided by the current cross section area (assuming its change is neglected here). The nodal point stress vector ð1 t L t T t^ 21 t t F 5 B 11dr 5 σ 11 : (11.45c) Σ t t L t 1 21 2 Mass matrix and external load vector These are same as the those given in Eqs. (11.37)(11.38a) and (11.38b) (with the subindex 0 deleted), with the integration over the original reference configuration: 0 0 L ρS 1 1=2 M5 ; (11.46) 1=2 1 3 t1Δt t1Δt 0 0 t T1 T1 t 1 Δt R5 ð u 1 uÞ: (11.47) 5 2 ρf g t1Δt 0 1 2 p 0 Dynamic equation Eqs. (11.24a)(11.24c) now takes the following forms:  t 0 0 1=2 t 1 Δt u€1 L ρS 1 C σ 1 21 1 1 t 11 1 t 3 1=2 1 u€2 1 L 21 L 21

21 1



t 1 Δt u1 21 t T1 σ11 ; 2 5 u2 1 2 t 1 Δt pðtÞ (11.48a)

422

FluidSolid Interaction Dynamics

0

L 0 ρS 1 3 1=2

1=2



t 1 Δt



1

u€1 u€2



t  C 1 σ 11 1 u1 21 21 0 0 1 1 ρ g f t t 1 1 0 1 u2 L 21 L 21 t 1 Δt 0 0 t u1 21 t T1 5 2 ρf g σ11 1 : 2 0 1 u2 1 0 

1

(11.48b)

Introducing the displacement boundary condition At node 1, Eqs. (11.48a) and (11.48b) become 0 0 2 L 0 ρS t 1 Δt u€2 1 t K u2 5 2 t 1 Δt pðtÞ 2 t σ 1 L ρS ω u^eiωðt1ΔtÞ ; 11 3 6

(11.48c)

0 0 L 0 ρS t 1 Δt L ρS ω2 u^ iωðt1ΔtÞ e ; u€2 1 ðt K 1 ρf gÞu2 5 2 ρf gt u2 2 t σ11 1 6 3

(11.48d)

0

0

where the tangent stiffness is now in the form t

K5

t

σ 11 1 C : t L

(11.48e)

The reaction force at node 1 can also be calculated by Eq. (11.41e).

11.1.6.5 Iteration equation Considering the UL formulation in Eq. (11.48c) and using the implicit iteration in Eq. (11.30e), for this example, we obtain the iteration equation t

0 0 2 ^ iωðt1ΔtÞ t 1 Δt ðI 2 1Þ 2 t 1 Δt σ ðI 2 1Þ 1 L ρS ω u e K^ ΔuðIÞ 5 2 p 11 2 6 # " 0 0 L ρS 4 t 1 Δt ðI 2 1Þ t  4 t 40 L 0 ρS t σ 11 1 C 40 L 0 ρS t ^ t t _ € ; K 5 1 2 u 2 u 2 K 1 5 t : u u 2 2 2 2 Δt2 Δt 2 3 3Δt2 L 3Δt2

(11.49a)

With its convergence conditions, we simply use :ΔuðIÞ 2 :

: :uðI21Þ 2

# εD :

(11.49b)

When the convergence is reached, we can calculate the values of variables at time t 1 Δt using the following equations:  4 4 t 1 Δt ðIÞ t 1 Δt t 1 Δt ðI 2 1Þ t 1 Δt t 1 Δt ðIÞ t t € € u2 5 t 1 Δt u ðIÞ 5 u 1 Δu ; 5 5 u 2 u u u u_ 2 t u€2 ; 2 2 2 2 2 2 2 2 Δt2 Δt 2  Δtt 1 Δt t 1 Δt ^ iωðt1ΔtÞ ; u_2 5 t u_2 1 u€2 1 t u€2 ; t 1 Δt L 5 t 1 Δt u2 2 t 1 Δt u1 5 t 1 Δt u2 2 ue (11.49c) 2 t 1 Δt

σ11 5

Cðt 1 Δt L 2 LÞ t 1 Δt p 5 ρf gt 1 Δt u2 : ; t 1 Δt L

Based on these results, we can enter into the next time period from t 1 Δt to t 1 2Δt. For the FSI case with the radiation and free surface wave considered, the water pressure applied on the rod cannot be given by the simple formulation in Eq. (11.49c). In this case, we must consider fluid dynamic analysis, which will be completed after a discussion on CFD formulation in the next subsection of this chapter.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

11.2

423

Updated arbitrary LagrangianEulerian formulations in computational fluid dynamics

As we have mentioned, for nonlinear FSI analysis, an ALE reference system is needed to keep the validity of compatibility conditions on the FSI interfaces. In Sections 3.1.3 and 3.1.5, we have briefly presented the fundamental concepts and formulations of ALE and updated ALE systems, which were also used to establish the governing equations in the integration and differential forms for continuum mechanics in Sections 3.5.9.3 and 3.5.9.4. Here, based on this knowledge, we will construct ALE/updated ALE formulations in CFD, so that the numerical model in CFD can be combined with the UL model in solid mechanics to deal with nonlinear FSI dynamic problems. For readers to more deeply understand ALE system and its applications, we give a short review on its creation and development as follows. More detailed review papers with listed historical references may be found in the review papers, such as Donea and Belytschko (1992), Huerta et al. (1992), Huerta (1993), and Zhang et al. (1997). Furthermore, an in-depth survey of ALE methods, including both conceptual aspects and numerical implementation details in view of the applications in large deformation material response, fluid dynamics, nonlinear solid mechanics, and coupled FSI problems, by Donea et al. (2004) (https://doi.org/ 10.1002/0470091355.ecm009) can be read on the webpage http://onlineliberary.wiley.com. Interested readers may refer to these important publications.

11.2.1

History and development of arbitrary LagrangianEulerian description

As discussed in Chapter 3, Fundamentals of continuum mechanics, in computational structure mechanics based on Lagrangian algorithms, each individual node of the computational mesh with its associated material particle flows during motion, which allows for easy tracking of moving boundaries, such as the free surfaces of water and interfaces between different materials. Also facilitated is the treatment of materials with history-dependent constitutive relations. Its weakness is its inability to follow large distortions of the computational domain without recourse to frequent remeshing operations. Compared with the Lagrangian algorithms, in CFD based on Eulerian algorithms, the computational mesh is fixed in the space, and the fluid moves with respect to the grid, in which large distortions in the continuum motion are relatively easily handled but generally at the expense of precise interface definition and the resolution of flow details. To overcome the shortcomings of purely Lagrangian and purely Eulerian descriptions, scientists and engineers have developed a technique to combine the best features of both approaches, which is known as the ALE description. In this ALE approach, the nodes of the computational mesh are either allowed to move within the continuum in the normal Lagrangian model, or are fixed in Eulerian manner, or are moved at some arbitrary specified velocities that provide a continuous rezoning capability. As a result, the computational mesh offered by an ALE description has freedom in moving during the numerical process, so that large distortions of the continuum can be handled much more easily than by a purely Lagrangian method and with higher resolution than that afforded by a purely Eulerian approach. The ALE coordinate system was originally proposed by Noh (1964), who adopted a finite difference (FD) scheme to solve 2D problems with moving boundaries in fluid dynamics, in which the mesh points of fluids are allowed to be in motion with the body or to be fixed in the space. Mathematically, as discussed in continuum mechanics by Truesdell (1991, 1966), this is essentially to arbitrarily choose a referential configuration to describe the motions of the continuum. Generally, the motions of mesh points on the referential configuration may be along only one direction with other directions fixed. For example, to deal with free surface water wave problems, the mesh points are allowed to move vertically but are fixed on the horizontal water level. Later, Trulio (1966) described the theory and structure of the ALE method used in the AFTON codes for the numerical solution of NavierStokes fluid equations, while Amsden and Hirt (1973) presented an ALE computer program for fluid flow at all speeds. Its further development by Hirt et al. (1974, 1997) is particularly noteworthy in its application to arbitrary finite difference meshes and permits flows at all speeds to be treated. Further refined numerical schemes of study to calculate 3D fluid flows around sharp interfaces were reported by Pracht (1975), Chan (1975), and Stein et al. (1977), and discussions of the analysis of free surface fluid flows and moving boundaries were presented by Ramaswamy and Kawahara (1986, 1987) and by Floryan and Rasmussen (1989). These contributions to the development of ALE techniques adopting FD techniques for the description of nonviscous fluid dynamics have produced powerful numerical tools (see, for example, Amsden and Hirt, 1973; Amsden et al., 1980, 1981). For viscous flows, Hughes et al. (1981) and Ramaswamy and Kawahara (1987), respectively, established ALE FE formulation for incompressible viscous flows and unsteady, convective, incompressible viscous surface fluid flows. Huerta and Liu (1988a,b, 1989), adopting the PetrovGalerkin FE approach based on an ALE description, presented excellent results for many engineering problems involving viscous flows with large free surface motions.

424

FluidSolid Interaction Dynamics

In nonlinear solid mechanics, contact problems involve the unknown contact areas between two solid bodies. To simulate this type of problem, FE approaches based on Lagrangian descriptions use very small and accurate meshes to keep the FE nodes of each body on the contact areas always in the contacted condition, which is very difficult due to large deformations, complex geometry, slippage, material nonlinearity, and different friction effects. The ALE method provides a very good way to address this difficulty, so that many publications on this subject have been reported. Haber and Abel (1983) noted the total displacement as a summation of the mesh displacement and the material displacement relative to the mesh point, kept the same mesh displacement at the contact point, but allowed different relative material displacements along the tangent direction of the contact interface. This technique essentially is an ALE approach, and it has been widely used to deal with friction contact problems (Haber and Abel, 1983; Haber and Koh, 1985; Haber and Hariandja, 1985; Haber, 1984); metal forming process (Huetink, 1982; Huetink et al., 1987, 1989, 1990; Schreurs et al., 1982, 1986; Liu et al., 1986, 1987, 1991), materials with memory and friction (Chen et al., 1989), large deformation analysis of elasticviscoplastic solids (Ghosh and Kikuchi, 1991), FEM of incompressible hyperelasticity (Yamada and Kikuchi, 1993), impact dynamics (Liu et al., 1988; Benson, 1987, 1989), and large sloshing (Liu and Huang, 1994). In fracture mechanics based on the ALE method, Haber and Koh (1985) allowed the crack length to be a continuous extension and derived explicit expressions for energy release rates, providing a very good stress strength factor (Haber, 1984). The ALE method essentially models engineering problems using a moving reference coordinate system in which the moving velocity is independent of the material motions and can be defined as arbitrary. This method combines the Lagrangian and Eulerian descriptions, overcoming their individual shortages, and therefore becomes a very important analysis method for nonlinear dynamics. Since the meshes move independently, the corresponding numerical process is more complex. To improve numerical efficiency, many publications have been reported, which intend to develop effective mesh motion algorithms, various solution strategies, reduced convective effects, etc.: for example, Ponthot (1988, 1989) for effective mesh management; Benson (1989) and Benson and David (1992) for vectorization techniques; Cheng and Kikuchi (1989) and Giuliani (1982) for mesh rezoning treatments; and the equipotential zoning techniques by Winslow (1963) and Winslow and Barton (1982). For nonlinear FSI problems in this book, in order to retain the validity of compatibility conditions on the coupling interfaces in nonlinear FSI problems, a lot of publications have investigated, developed, and discussed ALE numerical models using FEM both in the solid and the fluid domains: for example, Donea et al. (1977, 1980), Belytschko and Kennedy (1978), Chang et al. (1979), Kennedy and Belytschko (1981), Donea (1980, 1983), Hughes et al. (1981), Liu (1981a,b), Liu and Ma (1982a,b), Liu and Belytschko (1983), Liu and Gvildys (1986), Donea et al. (1982), Belytschko et al. (1980, 1982), Huerta and Liu (1988b), Ward et al. (1988), Nitikitpaiboon and Bathe (1993), Nomura and Hughes (1992), Nomura (1994), Bathe et al. (1995), and Wang (2000). In this model, the developed powerful techniques of FEM (see, e.g., Bathe, 1996; Zienkiewicz and Taylor, 1989, 1991) are used to describe the dynamics of the solid in the solid domain, but the benefit of adopting the powerful CFD to describe the dynamics of the fluid is reduced. To address this problem, Xing et al. (2002, 2003, 2005) using the ALE scheme developed a mixed FEFD method to simulate nonlinear FSI problems. The most important benefit of this mixed method is to combine the powerful FE software for solid structure analysis and the powerful FD software for fluid CFD simulation through the proposed partition iteration approach to solve nonlinear FSI problems. In principle, the idea of mixed FEFD numerical method for nonlinear FSI problems based on the partition iteration approach can be extended to construct various mixed schemes; for instance, various effective CFD methods may be chosen to match an FE or BE scheme in solids. Therefore the chapter title includes “mixed FECFD” approach to discuss other possible methods; it is not restricted to the FD formulations proposed previously. Toward this purpose, we discuss further some CFD schemes in ALE formulations to deal with nonlinear FSI problems. However, for readers, especially students who are not familiar with CFD methods, to better understand this method, we give some basic discretization techniques before tackling the FSI schemes. This fundamental knowledge can be found in many publications, such as Hirsch (1988, 1990) and Chung (2002).

11.2.2

Basic discretization techniques of computational fluid dynamics

11.2.2.1 Finite difference method Basic concept The FD approximation is the oldest method used to obtain numerical solutions of differential equations, and the idea corresponds to an estimation of a derivative by the ratio of two differences. For a function f ðxÞ, the derivative at point x is defined by

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

fx 

@f f ðx 1 ΔxÞ 2 f ðxÞ 5 lim ; Δx-0 @x Δx

425

(11.50a)

which is an approximation of the exact value of fx if Δx is small but finite. The error of this approximation is measured by the truncation error, which is dependent on Δx and goes to zero for Δx tending to zero. The power of Δx in the truncation error is called the order of the difference approximation, which can be determined by a Taylor series of f ðx 1 ΔxÞ at point x in the form f ðx 1 ΔxÞ 5 f ðxÞ 1 Δxfx ðxÞ 1

Δx2 fxx ðxÞ 1 ?; 2

(11.50b)

from which it follows f ðx 1 ΔxÞ 2 f ðxÞ Δx 5 fx ðxÞ 1 fxx ðxÞ 1 ? 5 fx ðxÞ 1 0ðΔxÞ: Δx 2

(11.50c)

This approximation for fx is said to be of the first order in Δx, indicating that the truncation error 0ðΔxÞ goes to zero like the first power of Δx. Finite difference approximations for first-order derivative Considering a 1D space with a space discretization represented by N discrete mesh points xI ; I 5 1; 2; . . .; N; in which at point xI ; the value of function f ðxI Þ 5 fI , the following FD approximations can be defined for the first derivatives ðfx ÞI  ð@f =@xÞx5xi . Forward difference ð f x ÞI 5

fI11 2 fI 1 0ðΔxÞ: Δx

(11.51a)

ð fx ÞI 5

fI 2 fI21 1 0ðΔxÞ: Δx

(11.51b)

fI11 2 fI21 1 0ðΔx2 Þ: 2Δx

(11.51c)

Backward difference

Central difference ð fx ÞI 5

The central difference has the second order of accuracy that can be verified by a Taylor expansion as follows: f ðx 2 ΔxÞ 5 f ðxÞ 2 Δxfx ðxÞ 1

Δx2 fxx ðxÞ 1 ?; 2

(11.52a)

from which, when this is subtracted from Eq. (11.50b), it follows f ðx 1 ΔxÞ 2 f ðx 2 ΔxÞ 5 fx ðxÞ 1 0ðΔx2 Þ; 2Δx

(11.52b)

that is, Eq. (11.51c). Based on the central difference, the forward and backward formulations can be considered a corresponding central difference with respect to the middle point: xI11=2 5

xI 1 xI11 ; 2

xI21=2 5

xI21 1 xI ; 2

(11.53a)

so that we have the second-order approximations for the middle points in Eq. (11.53a) in the forms fI11 2 fI 1 0ðΔx2 Þ; 2Δx fI 2 fI21 1 0ðΔx2 Þ; ð fx ÞI21=2 5 2Δx

ð fx ÞI11=2 5

(11.53b)

but involving only the same mesh points I and (I 1 1). This is an important property that is often used in computations due to its compactness.

426

FluidSolid Interaction Dynamics

Finite difference approximations for second-order derivative The FD approximation of higher-order derivatives can be obtained by repeated application of first-order formulas. For example, a second-order approximation to the second derivative ð fxx ÞI  ð@2 f =@x2 Þx5xi is obtained by 0 1   @2 f ð fx ÞI11 2 ð fx ÞI ð f ÞI11 2 2ð f ÞI 1 ð f ÞI21 5 ð fxx ÞI  @ 2 A 5 1 0 Δx2 ; (11.54a) 2 @x Δx Δx I

whereas backward approximations are used for the first-order derivatives, that is, ð fx ÞI11 5

ð f ÞI11 2 ð f ÞI ; Δx

ðfx ÞI 5

ð f ÞI 2 ð f ÞI21 : Δx

(11.54b)

Using a Taylor expansion, we obtain   ð f ÞI11 2 2ð f ÞI 1 ð f ÞI21 Δx2 @4 f 5 ð f Þ 1 1 ?; xx I Δx2 12 @x4

(11.54c)

which confirms the second-order accuracy for the second-order derivative approximation in Eq. (11.54a). Difference formulas with multiple points A difference formula with multiple points can be constructed by using the Taylor expansion approach. For example, a one-sided, second-order difference formula for the first-order derivative ( fx)I involving three upstream points I2, I1, and I can be assumed in the form ð fx ÞI 5

αfI 1 βfI21 1 γfI22 1 0ðΔx2 Þ; Δx

(11.55a)

where the coefficients α, β, and γ are sought from the following Taylor expansions: fI21 5 fI 2 Δxð fx ÞI 1

Δx2 Δx3 ðfxx ÞI 2 ð fxxx ÞI 1 ?; 2 6

fI22 5 fI 2 2Δxð fx ÞI 1 2Δx2 ð fxx ÞI 2

ð2ΔxÞ3 ðfxxx ÞI 1 ?; 6

(11.55b)

from which it follows αfI 1 βfI21 1 γfI22 5 ðα 1 β 1 γÞfI 2 Δxðβ 1 2γÞðfx ÞI 1

Δx2 ðβ 1 4γÞðfxx ÞI 1 0ðΔx3 Þ: 2

(11.55c)

Therefore, identifying with Eq. (11.55a), we obtain the three conditions α 1 β 1 γ 5 0;

β 1 2γ 5 2 1;

β 1 4γ 5 0;

(11.55d)

and Eq. (11.55a) becomes ð f x ÞI 5

3fI 1 4fI21 1 fI22 1 0ðΔx2 Þ: 2Δx

(11.55e)

General operator methods for finite difference formulas A general theory based on operator definitions can be read in the book by Hildebrand (1956), in which the following operators are defined: E21 fI 5 fI21 :

Displacement operator E: EfI 5 fI11 ; 1

Forward difference operator δ : 2

1

δ fI 5 fI11 2 fI : 2

(11.56a) (11.56b)

Backward difference operator δ : δ fI 5 fI 2 fI21 :

(11.56c)

Central difference operator δ: δfI 5 fI11=2 2 fI21=2

(11.56d)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

ðfI11 2 fI21 Þ : Central difference operator δ: δfI 5 2   fI11=2 2 fI21=2 Averaging operator μ: μfI 5 : 2 @f Differential operator D: Df 5 fx  : @x

427

(11.56e) (11.56f) (11.56g)

From these definitions, we have the following obvious relations: δ1 5 E 2 1;

δ2 5 1 2 E21 5 E21 δ1 ;

δ1 δ2 5 δ2 δ1 5 δ1 2 δ2 5 δ 2 ; ðE 2 E21 Þ ; δ 5 E1=2 2 E21=2 5 2

μ5

ðE1=2 1 E21=2 Þ ; 2

μ2 5 1 1 δ2 =4;

(11.57)

δ2 5 δ1 δ1 5 E2 2 2E 1 1;

δ3 5 ðE21Þ3 5 E3 2 3E2 1 3E 2 1: Based on these operators, we can derive the key to the operator technique for generating FD formulas, which is the relation between the derivative operator D and the displacement operator E and which can be obtained from the Taylor expansion: # " ðΔxDÞ2 ðΔxDÞ3 1 1 ? f ðxÞ; Ef ðxÞ 5 1 1 ΔxD 1 2! 3! (11.58) Ef ðxÞ 5 eΔxD f ðxÞ;

E 5 eΔxD ;

ΔxD 5 ln E:

As a result, we can obtain the FD formulas as follows. Forward differences ΔxD 5 ln E 5 lnð1 1 δ1 Þ 5 δ1 2

δ12 δ13 δ14 1 2 1? 2 3 4

(11.59a)

The order of accuracy of the approximation increases with the number of terms kept on the right-hand side of Eq. (11.59a), and the first neglected term gives the truncation error. For example, if the first two terms are kept, we obtain the second-order formula ð fx ÞI 5 DfI 5

ðδ1 2 δ12 =2ÞfI 0ðδ13 ÞfI 1 Δx 3Δx

5

2δ1 fI 2 ðE2 2 2E 1 1ÞfI Δx2 1 fxxx 2Δx 3

5

2fI11 2 2fI 2 ð fI12 2 2fI11 1 fI Þ Δx2 1 fxxx 2Δx 3

5

2 3fI 1 4fI11 2 fI12 Δx2 1 fxxx ; 2Δx 3

where the last term on the right-hand side gives the truncation error. Backward difference δ22 δ23 δ24 1 1 1? ΔxD 5 ln E 5 2 lnð1 2 δ2 Þ 5 δ2 1 2 3 4 Central difference δfI 5 fI11=2 2 fI21=2 5 ðE ΔxD 5 2 sin h21 δ 5 δ 2

1=2

2E

21=2

Þ fI 5 ðe

δ3 3δ5 1 1? 24 640

ΔxD=2

2ΔxD=2

2e

(11.59b)

(11.60)

  ΔxD Þ fI 5 2 sin h fI ; 2 (11.61)

428

FluidSolid Interaction Dynamics

Higher-order derivatives The operator approach can be used to generate FD formulas for higher-order derivatives. Based on Eqs. (11.59) (11.61), we obtain the following formulas for the nth-order derivatives; see for instance, Ames (1977) and Hirsch (1988). Forward difference  n  @ f 1 5 Dn fI 5 n ½ lnð11δ1 Þn fI n @x I Δx # " (11.62a) 1ðn11Þ 1ðn12Þ 1ðn13Þ 1 nδ nð3n 1 5Þδ nðn 1 2Þðn 1 3Þδ 1n 1 2 1 ? fI 5 n δ 2 Δx 2 24 48 Backward difference  n  @ f ð21Þn n 5 D f 5 ½ lnð12δ1 Þn fI I @xn I Δxn # " 1 nδ1ðn11Þ nð3n 1 5Þδ1ðn12Þ nðn 1 2Þðn 1 3Þδ1ðn13Þ 1n 1 2 1 ? fI 5 n δ 2 Δx 2 24 48 Central difference

0 1n 2 3 1 n 2 4 n @ f 2 δ δ nδ nð5n 1 22Þδ @ A 5 Dn fI 5 @ sin h21 A fI 5 n 41 2 1 1 ?5fI @xn Δx 2 Δx 24 64 3 90

(11.62b)

0

(11.62c)

I

0 1n 2 3 n 2 4 n 2 @ f μ 2 δ μδ ðn 1 3Þδ ð5n 1 52n 1 135Þδ @ @ A 5 D n fI 5 sin h21 A fI 5 n 41 2 1 1 ?5fI @xn 2 Δx 24 5760 ð110:25δ2 Þ1=2 Δx 0

1

(11.62d)

I

Here, Eq. (11.62c) generates FD formulas with the function values at the integer mesh points for n even, but n uneven involves half-integer mesh points. However, Eq. (11.62d) provides the reverse results. From these general formulas, we can obtain the corresponding FD formula for a derivative order n. For example, second-order formulas, n 5 2, are 1 11δ14 5δ15 ð fxx ÞI 5 2 δ12 2 δ13 1 2 1 ? fI ; (11.63a) Δx 12 6 1 11δ14 5δ15 ð fxx ÞI 5 2 δ12 1 δ13 1 1 1 ? fI ; (11.63b) Δx 12 6 1 δ4 δ6 δ8 ð fxx ÞI 5 2 δ2 2 1 2 1 ? fI ; (11.63c) Δx 12 90 560 1 5δ4 259δ6 ð fxx ÞI 5 2 δ2 2 1 1 ? fI : (11.63d) Δx 24 5760 These equations define four families of difference operators for the second-order derivatives to various orders of accuracy. If we keep the first two terms, we obtain the following FD formulas: Forward: second-order accuracy ð fxx ÞI 5

1 11 ð2fI 2 5fI11 1 4fI12 2 fI13 Þ 1 Δx2 fxxxx : Δx2 12

(11.64a)

1 11 ð2fI 2 5fI21 1 4fI22 2 fI23 Þ 2 Δx2 fxxxx : 2 Δx 12

(11.64b)

Backward: second-order accuracy ð fxx ÞI 5

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

429

Central: integer points, fourth-order accuracy ð fxx ÞI 5

1 Δx4 @6 f ð2 fI12 1 16fI11 2 30fI 1 16fI21 2 fI22 Þ 1 : 2 12Δx 90 @x6

(11.64c)

Central: half-integer points, fourth-order accuracy ð fxx ÞI 5

 259Δx4 @6 f 1  : 25fI15=2 1 39fI13=2 2 34fI11=2 1 34fI21=2 1 39fI23=2 2 5fI15=2 1 2 48Δx 57; 600 @x6

(11.64d)

Implicit finite difference formulas and its general generation In the previous FD formulations, the derivatives of the function at a point I are directly calculated by the values of the function at several mesh points and are considered explicit FD formulas. For a higher-order accuracy, when the derivatives at difference mesh points are related to one another, the derivatives at different mesh points appear simultaneously in the generated FD formulas, so that the derivatives at mesh points cannot be represented explicitly based on the function values at the mesh points. This type of FD formula is considered an implicit FD formula, such as an implicit upwind algorithm for computing turbulent flows on unstructured grids by Anderson and Bonhaus (1994). The general generation method of an implicit FD formula is similar the one discussed for Eqs. (11.55a)(11.55e) but needs to include the derivatives involved. Detailed information can be found in Peyret and Taylor (1982). Here we discuss an example to illustrate this method. Consider a two-point expression in the form α1 fI11 1 α0 fI 1 β 1 ð fx ÞI11 1 β 0 ð fx ÞI 1 γ 1 ð fxx ÞI11 1 γ 0 ð fxx ÞI 5 0:

(11.65a)

Developing all the variables in a Taylor series about point I, we have fI11 5 fI 1 Δxð fx ÞI 1

Δx2 Δx3 ð fxx ÞI 1 ð fxxx ÞI 1 ?; 2 6

ð fx ÞI11 5 ð fx ÞI 1 Δxð fxx ÞI 1

Δx2 ð fxxx ÞI 1 ?; 2

ð fxx ÞI11 5 ð fxx ÞI 1 Δxð fxxx ÞI 1

(11.65b)

Δx2 ð fxxxx ÞI 1 ?; 2

which, when substituted into Eq. (11.65a) with at least second-order accuracy for the second derivative, gives     α1 Δx2 α1 Δx3 β Δx2 1 γ 0 1 γ 1 1 β 1 Δx ð fxx ÞI 1 1 γ 1 Δx 1 1 ðα0 1 α1 Þ fI 1 ðα1 Δx 1 β 0 1 β 1 Þð fx ÞI 1 ð fxxx ÞI 5 0; 2 6 2   α1 Δx4 γ 1 Δx2 R5 1 ð fxxxx ÞI : 12 2 (11.65c) Letting the coefficients of Eq. (11.65c) vanish, we obtain α1 5 2 α0 5 1;

β0 5 β1 5

2 Δx ; 2

γ1 5 2 γ0 5

Δx2 ; 12

(11.65d)

so that we further obtain a two-point implicit FD formula with 1 1 1  ð fI11 2 fI Þ 2 ð fx ÞI11 1 ð fx ÞI 1 ð fxx ÞI11 2 ð fxx ÞI 5 0; Δx2 2Δx 12

(11.65e)

with its truncation error R5

Δx2 ð fxxxx Þ: 8

(11.65f)

430

FluidSolid Interaction Dynamics

Multidimensional finite difference formulas The partial derivatives of functions of several variables can be approximated by considering each variable separately. For example, in a 2D case, we define fIJ 5 f ðxI ; yJ Þ, indicate the FD operator in the space coordinate by subscript x or y, and have the following FD formulas with first-order accuracy:   @f fI11;J 2 fI;J 1 1 δ fI;J 1 0ðΔxÞ; 5 5 ð fx ÞIJ 5 @x IJ Δx x Δx   @f fI;J11 2 fI;J 1 1 δ fI;J 1 0ðΔyÞ; ð fy ÞIJ 5 5 5 (11.66) @y IJ Δy y Δy  2  @f fI11;J 2 2fI;J 1 fI21;J ð fxx ÞIJ 5 5 1 0ðΔx2 Þ: @x2 IJ Δx2 Nonuniform meshes For nonuniform meshes, the discretization of the equation can be performed after a transformation from the physical space (xi) to a Cartesian computational space (ξi ), which is defined by ξ i 5 ξ i ðxj Þ. Here, subscripts are the tensor index type used in this book. This transformation must be one to one and reversible to keep a point in the computational space corresponding to only one point in the physical space. All the formulas previously derived can be applied in the computational space. Fig. 11.3 shows a transformation for 2D case.

11.2.2.2 Finite volume method The finite volume method (FVM) is based on the integral form of the conservation laws of continuum mechanics discussed in Chapter 3, Fundamentals of continuum mechanics. This method was introduced and developed by McDonald (1971), MacCormack and Paullay (1972), and Rizzi and Inouye (1973). The method takes an arbitrary mesh with a large number of options to define the control volumes for the conservation laws, giving considerable flexibility to modify the shape and location of the control volumes, as well as to evaluate the fluxes through the control surfaces for computational accuracy. The direct discretization of the integral form of the conservation laws can ensure that the basic laws also remain conserved at the discrete level, which is a most fundamental property for numerical schemes. Conservative discretization The general integral form of the conservation laws in the Eulerian coordinate system is given by Eq. (3.83b), that is, ð ð ð @ðρFÞ 1 ðρFvÞUr dΩf 1 qUνdΓ f 5 rρdΩf ; (11.67a) @t Ωf Γf Ωf which, for the numerical implementation in CFD, can be rewritten as ð ð ð @ ρFdΩf 1 ðq 1 ρFvÞUνdΓ f 5 rρdΩf ; @t Ωf Γf Ωf

(11.67b)

where the Green theorem has been introduced to transform the volume integration of the second term on the left-hand side into its surface integration, and the partial derivative with respect to time t is moved to the front of the integration x2

s ξ2

Line

ξ2 Lines ξ1

x1

ξ1

FIGURE 11.3 Transformation between the physical space (left) and the computational space (right) in a 2D case.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

431

notation due to its local time derivative at a fixed space point. Using the notations given in Chapter 3, Fundamentals of continuum mechanics, Eq. (11.67b) can be recast into the so-called conservation form of the NavierStokes system of equations covering mass, momentum, and energy conservation: ð ð ð @ _ UdΩf 1 ðQi 1 Gi Þν i dΓ f 5 BdΩf ; (11.67c) @t Ωf Γf Ωf where

2

ρf

3

2

ρ f vi

3

6 7 7 ρ f vj ρf vj vi 1 pδij 5; 5; Qi 5 4 ρf ðe 1 vj vj =2Þ ρf ðe 1 vj vj =2Þvi 1 pvi 2 3 2 3 0 0 6 7 6 7 ρf fj 2 σ~ ij Gi 5 4 5; B 5 4 5: ρf ð fj vj 1 θÞ 2vj σ~ ij 1 hi

_ 6 U54

(11.67d)

In this matrix equation, the first line represents the mass conservation law, while the second and third lines denote _ the momentum and energy conservation laws, respectively. The time variation of the variable vector U inside the volume depends only on the surface values of the fluxes. Hence for an arbitrary subdivision of the volume Ωf , we can write the conservation law for each subdomain and obtain the global conservation law by a summation of all subdomain conservation laws. As shown in Fig. 11.4, the volume Ωf is divided into two subdomains Ωð1Þ and Ωð2Þ f f , for each of which Eq. (11.67c) becomes ð ð ð @ _ ðIÞ ðIÞ ðIÞ UdΩf 1 ðQi 1 Gi Þν i dΓ f 5 BdΩðIÞ (11.68) f ; ðIÞ ðIÞ @t ΩðIÞ Γf Ωf f where I 5 1,2. In this equation, the normal vector ν ðIÞ i is defined with its positive direction pointing from inside to outside ð2Þ the domain on the common boundary ν ð1Þ i 5 2 ν i , as shown in Fig. 11.4. As a result, the summation of Eq. (11.68) for the two subdomains should cancel the surface flux integrations on their common boundary. This essential property has to be satisfied by the numerical discretization of the flux contributions in order for a scheme to be conservative. For a discretization, if the resulting equation still contains flux contributions from inside to the total cell, this discretization is said to be nonconservative, and the internal flux contributions appear as numerical internal volume sources. In computational mathematics, the LaxWendroff theorem (Lax and Wendroff, 1960) guarantees that if a conservative numerical scheme for a hyperbolic system of conservation laws converges, then it converges toward a weak solution of the original equation. Fundamental concepts of finite volume method Discrete equation and its constrains _ ðIÞ For a control volume ΩðIÞ f and the associated variable U , Eq. (11.68) can be rewritten by the discrete form @ _ ðIÞ ðIÞ  X ðIÞ ðIÞ ðIÞ U Ωf 1 ðQi 1Gi ÞðIÞ ν ðIÞ (11.69) i Γ f 5 B Ωf ; @t ðIÞ Γf

ð1Þ ð2Þ FIGURE 11.4 Volume Ωf is divided into the two subdomains, Ωfð1Þ and Ωð2Þ f , whose the common boundary is the outside normal ν i 5 2 ν i .

432

FluidSolid Interaction Dynamics

where the sum of the flux terms refers to all the external sides of the control cell ΩðIÞ f . The following constraints on the chosen ΩðIÞ volumes for a conservative FVM have to be satisfied: f 1. Their sum should cover the whole domain Ωf : ðIÞ 2. Adjacent ΩðIÞ f may overlap if each internal surface Γ f is common to the two volumes. 3. Fluxes along a cell surface have to be computed by formulas independent of the cell in which they are considered. This ensures that the conservative property is satisfied, since the flux contributions of internal boundaries will cancel when the contributions of the associated finite volumes are added. Features of finite volume method equation The FVM Eq. (11.69) has the following features distinguishing it from FD and FEM: _ ðIÞ

1. Point I locations in the control volume ΩðIÞ f do not appear explicitly, so that U f can be considered an average value _ of the flow variable U f over the control cell, and the first term in Eq. (11.69) represents the time change rate of the averaged flow variable over the control cell. ðIÞ 2. The mesh coordinates appear only in the calculations of cell volume ΩðIÞ f and boundary areas Γ f . 3. If no source term on the right-hand side of Eq. (11.69), the FVM equation implies that the variation of the averaged _ ðIÞ value U f of the variable over a time interval Δt is balanced by the sum of the fluxes exchanged between neighboring cells. It is convenient to program the method by sweeping through the cell surface to calculate the flux on the common boundary of two adjacent cells, adding this flux to one cell but subtracting it from the other. This automatically guarantees global conservation. 4. FVM allows a natural introduction of boundary conditions defined by the physical problem, such as some normal components vanishing on the solid boundary. Mesh and control volume definitions FVM can handle any type of mesh and has the same flexibility as FEM. Two types of meshes can be considered. 1. An FD-type mesh, in which all mesh points lie on the intersection of two or three families of lines. Generally these lines are considered curvilinear coordinate lines and are designated as structural meshes. 2. An FE-type mesh, in which the mesh points cannot be identified with coordinate lines in general cases. Therefore they cannot be represented by a set of integers but have to be numbered individually in a certain order. This type of mesh is designated as unstructured. For the selected mesh, there are two ways to define the variables: 1. A cell-centered FVM, which sets the variables associated with a cell with averaged values over the cell, which are considered representative of some point inside the cell, such as its central point. In this case, the mesh cell is just the control volume. 2. A cell-vertex FVM, which attaches the variables to the mesh points, that is, the cell vertices. This case provides greater flexibility in defining the control volume around a mesh point; for example, it is very often used to select a control volume covered by half-mesh-cell surfaces around a vertex. Evaluation formulas for finite volume method It is often necessary in FVM to complete the evaluations involving the mesh geometries and variables, of which the formulas are as follows. General integration formulas In the NavierStokes equations for viscous flows, the viscous flux components are functions of the velocity gradients, for which we have to estimate approximate values on the cell faces. A general procedure can be derived by the Green theorem; for example, for a scalar function U in an arbitrary volume Ω of boundary Γ , the Green theorem confirms ð I U;i dΩ 5 Uν i dΓ ; (11.70a) Ω

Γ

from which, the averaged gradients can be defined as  a ð I @U 1 1 1 X ðJÞ  U;i dΩ 5 Uν i dΓ 5 U ðν i Γ ÞðJÞ ; @xi Ω Ω Ω Ω Γ Ω J where superindex J identifies the boundary face with the related variables of the volume Ω.

(11.70b)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

433

Evaluation of cell face areas In Eq. (11.70b), the positive direction of the unit normal vi points from inside to outside the volume. Therefore the , cell face area ðν i Γ ÞðJÞ is a vector shown by Fig. 11.5. The face area vector Γ ABCD can be evaluated by using the vector cross-multiplication formula given by Eq. (2.10) in Chapter 2, Cartesian tensor and matrix calculus, that is, , 1 Γ ABCD 5 AC 3 BD 5 ðAB 3 BC 1 CD 3 DAÞ; 2

which can be written in the tensor form, such as ,  Γ ABCD 5 eijk ðACÞj ðBDÞk ; ðACÞj 5 xCj 2 xAj ; i

B ðBDÞk 5 xD k 2 xk :

(11.71a)

(11.71b)

Here, the superindex identifies the vertex point coordinates. Substituting U 5 1 into Eq. (11.70a), we obtain I Uν i dΓ 5 0; (11.71c) Γ

which may be used to check the result evaluated by Eqs. (11.71a)(11.71c). Evaluation of control cell volumes Replacing U by xi in Eq. (11.70a), we obtain ð ð I xi;i dΩ 5 3dΩ 5 3Ω 5 xi ν i dΓ ; Ω

Ω

(11.72a)

Γ

so that the control cell volume is evaluated by I 1 1 X ðJÞ 1 X ðJÞ ðJÞ Ω5 xi ν i dΓ 5 xi ðν i Γ ÞðJÞ 5 x Γ ; 3 Γ 3 J 3 J i i

(11.72b)

,

with Γ ðJÞ i 5 Γ i ðJÞ determined by Eq. (11.71b) for face J. Evaluation of fluxes through cell faces In Eq. (11.69), the second term on the left-hand side defines the fluxes through cell faces and involves the value of cell variables on the cell faces. The values of variables on the cell faces depend on the discretization scheme adopted. G

G

For cell-centered FVM, the value of a variable is defined at the cell center, so that the averaged values of variable values on the two adjacent cells are chosen to evaluate the fluxes through cell faces. For cell-vertex FVM, the variable’s value is defined at each vertex of the cell face, so that the averaged value of the variables at all vertices of a cell can be used to evaluate the fluxes through this cell face.

ΓDCGH G

H

C

ΓEAGH D F

ΓBFGC

E

B ΓABCD A

FIGURE 11.5 Three-dimensional hexahedral control volume with its face area vectors.

434

FluidSolid Interaction Dynamics

11.2.3

Consistency, stability, and convergence of numerical schemes

For numerical schemes that can be adopted to simulate practical problems governed by the partial differential equations, we have to impose some additional conditions to guarantee that: 1. the discretized equation of the scheme provides an acceptable approximation to the original differential problems; 2. the numerical solution obtained in each step can stably approach the exact solution of the discrete equation; and 3. the numerical solution approaches the exact solution of the differential equation. These three requirements are defined as consistency, stability, and convergence conditions, respectively, and they cover different aspects of the relations between the discretized equations, the numerical solution, and the exact analytical solution of the differential equation. The consistency condition for condition (1) defines a relation between the differential equation and its discrete formulation, which concerns the structure of numerical formulation; the stability condition (2) establishes a relation between the computed solution and the exact solution of the discretized equations, while the convergence condition (3) connects the computed solution to the exact solution of the differential equation. We illustrate these three concepts by an example equation called a linear convection equation with its initial and boundary conditions for a . 0 @u @u 1 a 5 0; @t @x uðx; 0Þ 5 f ðxÞ; 0 # x # L; uð0; tÞ 5 gðtÞ; t $ 0;

(11.73)

which, if the variable u 5 ρ (mass density) and a 5 v (flow velocity), describes the transport of mass ρ by a flow of velocity v. Also, this equation is a typical first-order hyperbolic equation in (x,t) being viewed to describe a wave of amplitude u propagating with a wave speed a.

11.2.3.1 Consistency We choose a central difference formula for the discretization of the space derivative ux and a forward difference formula for ðux ÞI , so that we obtain the discretization equations ðut ÞI 5

 un11 2 unI 2 a n I uI11 2 unI21 ; 5 2Δx Δt

2a ðuI11 2 uI21 Þ; 2Δx

(11.74a)

where n indicates the time step, and unI 5 uðxI ; nΔtÞ: This is an explicit scheme because only one unknown at level (n 1 1) is included. Based on the following Taylor expansions: 5 unI 1 Δtðut ÞnI 1 un11 I

Δt2 ðutt ÞnI 1 ?; 2

unI11 5 unI 1 Δxðux ÞnI 1

Δx2 Δx3 ðuxx ÞnI 1 ðuxxx ÞnI 1 ?; 2 6

unI21 5 unI 2 Δxðux ÞnI 1

Δx2 Δx3 ðuxx ÞnI 2 ðuxxx ÞnI 1 ?; 2 6

(11.74b)

substituted into Eq. (11.74a), we can obtain a truncation error εT between the values of the differential equation and the discretized equation and at time step n and at space point I in the form 2 3 n11 n   u 2 u a 2 Δt Δx2 I unI11 2 unI21 5 5 ðutt ÞnI 2 εT 5 ðut 1aux ÞnI 2 4 I 1 aðuxxx ÞnI 1 0ðΔt2 ; Δx4 Þ: (11.74c) 2Δx 2 Δt 6 It confirms that this truncation error vanishes when Δt and Δx tend to zero; therefore the scheme in Eq. (11.74a) is consistent. The accuracy of the scheme is the first order in time and second order on space. The condition for consistency may be restated: a scheme is consistent if the truncation error tends to zero for Δt and Δx tending to zero.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

435

11.2.3.2 Stability The stability condition is a requirement solely on the numerical scheme and does not involve the differential equation. If we define the error ε as the difference between the computed solution u and the exact solution u of the discretized equation, that is, εnI 5 unI 2 unI ; the stability condition can be written as



lim εnI # K;

n-N

(11.75a)

at fxed Δt;

(11.75b)

where K is a positive constant and independent of n. This defined condition does not ensure that the error will not become unacceptably large at fixed intermediate times tn 5 nΔt; so that a more general definition of stability is introduced by Lax and Richtmyer (1956) and further developed by Ritchtmyer and Morton (1967), which is based on the time behavior of the solution itself instead of the error behavior. It states that any component of the initial solution should not be amplified without bound, which can be represented as :An : # K;

0 , Δt , τ;

0 # nΔt # T;

(11.76a)

for fixed values τ and T and for all n. Here, A is an operator, a function of time step Δt and mesh size Δx, which connects the values of variable vectors Un and Un11 by the transformation Un11 5 AUn ;

n 5 0; 1; 2; . . .:

(11.76b)

11.2.3.3 Convergence ~ tÞ of the differential The convergence condition requires the numerical solution unI approaching the exact solution uðx; equation at any point xI 5 IΔx and time tn 5 nΔt when Δx and Δt tend to zero. This can be represented as

n

~ ε~ nI 5 unI 2 uðIΔx; nΔtÞ; lim

ε~ 5 0; at fixed xI 5 IΔx and tn 5 nΔt: (11.77a) Δt-0;Δx-0 I Using the operator developed by Ritchtmyer and Morton (1967), the convergence condition can also be rewritten as lim

Δt-0;Δx-0

~ :½AðΔtÞn Uð0Þ 2 UðtÞ: 5 0;

at fixed t 5 nΔt:

(11.77b)

The conditions of consistency, stability, and convergence are related to one another, and the precise relation is given by the fundamental Equivalence Theorem of Lax, a proof that can be found in Ritchtmyer and Morton (1967). The theorem states that, for a well-posed initial value problem and a consistent discretization scheme, stability is the necessary and sufficient condition for convergence. This theorem shows that in order to numerically solve an initial value problem, two tasks have to be performed: (1) analyze the consistency condition to determine the order of accuracy of the scheme and its truncation error, and (2) investigate the stability of the scheme. From these two steps, convergence can be defined with no additional analysis.

11.2.3.4 von Neumann method for stability analysis Various methods have been developed for the analysis of stability, but nearly all of them are limited to linear problems, in which the von Neumann method (see, e.g., Cranck and Nicholson, 1947; Charney et al., 1950; Hirsch, 1988) is the most widely applied technique for stability analysis, allowing an extensive investigation of the behavior of the error as a function of the frequency content of the initial data and of the solution. Here, we illustrate this method to investigate the stability of Eq. (11.74a) as follows. Error equation and its Fourier series As mentioned, in Eq. (11.75a) the error εnI is the difference between the computed solution unI and the exact solution unI of the discretized equation. Substituting Eq. (11.75a) into Eq. (11.74a), we obtain   un11 2 unI εn11 2 εnI 2a  n a  n I εI11 2 εnI21 ; 1 I 5 uI11 2 unI21 2 2Δx 2Δx Δt Δt

(11.78a)

436

FluidSolid Interaction Dynamics

from which, when considering unI satisfying Eq. (11.73) exactly, it follows  εn11 2 εnI a  n I εI11 2 εnI21 : 52 2Δx Δt

(11.78b)

Therefore the errors εnI evolve over time in the same way as the numerical solution unI . Introducing an error vector ε consisting of the errors at all space points as we define the variable vector U for Eq. (11.76b), we have εn11 5 Aεn ;

n 5 0; 1; 2; . . .:

(11.79)

As shown in Fig. 11.6, we consider the domain (0,L) of the problem in Eq. (11.73) as the half period of the error ε, so that this error εnI can be represented by the following Fourier series: εnI 5

N X

EnJ ejκJ IΔx 5

J52N

L Δx 5 ; N

N X

EnJ e jIφJ ;

φJ 5 κJ Δx 5

J52N

xI 5 IΔx;

κJ 5 Jκmin ;

κ5

2π ; λ

Jπ ; N

I 5 0; 1; 2; . . .; N; λmax 5 2L;

(11.80)

κmin 5

2π : λmax

Here, λ and κ represent the wavelength and the wavenumber, respectively, of which the maximum wavelength is λmax 5 2L corresponding to the minimum wavenumber κmin 5 π=L; EJ and φJ , respectively, are the amplitude and the phase angle of the Jth harmonic component of wavenumber κJ 5 Jπ=L. This Fourier series reveals the frequency components of the error, of which the region around φ 5 0 corresponds to the low frequencies, while the region close to φ 5 π is associated with the high-frequency range of the spectrum. The value of φ 5 π gives the highest frequency resolvable on the mesh, that is, the frequency of the wavelength 2Δx. Amplification factor For a single harmonic component EnJ e jIφJ , the error Eq. (11.78b), dropping the subscript J, becomes En11 2 En jIφ aEn  jðI11Þφ e 1 e 2 e jðI21Þφ 5 0; Δt 2Δx

(11.81a)

σEn jφ ½e 2 e2jφ  5 0; 2

(11.81b)

which, when divided by e jIφ , gives En11 2 En 1

σ5

aΔt : Δx

Here, the parameter σ is called the Courant number, which physically represents the ratio of the distance that a numerical disturbance is transmitted in the time Δt over the mesh size Δx. The stability condition in Eq. (11.77a) is satisfied if the amplitude of any error component En does not increase in time, that is,

n11

E

jAj 5

n

# 1; for all φ: (11.81c) E

Error ε

x –L

0

L

FIGURE 11.6 Fourier representation of the error on the domain (2L, L): the maximum wavelength λmax 5 2L with the corresponding minimum wavenumber κmin 5 π=L.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

437

The quantity A is defined as the amplification factor of the error, Aðφ; Δt; ΔxÞ 5

En11 ; En

(11.81d)

which is a function of time step Δt, mesh size Δx, and frequency. Unconditionally unstable example For the scheme defined by Eq. (11.74a), Eq. (11.81b) becomes A 2 1 1 jσ sin φ 5 0: The satiability condition in Eq. (11.81c) requires the modulus of the amplification factor A qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jAj 5 1 1 σ2 sin2 φ

(11.81e)

(11.81f)

is not larger than 1, which is clearly never satisfied. Hence, the centered scheme in Eq. (11.74a) for the convection equation with forward difference in time is unconditionally unstable. Conditional stability example If we adopt the forward difference in time and the explicit first-order upwind scheme (backward difference) in space to solve the convection Eq. (11.73), the numerical equation is  un11 2 unI 2 a n I uI 2 unI21 ; 5 Δx Δt for which the error equation is given by  n11   E 2 En e jIφ 1 σEn e jIφ 2 e jðI21Þφ 5 0;

(11.82a)

(11.82b)

Or, after division by En e jIφ , A 5 1 2 σ 1 σe2jφ 5 1 2 σ 1 σ cos φ 2 jσ sin φ:

(11.82c)

The real and imaginary parts of the amplification factor are AR 5 1 2 σ 1 σ cos φ; AI 5 2 σ sin φ; which geometrically describe a circle of radius σ with its center at ð1 2 σ; 0Þ in the complex plane, that is, 2 AR 2ð12σÞ 2 AI 1 5 1: σ σ

(11.83d)

(11.83e)

Fig. 11.7 shows this circle, from which it is clearly seen that the scheme is stable for 0 , σ # 1;

(11.83f)

so that the scheme in Eq. (11.82a) is conditionally stable. This condition is called as the CourantFriedrichsLewy condition (Courant et al., 1928 or Courant et al., 1967). Unconditional stability example Considering the implicit scheme to solve the convection equation in Eq. (11.73), the numerical equation is defined by  un11 2 unI 2 a  n11 I uI11 2 un11 5 I21 ; 2Δx Δt

(11.84a)

in which, for the time step n 1 1, the backward difference in time is adopted. The error equation of this scheme is  εn11 2 εnI a  n11 I εI11 2 εn11 52 I21 : 2Δx Δt

(11.84b)

438

FluidSolid Interaction Dynamics

A1 Region of instability A 1 σ

φ

1–σ

AR

FIGURE 11.7 Geometric representation of the complex amplification factor A in the complex place.

After the introduction of a harmonic component in the form En ejIφ , Eq. (11.84b) gives ðEn11 2 En Þe jIφ 1 σEn11 ½e jφ 2 e2jφ e jIφ =2 5 0;

(11.84c)

pffiffiffi A 2 1 1 jσA 2 sin φ 5 0;

(11.84d)

or

The modulus of the amplification factor A is

A5

1 : ð1 1 jσ sin φÞ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð1 1 jσ sin φÞð1 2 jσ sin φÞ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 5 , 1; 1 1 σ2 sin2 φ

pffiffiffiffiffiffiffiffi jAj 5 AA 5

(11.84e)

for all values of σ, and therefore the implicit scheme in Eq. (11.84a) is unconditionally stable. General von Neumann’s formation Referring to Eqs. (11.76a) and (11.76b), we consider it a time integration scheme, that is, Un11 5 AðφÞUn ;

n 5 0; 1; 2; . . .;

(11.85a)

in which the vector Un denotes the value of variable at time step n, and the matrix AðφÞ is a discretization operator, a function of phase angle φ to identify the scheme. We also call this matrix an amplification matrix. The stability condition in Eq. (11.76a) requires that the matrix ½AðφÞn remains uniformly bounded for all values of φ. The matrix A is assumed to be an N 3 N matrix with N eigenvalues λI ðI 5 1; 2; 3; . . .; NÞ, defined by the eigenvalue equation AϕI 5 λI ϕI ;

detjA 2 λIj 5 0;

(11.85b)

where I is an N 3 N unit matrix. All the eigenvectors φI corresponding to each eigenvalue λI ðI 5 1; 2; 3; . . .; NÞ can be represented in an orthogonal matrix:  Φ 5 ϕ1 ϕ2 ? ϕN ; ΦT 5 Φ21 ; (11.85c)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

439

which transform matrix A into a diagonal matrix in the form ΦT AΦ 5 diagðλI Þ 5 Λ;

A 5 ΦΛΦT ;

AΦ 5 ΦΛ:

(11.85d)

n

Based on this property, the matrix A is represented by An 5 ΦΛn ΦT :

(11.85e)

The spectral radius of the matrix A is defined as the largest modulus of its eigenvalues, that is, ρðAÞ 5 Max jλI j;

(11.85f)

I

so that if the spectral radius of the amplification matrix is not larger than 1, that is, ρðAÞ 5 Max jλI j # 1; I

(11.85g)

the scheme is stable, since the matrix Λn is finite with n increases, and the matrix An is finite as well. Furthermore, if we require the spectral radius to be less than 1, that is, ρðAÞ 5 Max jλI j , 1; I

(11.85h)

we will have lim Λn 5 0;

n-N

lim An 5 0:

n-N

(11.85i)

Therefore any calculation errors generated in the integration process will tend to zero as the time step increases.

Spectral analysis of numerical errors Referring to Eqs. (11.76a) and (11.76b), we can represent the variable unI in the solution vector Un in a Fourier series: unI 5

N X J52N

N X

vnJ ejIκJ Δx 5

vnJ ejIφJ ;

(11.86a)

J52N

where vnJ and φJ 5 κJ Δx are the amplitude and the phase angle of the Jth harmonic component of the variable unI . The amplification factor is now given by the ratio of the amplitude, omitting subscript J, that is, A5

vn11 5 Aðφ; Δt; ΔxÞ: vn

(11.86b)

The amplification factor A allows an evaluation of the frequency distribution of the discretization errors generated by the numerical scheme, which defines the numerical values of the time evolution of the solution; the amplitude vn of the harmonic component with wavenumber κ can be written in the form 2jωtn 2jωnΔt ^ ^ 5 vðκÞe ; vn 5 vðκÞe

n 5 0; 1; 2; . . .;

(11.86c)

^ where ω 5 ωðκÞ is a complex frequency corresponding to the real wavenumber κ. The amplitude vðκÞ can be determined by the Fourier decomposition of the initial condition of the problem, such as in Eq. (11.73) for the convection problem: ð 1 L ^ 5 vðκÞ f ðxÞe2jκx dx: (11.86d) 2L 2L Therefore the harmonic components κ in Eq. (11.86a), when Eq. (11.86c) is introduced, becomes  n 2jωnΔt jIκΔx ^ e ; uI κ 5 vðκÞe

(11.86e)

and the amplification factor can be denoted as A 5 e2jωΔt :

(11.86f)

440

FluidSolid Interaction Dynamics

There is an error of this amplification factor in the numerical scheme compared with the exact amplification function of the differential equation, such in the same type of function of the exact variable denoted by over waves (B): ~ : A~ 5 e2jωΔt

(11.86g)

This allows us to investigate the frequency property of the numerical errors. The complex frequency ω can be written as ω 5 ξ 1 jη;

(11.86h)

from which, when substituted into Eq. (11.86f), it follows A 5 eηΔt e2jξΔt 5 jAje2jϕ ;

ϕ 5 ξΔt;

jAj 5 eηΔt :

In a similar way to dealing with the exact amplification function in Eq. (11.86g), we obtain





~ ~

A~ 5 eη~ Δt : A~ 5 eη~ Δt e2jξΔt 5 A~ e2jϕ~ ; ϕ~ 5 ξΔt;

(11.86i)

(11.86j)

The error in amplitude, called the diffusion or dissipation error, is defined by the ratio of the computed amplitude to the exact one: εD 5

j Aj ; eη~ Δt

(11.86k)

The error in the phase, the dispersion error, is defined as the difference εϕ 5 ϕ 2 ϕ~

(11.86l)

ϕ εϕ 5 : ϕ~

(11.86m)

for parabolic problems, and, for convection problems,

For example, for the hyperbolic problem in Eq. (11.73), the exact solution is a single wave propagating with the velocity a, so that ϕ~ 5 κaΔt: Before ending this subsection, we should mention the following points: 1. This subsection discussed the stability problem by investigating several simple numerical methods in a 1D case. The basic ideas can be extended to 2D or 3D cases. The main aim of this book is to discuss FSI, not detailed CFD methods, so this content is omitted. Interested readers may refer to the many published books on CFD methods to get these details, such as Hirsch (1988) and Chung (2002), or to the more original publications cited in this book. 2. The von Neumann method is valid for linear problems. For nonlinear problems, it can provide the stability condition at a local point, around which the nonlinear problem is linearized by using the method given in Section 11.2.3.5. It has been proven that von Neumann analysis can provide a necessary condition for the stability of nonlinear problems (Ritchtmyer and Morton, 1967). 3. In Appendix A, based on the von Neumann method, we investigate the spectral radius and the corresponding stability and accuracy of five time integration methods used in numerical analysis. This discussion provides more examples for readers to understand how to use the von Neumann method for stability analysis of numerical schemes. 4. For nonlinear problems, an energy flow method for stability analysis of generalized nonlinear partial differential equations has been developed in Xing (2015). This method, based on the equations of a nonlinear system represented in the phase space, defines a generalized energy flow variable, a scalar, through which the stability and other nonlinear behavior of the system are investigated.

11.2.3.5 Nature of flow equations with its quasi-linear form We consider a general system of n first-order partial differential equations for the n unknown functions u j in an mdimensional space xk (k 5 1,2,. . .,m), represented in the conservation form occupying the summation convention of repeated superscripts or subscripts: @Fik 5 Ri ; @xk

k 5 1; 2; . . .; m; i 5 1; 2; . . .; n:

(11.87a)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

441

By introducing a Jacobian matrix Jk, of which its element is Jijk 5

@Fik ; @uj

(11.87b)

we obtain the quasi-linear form of Eq. (11.87a) as follows: Jijk

@uj 5 Ri ; @xk

k 5 1; 2; . . .; m; i; j 5 1; 2; . . .; n:

(11.87c)

This equation can be rewritten in a matrix form Jk

@U 5 R; @xk

k 5 1; 2; . . .; m;

(11.87d)

where U and R are the unknown variables vector and the nonhomogeneous source vector, respectively, while Jk are (n 3 n) matrices. If the homogeneous equation of Eq. (11.87d) vanishes, that is,  k ^ 5 0; (11.88a) J νk U there will exist a wave solution ^ jðνUxÞ 5 Ue ^ jðν k x Þ ; U 5 Ue k

(11.88b)

^ propagating in the normal direction of surface S with its front surface S(x) 5 Ct, on which the wave amplitude is U, (x) 5 Ct, ν k 5

@S=@xk

: @S=@xk

(11.88c)

This wavefront surface is called the characteristic surface of the wave, as shown in Fig. 11.8. Eq. (11.88c) defines the normal for the wavefront surface, which corresponds to the vanishing determinant of the system that is,

k

J ν k 5 0; (11.88c) and which can have, at most, n solutions with n characteristic surfaces. According to the nature of the n characteristic normals, Eq. (11.87a) has the following mathematical nature: G

G G G

Hyperbolic: if all the n characteristic normals are real with the n linear independent characteristic solutions of Eq. (11.88b) Elliptic: if all the characteristic normals are complex Hybrid: if some are real and others complex Parabolic: if there are fewer than n real characteristic normals Here we give two examples to illustrate these mathematical characteristics, as follows.

FIGURE 11.8 Wavefront surface and its unit normal vector.

442

FluidSolid Interaction Dynamics

Two first-order equations in two-dimensional space Considering a system in the 2D space, @u @v 1 c 5 q1 ; @x @y @v @u b 1 d 5 q1 ; @x @y

a

or in its matrix form



a 0

0 0 @ u 1 d b @x v

Therefore with x1 5 x; x2 5 y; we have



a J 5 0 1

0 ; b

(11.89a)

q c @ u 5 1 : q2 0 @y v

0 J 5 d 2

leading to the determinant in Eq. (11.88c) becoming





a 0

aν =ν 0 c

1 ν 2

5

1 2 ν

0 b 1 d d 0

c ; 0

c

5 0; bν 1 =ν 2

(11.89b)

(11.89c)

(11.89d)

from which the characteristic normals follow:  2 ν1 cd 5 : ab ν2

(11.89e)

If cd/ab . 0, the system is hyperbolic, for example, c 5 d 5 a 5 b 5 1 with a vanishing right-hand side, we obtain the well-known wave equation @2 u @2 u 2 5 0: @x2 @y2

(11.89f)

If cd/ab . 0, the system is elliptic, for instance, a 5 b 5 1 and c 5 2 d 5 2 1; leading to the Laplace equation @2 u @2 u 1 5 0; @x2 @y2

(11.89g)

which is the standard form of elliptic equations describing diffusion phenomena. If b 5 0, there is only one characteristic normal v2 5 0, and the system is parabolic. For example, with a 5 1, b 5 0, c 5 2 d 5 2 1; and q1 5 0; q2 5 v; we obtain the standard parabolic equation @u @2 u 5 : @x @y2

(11.89h)

Stationary shallow-water equations These equations describe the spatial distribution of the height h of the free surface in a stream with velocity components u and v, in the gravitational field of acceleration g, which is in the following matrix form: 2 3 2 3 2 3 2 3 u h 0 h v 0 h h @ 4 g u 0 5 4 u 5 1 4 0 v 0 5 @ 4 u 5 5 0; (11.90a) @x @y 0 0 u v g 0 v v of which the three characteristic normals λ 5 ν x =ν y can be obtained as the solutions of equation:

2 3



uλ 1 v hλ h



4 gλ 5 5 0: uλ 1 v 0



g 0 uλ 1 v

(11.90b)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

From this equation, it follows

    ðuλ 1 vÞ u2 2 gh λ2 1 2uvλ 1 v2 2 gh 5 0;

and the three solutions λ

ð1Þ

5 2 v=u;

λ

ð2;3Þ

5

2 uv 6

443

(11.90c)

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ghðu2 1 v2 2 ghÞ : u2 2 gh

(11.90d)

pffiffiffiffiffi It can be seen that gh is a critical velocity and that the system is hyperbolic for the velocities of u2 1 v2 . gh. Otherwise, the system is hybrid, since the first solution λð1Þ is always real, of which the characteristic surface is the ð1Þ ð1Þ streamline, since the vector νð1Þ has components proportional to ν ð1Þ x 5 2 v and ν y 5 u; therefore ν Uv 5 0; implying that the normal of the characteristic surface is perpendicular to the flow velocity.

11.2.4 Updated arbitrary LagrangianEulerian formulations in finite difference method and finite volume method 11.2.4.1 Generalized numerical conservative laws in finite volume method Using the same matrix approach in Eqs. (11.67a)(11.67d) and the general conservative law in Eq. (3.86) for the updated ALE system, we can express the conservative laws for that system in the following matrix forms. Matrix integration form The general integral form of the conversation laws in the updated ALE coordinate system is given by Eq. (3.86a), from which follows its matrix form in the same structure as in Eq. (11.67c), that is, ð ð ð _ @ _ UdΩf 1 ðQ i 1 Gi Þν i dΩf 5 BdΩf ; (11.91a) @t Ωf Ωf Ωf where

2 6 Gi 5 4

0 2 σ~ ij

3 7 5;

2vj σ~ ij 1 hi 3 ρf _ 6 7 ρf vj U 54 5; ρf ðe 1 vj vj =2Þ 2

2 _ 6 Qi 5 4

ρf ðvi 2 v~i Þ ρf vj ðvi 2 v~i Þ 1 pδij

ρf ðe 1 vj vj =2Þðvi 2 v~i Þ 1 pvi 2 3 0 6 7 ρf fj B54 5:

3 7 5; (11.91b)

ρf ðfj vj 1 θÞ _

It should be noticed that the definition of vector matrix Q i is now different compared with Qi in Eq. (11.67d) for the Eulerian domain fixed in the space. For the updated ALE system, the domain is a moving one with the mesh mov_ ing velocity v~i involved, so that the matrix Q i is a function of the mesh velocity v~i . Also, if we consider that the mechanical energy conservation is automatically satisfied for momentum conservation, as given by Eq. (3.105a), the third lines in the matrices in Eq. (11.92b) can be reduced by eliminating the related mechanical variables and leaving only the thermo-variables. Discretized matrix form in finite volume method In a similar approach to deriving Eq. (11.68), we have the discretized formulation of the general conservation laws in order to be implemented in the numerical analysis based on the FVM, that is, ð ð ð @ _ _ ðIÞ UdΩðIÞ 1 ð Q 1 G Þν dΩ 5 BdΩðIÞ (11.91c) i i i f f f : ðIÞ ðIÞ @t ΩðIÞ Ωf Ωf f Numerical process For us to illustrate the numerical process more clearly, we discuss a 2D system as developed by Hirt et al. (1974, 1997), in which the fluid is assumed nonviscous, so that the viscous stress σ~ ij  0 in Eq. (11.91b). A typical mesh

444

FluidSolid Interaction Dynamics

5 3 B

2 A

4 6

1

D

C

9

8 7

FIGURE 11.9 Two-dimensional mesh structure, where the dashed lines enclose the integration volume for the momentum.

structure is shown in Fig. 11.9 in which the dashed lines enclose the momentum integration volume. The velocity variable vi is assigned at each vertex, and for each cell, the mass density ρ, the cell mass M, the internal energy e, and the cell volume Ω are assigned. These variables are advanced with time governed by Eqs. (11.91a)(11.91c). The numerical processes are described as follows. Input data: Mesh vertex coordinates: xðIÞ I 5 1; 2; . . .; 8; i 5 1; 2; i ; Mesh vertex velocities: vðIÞ ; I 5 1; 2; . . .; 8; i 5 1; 2; i Cell densities: ρα ; α 5 A; B; C; D; Cell internal energies: eα ; α 5 A; B; C; D: Initial calculations Cell volume Ωα : Such as ΩA 5 Aera of Δ123 1 Aera Δ134; Cell mass:

Vertex mass:

Cell total energy:

Mα 5 ρα Ωα ;

Such as M ð4Þ 5

(11.92b)

D 1X Mα ; 4 α5A

Such as EA 5 eA 1

(11.92a)

4 1X vðIÞ vðIÞ ; 8 I51 i i

Cell pressure is calculated by equation of state:

pα 5 f ðρα ; eα Þ:

(11.92c)

(11.92d)

(11.92e)

Here, the vertex mass is obtained by assuming that the mass in each cell and the cell kinetic energy in Eq. (11.92c) are equally shared among its four corner vertices. Lagrangian calculation: phase 1, explicit In this step, we assume that the fluid undergoes an imagined Lagrangian motion in order to advance the variables by the time progressing a step Δt, so that we obtain the values of the variables, which will be modified in the next phase by _ considering the mesh motions. In this step, the matrix Q i in Eq. (11.91b) to be used in the calculations is now given by 2 3 0 _ (11.93a) Q i 5 4 pδij 5: pvi Since we do not consider the viscous stress and thermo-effect, the matrix Gi vanishes. Therefore the advance of _ variables U is contributed to by the pressure gradients and the body force B. On considering these conditions, the

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

445

momentum law in Eq. (11.91c) at each vertex, such as vertex 4, is approximated in the following numerical approximate form: I @vð4Þ j M ð4Þ 5 2 pð4Þ ν j dΓ 1 M ð4Þ fj ; (11.94a) @t 1368 from which, when considering the outside normal vector of the volume of vertex shown in Fig. 11.9, _ð4Þ

ν 1 5 dx2 =dΓ ; ν 2 5 2 dx1 =dΓ , and approximating the time derivative, it follows that the temporary value v j of the velocity is calculated by I Δt _ð4Þ v 1 5 vð4Þ 1 2 pð4Þ dx2 1 Δtf1 1 M ð4Þ 1368  (11.94b)        Δt ð1Þ ð3Þ ð3Þ ð6Þ ð6Þ ð8Þ ð8Þ ð1Þ 5 vð4Þ 1 x 2 x x 2 x x 2 x x 2 x ; 1 p 1 p 1 p =2 1 Δtf p A B C D 1 1 2 2 2 2 2 2 2 2 M ð4Þ I Δt _ð4Þ 1 pð4Þ dx1 1 Δtf2 v 2 5 vð4Þ 2 M ð4Þ 1368  (11.94c)        Δt ð3Þ ð1Þ ð6Þ ð3Þ ð8Þ ð6Þ ð1Þ ð8Þ 5 vð4Þ 1 x 2 x x 2 x x 2 x x 2 x ; 1 p 1 p 1 p =2 1 Δtf p A B C D 1 2 1 1 1 1 1 1 1 1 M ð4Þ where a factor 1/2 is introduced for the middle integration along the boundary line of the volume of vertex 4, since the value of the pressure in a cell is assigned as a constant shared by the two half-cell areas, such as pA in cell A shared by Δ123 and Δ134. For viscous fluids, there exists the viscous stress σ~ ij which will do work to advance the cell energy, so that in this step, the work W ðIÞ for each vertex can be calculated using a similar integration as given by Eqs. (11.94a)(11.94d) based on the σ~ ij -related integration in the energy conservation law Eq. (11.91c). The work done by the pressure is not calculated in this step but will be calculated in the second phase. The advanced energy for each cell, such as for cell A, is adjusted to E A 5 EA 1

4 1X W ðIÞ : 4 I51

(11.94d)

If there is no viscous stress, as in our discussion here, this work calculation is not necessary. However, for problems involving shock waves, it is necessary to add to p an artificial viscous pressure q 5 2 λρrUv that is proportional to the velocity divergence. For the expanding cells, this artificial pressure is replaced by zero (Hirt et al., 1974, 1997). Therefore the work done by this artificial pressure in shock wave problems needs to be calculated in this step. Lagrangian calculation: phase 2, implicit First part: When an explicit calculation is sought, this step can be omitted. The aim of phase 2 is to obtain the new velocities that have been accelerated by the time-advanced pressure gradients. The time-advanced pressures are determined by the function of mass density and energy, as shown in Eq. (11.92e). When vertices are moved with their new velocities, which in turn are functions of the new pressures, these pressures are defined implicitly and must be determined by iteration. For example, the desired pressure marked by * of the cell A will be   (11.95a) pA 2 f ρA ; eA 5 0; for which the initial values of variables at iteration cycle n are approximated as   n pnA ΩA  n ΩA  n ρA 5 ρA  ; eA 5 eA 1 n 1 2 n : ΩA ρA ΩA

(11.95b)

Here, physically, the first equation implies the total mass ( 5 ρA ΩA ) conservation of cell A, and the increment ( 5 eA 2 enA ) of the energy equals the work done by the time-advanced pressure per unit mass. ΩA is the volume of the cell when its vertices are moved to xðIÞ 5 xðIÞn 1 vi Δt: i i

(11.95c)

446

FluidSolid Interaction Dynamics

This volume is computed in terms of coordinates shifted with velocities accelerated by the new pressures through formulae Eqs. (11.94a)(11.94d), in which the pressure p of each cell is replaced by p. The NewtonRaphson iteration to Eq. (11.95a) is considered an implicit equation for pA through Eqs. (11.94a)(11.94d), (11.95a), and (11.95b). _ The velocity v i values obtained by Eqs. (11.94a)(11.94d) in the previous step are used as initial guesses for the iteration. For each sweep, the following adjustments to each cell must be done: _

1. Compute Ω using the most updated values for velocities v i : 2. Compute new guesses for ρA and eA from Eq. (11.95b). 3. Make the pressure change calculation by means of   f ρA ; eA 2 pA ; ΔpA 5 γA

(11.95d)

where γ A is a relaxation factor, which is the rate

   Δ pA 2 f ρA ; eA ; γA 5 ΔpA

(11.95e)

to be determined numerically. 4. Adjust the current pressure pA by adding ΔpA , and adjust the velocities at the corners of the cell to reflect this pressure change. For example, for Eqs. (11.94a)(11.94d), due to ΔpA , the new velocities are   Δt ð3Þ ΔpA xð1Þ ; 2 2 x2 ð4Þ 2M   Δt _ð4Þ _ð4Þn ð3Þ ð1Þ Δp x 2 x v2 5 v2 1 : A 1 1 2M ð4Þ _ð4Þ v1

_ð4Þn

5 v1

1

(11.95f)

5. The mesh is repeatedly swept, and the preceding calculations are performed until no cell exhibits a pressure change violating the condition



Δp



p , ε; max

(11.95g)

where pmax is the actual or an estimated maximum pressure in the mesh and ε is a suitably chosen small number of order 1023 . For a given ε, the necessary iteration number to obtain convergence can be roughly estimated as follows. Based on the NavierStokes Eq. (3.133), the velocity change and the pressure change satisfy the approximate order as Δv Δp B ; Δt Δx

Δv

Δx BΔp; Δt

ðΔvÞ2 BΔp;

(11.95h)

which, when divided by the standard sound velocity and the corresponding pressure, gives the nondimensional relation



Δp

(11.95i) M 2 B



5 ε; p where M is the Mach number. Therefore the effective sound speed ceff is given by c2eff 5

v2 ; ε

(11.95j)

where v denotes a typical fluid speed. To satisfy the stability, the Courant number in Eq. (11.81b) must satisfied, that is,   ceff Δt # 1; (11.95k) Δx N where Δx is a typical cell dimension, and N is the necessary iteration number. The ratio Δt=N is the time for each iteration calculation process. Therefore the necessary iteration number will be of the order

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

N

vΔt ε0:5 Δx

:

447

(11.95l)

This result implies that the iteration number increases as either Δt increases or ε decreases but is independent of the actual material sound speed. Thus, once a tolerable level of pressure error ε has been chosen, the implicit scheme converges in a finite number of iterations, regardless of the actual material sound speed. Second part: The final values of v and p from the iteration in the first part are now the new Lagrangian values for the cycle. To complete the Lagrangian portion of a cycle, the cell energies must now be adjusted for the pressure work terms omitted in phase 1, and vertices must be moved with the fluid to their new positions. The energy for cell A in Fig. 11.9 is changed by adding the work done by the pressures I Δt  p v ν  dΓ ; (11.96a) EA 5 E A 1 4MA 1234 i i where the integration is defined along the boundary lines of cell A. The cell-edge pressures are obtained by a mass weighting scheme. For example, the pressure on edge 34, the line connecting nodes 3 and 4 in Fig. 11.9, is given by p34 5

MB pA 1 MA pB ; MB 1 MA

(11.96b)

which was recommended by Fromm (1961). The mesh vertices are moved with the fluid to their new locations 



xð4Þ 5 xð4Þn 1 Δtvð4Þ i i i :

(11.96c)

Here the new velocities are used to move the vertices since this makes the explicit Lagrangian portion of a cycle second-order accurate in time. The new densities of the cells can be obtained by their mass divided by the new cell volume, calculated based on the new location by Eq. (11.96c) of each cell, such as ρA 5

MA : ΩA

(11.96d)

The phase 1 and phase 2 calculations comprise an implicit Lagrangian method that is stable for any Courant number cΔt=Δx with the fluid speed of sound c. Phase 3: Rezoning calculations based on the arbitrary LagrangianEulerian system Lagrangian cell methods do not adequately describe flows undergoing large distortions in nonlinear FSI problems. Therefore we have to consider the mesh motions in an ALE system to eliminate the large distortions. For this purpose, the values of variables obtained in the Lagrangian calculations must be adjusted by the mesh velocity v~i . The adjustments associated with a shift in the position of a typical vertex, say 4 in Fig. 11.9, proceed as follows: 

xðIÞn11 5 xðIÞ 1 v~ðIÞ i i i Δt;

(11.97a)

where I 5 4 for vertex 4. With vertex 4 moving, the lines connecting it to its neighbors 1, 3, 6, and 8 sweep out volumes containing mass and total energy that must be exchanged among the adjacent cells. For example, if vertex 4 moves to the right, the grid line connecting 4 to 3 sweeps out of cell A and adds to cell B a volume equal to the summation of the areas of two small triangles in the form   i Δth ð4Þ  ð3Þ ð3Þ ΔΩ 5 xð4Þ v~1 x2 2 xð4Þ 1 v~ð4Þ : (11.97b) 2 2 1 2 x1 3 The associated mass change with this volume change can be calculated by an average weighting formulation, such as the mass subtracted from cell A and added to cell B is 1 MA 1 MB ΔM 5 ðΔΩ 1 αjΔΩjÞ 1 ðΔΩ 2 αjΔΩjÞ ; 2 ΩA 2 ΩB

(11.97c)

where α is called a donor cell weighting factor, with α 5 0, the flux is centered, and with α 5 1, the flux is the full donor cell. As discussed by Hirt et al. (1974), stability requires that this factor satisfies 0 , α , ΔΩ=Ω , 1 and stays as small as possible. Here, Ω is an average cell volume.

448

FluidSolid Interaction Dynamics

The total energy subtracted from cell A and added to cell B is 1 MA EA 1 MB EB ΔðMEÞ 5 ðΔΩ 1 αjΔΩjÞ 1 ðΔΩ 2 αjΔΩjÞ : 2 2 ΩA ΩB

(11.97d)

The shift of vertex 4 is also accompanied by a momentum exchange among vertices 1, 3, 6, and 8, since 4 is a corner of the control volumes for these vertices. For example, the line connecting vertices 4 and 2 sweeps out a volume   i Δth ð4Þ  ð2Þ ð4Þ ð4Þ ð2Þ ~ v~1 x2 2 xð4Þ x 2 x ΔΩ 5 1 v : (11.97e) 2 2 1 1 3 The mass in this volume is MA ΔΩ=ΩA , and the momentum it contains in the x1 direction is approximated as ΔðMv1 Þ 5

MA ð1Þ ðΔΩ 2 αjΔΩjÞvð3Þ 1 1 ðΔΩ 1 αjΔΩjÞv1 ; 2ΩA

(11.97f)

which must be subtracted from vertex 1 and added to vertex 3. The momentum change in the x2 direction is handled in the same way. When these exchanges among the cells and vertices surrounding the moved vertex have been completed for the cells A, B, C, and D, the mesh is ready for any other vertex to move to a new location. After all vertices have been moved to the positions desired, the total cell masses and energies are converted back to densities and specific energies, and vertex momenta are converted back to velocities. New specific internal energies e can be calculated by subtracting the average cell kinetic energy based on Eq. (11.92d). New cell pressure is computed from the state equation in Eq. (11.95a) in terms of the new values of ρ and e.

11.2.4.2 Generalized numerical conservative laws in finite difference method Differential matrix equation for finite difference method Applying the Green theorem to transfer the surface integration in Eq. (11.91a) into the volume integration, and considering that the integration volume is arbitrary, we obtain the differential matrix form of the conservation laws for finite difference method (FDM), that is, _

_ @U

@ðQ i 1 Gi Þ 1 5 B: (11.98a) @t x @xi _

In this equation, the variable matrices are given in Eq. (11.91b), in which the matrix Q i contains the mesh velocity v~i . For convenience in the following discussion on ALE formulations, this matrix is rewritten in the form 2 3 0 _ _ _ _ (11.98b) Q i 5 ðvi 2 v~i ÞU 1 pi ; pi 5 4 pδij 5; pvi so that Eq. (11.98a) becomes

 _ _ _

@ ðvi 2 v~ÞiU @U

@pi @Gi 1 1 1 5 B: @t x @xi @xi @xi

(11.98c)

_

Here, the vector U is the conservative variable, the vector B denotes the source vector, and the flux vector contains _ _ two components: the inviscid flux involving the vectors U and p i as well as the viscous term contributed by the viscous stress Gi. Numerical matrix equation for finite difference method In the updated ALE description, the current spatial coordinates are taken as the mesh coordinates; therefore the partial time derivative with x fixed in Eq. (11.98c) essentially is the partial time derivative with the mesh coordinate z fixed. To solve this equation using FDM, the discretization in time and space, as discussed in Section 11.2.2, are required. Actually, many well studied schemes for solving this equation are available in published form and are also effectively used in commercial CFD computer codes. For example, Hirsch (1990) and Chung (2002) give more details. The schemes include explicit and implicit ones and later require iteration calculations, similar to those previously discussed for FVM.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

449

Here, we do not wish to discuss these detailed solution approaches, since this book concerns the solutions of FSI problems, for which we seek only to confirm that any effective CFD schemes for fluids can be chosen together with FEA methods for solids to deal with nonlinear FSI problems based the proposed mixed FEFDM. As an illustration for FDM solution, we consider the fluid with no thermo-variables involved, so that the general momentum equation given by Eq. (11.98c) is   

@ ρðvi 2 v~i Þvj @ pδij @ðρvj Þ

@σ~ ij 5 ρf 2 2 1 : (11.99a) j @t x @xi @xi @xi In this equation, the subindex i obeys the summation convention discussed in Chapter 2, Cartesian tensor and matrix calculus. For this equation, the straightforward central discretization of second-order space accuracy in FD form leads to the semidiscretized system of ordinary differential equations in time. For example, in the case of a 1D incompressible fluid (ρ 5 1), Eq. (11.99a) reduces to

~ @v

@½ðv 2 vÞv @p @σ~ @v 2 1 ; σ~ 5 υ ; (11.99b) 5f 2

@t x @x @x @x @x for which the following first-order explicit central schemes

~ I11 2 ½ðv2 vÞv ~ I21 pI11 2 pI21 σ~ I11=2 2 σ~ I21=2 fI11 2 fI21 @vI

½ðv2 vÞv 2 1 1 ; 52 @t x 2Δx 2Δx Δx 2Δx

(11.99c)

where the viscous shear stress is defined by σ~ I11=2 5

υI11=2 ðvI11 2 vI Þ : Δx

The explicit scheme is obtained from a first-order accurate, forward time difference as follows:     Δt ½ðv2 v~ÞvtI11 2 ½ðv2 v~ÞvtI21 Δt ptI11 2 ptI21 t1Δt t vI 2 2 vI 5 2 2Δx 2Δx  t   t1Δt  t t t1Δt υΔt vI11 2 2vI 1 vI21 Δt fI11 2 fI21 ~ ; xt1Δt 1 1 5 xtI 1 vΔt: I Δx2 2Δx

(11.99d)

(11.99e)

~ so that the Courant number can be In the ALE numerical process, the fluid velocity relative to the mesh is v 2 v, ~ σ 5 ðv 2 vÞΔt=Δx, and the stability condition can be written as σ2 # 2β # 1;

β5

υΔt : Δx2

(11.99f)

For practical problems with complex geometry, the coordinate transformation as discussed in Section 3.1.3 is needed to transform the complex physical domain into the standard computational domain. For the updated ALE system, the mesh coordinate ξj in Eq. (3.19) is just the current special coordinate xi for the time period (t; t 1 Δt).

11.3 Mixed finite elementcomputational fluid dynamics solutions for nonlinear fluidsolid interaction problems In Sections 11.1 and 11.2, the FE formulations of solids in the UL system and the FVM/FDM formulations of fluids in the updated ALE system are given. Due to the interaction on the fluidstructure interface, the governing equations for solids and fluids cannot be solved separately. For solving the equations governing the solid dynamics, the fluid forces on the interaction interfaces must be given, while to solve the fluid equations, the solid motion on the interaction interface has to be known. The governing equations for solids and fluids have to be solved simultaneously. There are two approaches to solving FSI equations as follows.

11.3.1

Direct or simultaneous integration

The direct or simultaneous solution procedure sets up the integrated FSI equations and solves them simultaneously (Xing et al., 2002, 2003, 2005). This procedure is more suitable for solving linear FSI problems using available

450

FluidSolid Interaction Dynamics

numerical integration schemes with no interactive iteration process. As discussed in Chapter 7, Finite element models for linear fluidstructure interaction problems, for linear FSI analysis, the nonsymmetrical coupling equations for a mixed FE model can be symmetrized (Xing, 1984; Xing et al., 1991, 1996) in order to directly use the available computer codes for symmetrical FE equations.

11.3.2

Partitioned iteration

Another solution process is the partitioned iteration procedure (Xing et al., 2002, 2003, 2005), which is more suitable to integrating nonlinear FSI equations. In this method, the dynamics of the fluid, structure, and even electric unit, if electricmechanical coupling is involved, are solved separately, and the data are exchanged at every time step or iteration. As the flowchart shows in Fig. 11.10, the main steps of the iteration for the waterstructure interactions are as follows: 1. Initial calculations: Input initial conditions and physical parameters and calculate the initial values of all variables at the initial equilibrium state at time t 5 0. These are the conditions to be used in order to start calculations for the solution at time 0 1 Δt. 2. Fluid iteration calculations: Solve the fluid equation at time t 1 Δt based on the known solution of the FSI system at time t to derive the approximate fluid velocity and stress or pressure. 3. Solid iteration calculations: Use the approximate fluid stress or pressure describing the load on the FSI interface to calculate the corresponding solid displacement, velocity, and acceleration. 4. Repeat steps 2 and 3 until a suitably chosen convergence condition is reached; then go into the next time period from t 1 Δt to t 1 2Δt: Start Read geometry mesh. Physical parameters. Initial velocity, stress, and wave height Yes

t=t+Δt No

Update mesh

Solid displacement |Uk –Uk–1|≤ ε

Compute geometry Set boundary conditions Solve velocity and stress or pressure Calculate wave height

Calculate fluid forces

Solve displacement, velocity, and acceleration of structure

Velocity divergence ≤ ε

No

Yes No

Iteration finished ? Yes

Stop

FIGURE 11.10 Flowchart of partitioned solution procedure (Xing et al., 2002, 2003, 2005).

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

11.3.3

451

A simple example

To more directly understand the details of the mixed FECFD method for nonlinear FSI problems, we discuss the simple example shown in Fig. 11.2, which now includes the information on the water domain, as shown in Fig. 11.11. In this system, the fluid domain is considered a 2D water pond of depth 2h and width 2B. The elastic rod is inserted along the symmetric middle line of the pond. To distinguish the notations for coordinates, we replace the spatial coordinates (x, y) by ðx1 ; x2 Þ in the equations describing the fluid motions. We assume that the water is an incompressible and nonviscous fluid of mass density ρf 5 1 with a free surface where a nonlinear wave disturbance ζðx2 ; tÞ due to the gravitational acceleration g is considered. The two vertical walls and the base of the water pond are rigid, so that the water velocities perpendicular to the walls and the base vanish, and therefore we approximately consider the vertical motion of the water only, with horizontal motion neglected. Also, considering the symmetry of the system in its geometry and excitation forces, we discuss the right-hand half system in order to reduce the degrees of freedom in the numerical equations. The equations governing the water velocity vi ðxj ; tÞ and dynamic pressure pðxj ; tÞ in its motion are as follows.

11.3.3.1 Governing equations of water domain Mass conservation vi;i 5 0;

ðxj ÞAΩf :

(11.100a)

Momentum conservation @vi @vi  @p 1 vj 2 v~j 1 5 gδ1i ; @t @xj @xi

ðxi ÞAΩf :

(11.100b)

Boundary conditions On the two vertical walls and the rigid pond base, the normal velocity of water vanishes, that is, v2 5 0; v1 5 0;

x2 5 6 B; x1 5 L 1 h:

(11.100c)

On the free surface, a moving boundary for nonlinear cases, its geometrical position, and the corresponding outside normal are changed with time and at different space points, while the force on it always equals the atmospheric pressure assumed to be 0. Therefore we have to satisfy the following boundary conditions. Free surface is described by the equation

FIGURE 11.11 Nonlinear rodwater interaction system; the positive wave height is in the x direction.

452

FluidSolid Interaction Dynamics

ζ ðxi ; tÞ 5 x1 2 ζ ðx2 ; tÞ 5 L 2 h;

ζ ðx2 ; 0Þ 5 0;

1 ν 1 5 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 1 ζ 2;2

ζ ;2 ν 2 5 qffiffiffiffiffiffiffiffiffiffiffiffiffiffi; 1 1 ζ 2;2

(11.100d)

on which the dynamic and kinematic conditions are required: Dynamic condition: p 5 0; Kinematic condition:

05

Dζ @ζ @ζ 5 1 ðvi 2 v~i Þ : Dt @t @xi

(11.100e) (11.100f)

At the end point of the rod, that is, point t x2 , a moving point of the rod, the velocities of the rod and the water must be consistent, and the water pressure must equal the negative node force t T 2 of the rod at node 2, as shown in Eq. (11.39a). On the wet rod surface (line x2 5 0), the velocity in the y direction vanishes, also due to symmetricity. There are no horizontal flows in the y direction on the symmetric central line; therefore we have Kinematic condition: Dynamic condition:

v1 5 t u_2 ; x1 5 t x2 ;   p t x2 ; 0; t 5 2 t T 2 ;

Wet rod surface and symmetric line:

v2 ðx1 ; 0; tÞ 5 0:

(11.100g) (11.100h) (11.100i)

In this set of equations, the water velocity vi ðxj ; tÞ, the pressure pðxj ; tÞ, and the wave height ζðx2 ; tÞ are unknown variables to be determined.

11.3.3.2 Differential matrix formulations Artificial compressibility Eqs. (11.100a) and (11.100b) represent the conservation equations of mass and momentum, respectively, between which there exists no coupling due to the incompressible water. The fluid velocity obtained by solving Eq. (11.100b) has to satisfy Eq. (11.100a), which is considered a constraint condition on Eq. (11.100b). To overcome this restriction, the artificial compressibility relation (see, e.g., Hirsch, 1995) is introduced: @p 5 2 βvi;i ; @t

(11.101)

where β is an artificial compressibility constant, and t a pseudo-time variable. A relatively large value of Δt is chosen to approximate an incompressible fluid. Based on Eq. (11.101), the equations to be solved in terms of FDM are represented in the matrix form 2 3 2 3 2 3 2 3 p βv 0 @ 0 @4 5 @ 4 i5  4 vi 5 1 4 gδ1i 5: vi 5 2 p 2 vj 2 v~j (11.102) @t @xi @xj 0 ζ ζ 0 This equation with the boundary conditions given in Eqs. (11.100c)(11.100i) can be solved with the FDM schemes discussed in this chapter. To start the numerical calculations, the following initial conditions of the system are required: Pressure: Velocity: Wave height: Right-half fluid domain: FSI interface:

pðxi ; 0Þ 5 gðx1 2 L 1 hÞ; vi ðxj ; 0Þ 5 0; ζðx2 ; 0Þ 5 0; Ωf ðxi ; 0Þ; L 2 h , x1 , L 1 h; 0 , x2 , B; 0 x1 5 0 x2 5 L:

(11.103)

Here, 0 x2 5 L is the initial position of the end of the rod. The water pressure at the initial time equals the static water pressure.

11.3.3.3 Integrated dynamic equations From Eqs. (11.91a) and (11.91b), the integrated dynamic equations for this simple problem are given by ð ð ð  @ vi 2 v~i 0 1 vj dΩf 1 ν i dΓ f 5 dΩf ; @t Ωf Γ f ðvi 2 v~i Þvj 1 pδij Ωf gδ1j

(11.104)

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

453

in which the first line gives the mass conservation of the incompressible fluids, and the second line is the momentum conservation of the system. The energy conservation equation is not involved because no thermos-effect is considered in this simple example. Based on this equation, we can establish an approximate numerical equation in the next subsection.

11.3.3.4 Numerical equations Based on Eq. (11.104) and the method given in Section 11.2.4.1, the numerical equations can be established. However, to derive a hand-solved equation for this simple problem, we divide the right-half water domain into two subdomains with six vertices, as shown in Fig. 11.12, in which the details at the end of rod shown in Fig. 11.13 are neglected. Generally, we can use the interpolation method as used in FEA (Bathe, 1996) to represent the current position coordinates and the velocity and pressure distribution in the water by the interpolation functions and the vertex values, such as for the mesh in Fig. 11.11. These are given in the following form: xi 5

6 X

H ðIÞ ðr; sÞxðIÞ i ;

I51

vi 5 p5

6 X

H ðIÞ ðr; sÞvðIÞ i ;

(11.105)

I51 6 X

H ðIÞ ðr; sÞpðIÞ ;

I51

where r and s are local coordinates with the values at each vertex shown in Fig. 11.12. The interpolation functions H ðIÞ are given by H ð1Þ 5 0:25ð1 2 rÞð1 2 sÞ 2 0:5ð1 2 r 2 Þð1 2 sÞ; H ð2Þ 5 0:25ð1 2 rÞð1 1 sÞ 2 0:5ð1 2 r 2 Þð1 1 sÞ; H ð3Þ 5 0:5ð1 2 r 2 Þð1 2 sÞ; H ð4Þ 5 0:5ð1 2 r 2 Þð1 1 sÞ;  xðIÞ 1 2L ; H ð5Þ 5 0:25ð1 1 rÞð1 2 sÞ 2 0:5ð1 2 r 2 Þð1 2 sÞ; r 5 h   xðIÞ 2 2L H ð6Þ 5 0:25ð1 1 rÞð1 1 sÞ 2 0:5ð1 2 r 2 Þð1 1 sÞ; s 5 : B

(11.106)

FIGURE 11.12 Two finite volumes of the right-half water domain in Fig. 11.11, in which the details at the end of the rod shown in Fig. 11.13 are neglected.

454

FluidSolid Interaction Dynamics

1 2

4 3

A

5

6

FIGURE 11.13 Details of mesh vertices at the end of the rod for the boundary integrations.

However, for this simple example, we do not use these generalized interpolations to define the variable distributions in the water domain. Rather, considering the boundary conditions and the symmetry of the problem, we choose the following vertex values of the velocity and the pressure in the water. Then we select one with which to obtain the linear distribution of a variable on a boundary line between the two vertices. These approximations are described as follows. The wave height and velocity on free surface as a material surface t

ζ ðx2 ; t Þ 5 t ζ 1 1

t ð1Þ x1

 x2 t ζ2 2 tζ1 ; B t ð2Þ x2

5 L 2 h 1 ζ 1; t

ζ ð1Þ 5 t ζ 1 ;

5 L 2 h 1 ζ 2; t

ζ ð2Þ 5 t ζ 2 ;

t _ ð1Þ x1

5

t_

ζ 1 5 vð1Þ 1 ;

(11.107a) t _ ð2Þ x1

5

t_

ζ 2 5 vð2Þ 1

The mesh velocity The free surface and the wet interface are moving boundaries of the system, but the base and two walls are fixed in the space. Therefore, to follow the fluid motions on the free surface and the vertical motion of the wet interface, we consider the free surface and the wet interface with the vertical motion of line 3A4 as the material surfaces but the base and two walls as the Eulerian boundaries, so that the mesh velocity is given by vi 5 v~i ; on free surface; ð3Þ ðAÞ ðAÞ ðAÞ t vð3Þ xð3Þ 1 5 v~1 5 v1 5 v~1 ; 1 5 x2 5 x1 ; vð4Þ 1

5 v~ð4Þ 1 ;

v~ð5Þ 1

5 0 5 v~ð6Þ 1 ;

v~ðIÞ 2

5 0;

(11.107b)

ðI 5 1; 2; 3; A; 4; 5Þ:

From Eqs. (11.107a) and (11.107b), we can obtain the following vertex velocities. The fluid velocitizes at vertices ð1Þ

t_ t_ vð1Þ 1 5 x1 5 ζ1; ð2Þ vð4Þ 1 5 v1 =2;

vðIÞ 2 5 0;

ð2Þ

t_ t_ vð2Þ 1 5 x1 5 ζ 2;

vð5Þ 1 5 0;

t _2 5 vðAÞ vð3Þ 1 5 u 1 ;

vð6Þ 1 5 0;

ðI 5 1; 2; 3; 4; 5; 6Þ;

(11.107c)

ð1Þ vðAÞ 2 5 γv1

Here, the horizontal velocity on the symmetrical line vanishes due to the symmetry consideration, but its value along ð1Þ line 3A at the end of the rod is not zero, so vðAÞ 2 5 γv1 is assumed. The parameter γ . 0 is chosen, since for a positive ð1Þ ð2Þ v1 , the rod moves down and pushes water to move to the right. The velocity vð4Þ 1 5 v1 =2 is obtained by the linear interð2Þ ð6Þ polation based on the velocities v1 at vertex 2 and v1 5 0 at vertex 6. The water pressure Generally, the interpolation formulation in Eq. (11.105) can provide the distribution of the pressure, but here we assume that the pressure consists of its static and a dynamic pressure proportional to the wave height in the form   (11.107d) p 5 g x1 2 L 1 h 2 t ζ ; in t Ω: Based on the motions of fluid and the mesh structure given by Eqs. (11.107a)(11.107d), we can use the integrated form of the conservation laws in Eq. (11.104) to generate the approximate numerical equations as follows.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

455

Mass conservation Referring to Fig. 11.12 with the details on the boundary in Fig. 11.13 and considering the symmetry of the system, the first integration in Eq. (11.104) over the half volume on the right-hand side is calculated by ð I @t Ω 5 2 ðvi 2 v~i Þν i dΓ f 5 2 ðvi 2 v~i Þν i dΓ f : (11.108a) @t Γf 56421A3 Based on the following line integrations along the boundary lines, ð6 ð6 v1 ν 1 dΓ f 5 0 3 ð2 1Þdx2 5 0; 5 5 ð2 ð2 v2 ν 2 dΓ f 5 0 3 ð1Þð 2dx1 Þ 5 0; 6 6 ð1 ðA 0ν i dΓ f 5 0; 0ν 2 dΓ f 5 0; 2 1 ð3 ð3 ðvi 2 v~i Þν i dΓ f 5 0ν 1 dΓ f 5 0: A

(11.108b)

A

In the calculations of these line integrations, we have introduced the boundary conditions given in Eqs. (11.107a) (11.107d), as well as the zero components of the unit outside normal vectors on the integrated boundaries. Using the assumed wave height in Eq. (11.107a), the volume of the water equals the summation of the area of a trapezoid and the area of a rectangle under the rod end, as shown in Fig. 11.13, that is,           2h 2 t ζ 1 1 2h 2 t ζ 2 h 2 t u2 ð1 1 4BÞh 2 B t ζ 1 1 t ζ 2 2 t u2 t 3 B 2 1=2 1 Ω5  ; (11.108c) 2 2 2 where the length of line 3A is 1/2, that is, the half width of the rod. When the results in Eq. (11.108a) and (11.108b) are substituted into Eq. (11.108a), we obtain   2 B t ζ_ 1 1 t ζ_ 2 2 t u_2 @t Ω 2 Bðv1 1 v2 Þ 2 t u_2 5 5 5 0: (11.108d) @t 2 2 Here, B 2 1=2  B has been introduced. Momentum conservation The second line in Eq. (11.104) gives the integrated momentum conservation law: ð ð ð  @ vj dΩf 1 ðvi 2 v~i Þvj 1 pδij ν i dΓ f 5 gδ1j dΩf ; @t Ωf Γf Ωf

(11.109a)

from which, when using the velocity and pressure on each boundary line based on the vertex values given in Eqs. (11.107a)(11.107d), as well as completing the related boundary integrations and the volume integration in terms of a lump-mass (lump-volume) approximation considering the mass volume being average-shared by its vertices of each subdomain, the numerical momentum conservation equations can be derived. The volume integration is given by " " ð1Þ # # ð ð3Þ ð4Þ t_ v1 v1 1 vð2Þ @ Bh ζ 1 1 2t ζ_ 2 1 2t u_2 1 1 2v1 1 2v1 0 @ 0 @ dΩf 5 Ω 5 Ω ; 0Ω 5 ; (11.109b) t ðAÞ @t Ωf v2 @t @t 3 2γ ζ_ 1 2v2 and the boundary integrations for velocity take the form ð I v1 v ðvi 2 v~i Þ ðvi 2 v~i Þ 1 ν i dΓ f 5 0: ν i dΓ f 5 v v2 2 Γf 56421A35 The boundary integration of the pressure is given by

(11.109c)

456

FluidSolid Interaction Dynamics

ð Γf

ν1 ν2





I pdΓ f 5 5621A35

ð6 ð2 ð5 1 0 0 pdΓ f 5 pdx2 2 pdx1 2 pdx1 ; 0 1 ν2 5 6 1 1 ν1

(11.109d)

where the integration on line A3 is neglected since the half rod width is much less than B. Based on the pressure in Eq. (11.107d) and the wave height in Eq. (11.107a), we obtain the following integrations:   ð6 ð6 ð6 t   1 1  1 x2 t ζ2 1 tζ1 1 t t t ζ 2 ζ 1 dx2 5 gB 2h 2 pdx2 5 g L 1 h 2 L 1 h 2 ζ dx2 5 g 2h 2 ζ 1 2 ; B 2 2 0 5 0 5 0 5 0 (11.109e) t   ð2  ð2    L2h2 ζ 2 3t ζ 22 0 0 0 0 t t t pdx1 5 g x1 2 L 1 h 2 ζ 2 dx1 5 g x1 =22L1h2 ζ 2 x1 5 2 2gh h 2 ζ 2 2 ; 1 1 4h 6 1 6 1 L1h (11.109f)  L1h  ð5 ðA    t 2  0 0 0 3 ζ1 0 pdx1 5 g x1 2 L 1 h 2 t ζ 1 dx1 5 g x1 =22L1h2 t ζ 1 x1 5 2gh h 2 t ζ 1 2 : 4h 1 1 1 1 1 L2h2t ζ 1 1 (11.109g) Substituting these results into Eq. (11.109d), we obtain ! ! ð t ν1 1 0 ζ2 1 tζ1 3t ζ 22 t pdΓ f 5 gB 2h 2 1 2gh h 2 ζ 2 2 2 4h 0 1 Γ f ν2 3 2 t ζ2 1 tζ1 1 2 ! 7 6 4h 7 6 0 3t ζ 21 t 6 t 2 t 2  7 2 2gh h 2 ζ 1 2 5 2gBh6 t 7: t 4h 1 4 ζ1 2 ζ2 3 ζ 1 2 ζ 2 5 1 4hB B The body force integration is

ð Ωf

ð gδ1j dΩf 5

g Ωf

1 g dΩf 5 60 Ω : 0 0

(11.109h)

(11.109i)

Substituting the results into Eq. (11.109a) and noting 0 Ω 5 Bh=3 in Eq. (11.109b), we obtain the numerical momentum conservation equation as 3 2 t ζ2 1 tζ1 7 " # 6 4h t t 7 6 t @ ζ_ 1 1 2 ζ_ 2 1 2 u_2 7 6 : (11.110) 5 6g   6 t 2 t 2 7 7 6 tζ 2 tζ @t 3 ζ 2 ζ 2γ t ζ_ 1 2 1 5 1 4 2 1 4hB B Combining Eqs. (11.110) and (108d), we have the numerical equations for the fluid flow, in which there are three equations and four unknown variables: two wave heights t ζ 1 and t ζ 2 , one parameter γ to define the velocity t ðAÞ v 2 5 γ t ζ_ 1 5 t vA , as well as the rod velocity t u_2 at the rod end. For the solution of the problem, the solid equation in Eq. (11.48c) must be used.

11.3.3.5 Mixed finite elementcomputational fluid dynamics equations Based on the FE Eq. (11.48c) for the rod under the UL description and Eq. (11.110) with Eq. (11.108d) for the fluid in the updated ALE system, we obtain the mixed FECFD numerical equations for this simple nonlinear FSI problem. In order to use the partitioned iteration approach proposed in Section 11.3.2, the FSI equations can be modified by

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

457

eliminating t ζ_ 2 in the left-hand term in Eq. (11.110) in the base of Eq. (11.108d) and moving the term of the rod motion t u_2 on the right-hand side, so that we can rewrite the finalized FSI equations in the following form: 2 3 t 1 Δt ζ 2 1 t 1 Δt ζ 1 2 3 " # 2 2 6 7 2h 6 7 t 1 Δt 22 @ t 1 Δt ζ_ 1  71 B 5; t 1 Δt ζ_ 2 5 2 t 1 Δt u_2 =B 2 t 1 Δt ζ_ 1 ; u€2 4 5 3g6 6 t 1 Δt ζ 2 t 1 Δt ζ @t t 1 Δt vA 3 t 1 Δt ζ 22 2 t 1 Δt ζ 21 7 4 5 2 1 0 1 4hB B (11.111) L ρS t 1 Δt L ρS ω u^ iωðt1ΔtÞ e ; u€2 1 t K u2 5 2 t 1 Δt pðtÞ 2 t σ 11 1 6 3

0 0

0 0

2

where the pressure t 1 Δt pðtÞ can be calculated by Eq. (11.107d) at point ð t 1 Δt x1 5 L 1 t 1 Δt u2 ; height t 1 Δt ζ 5 t 1 Δt ζ 1 at point 1, that is,   t 1 Δt p 5 g h 1 t 1 Δt u2 2 t 1 Δt ζ 1 :

(11.112) x2 5 0Þ with wave (11.113)

In this example, since we have used the assumed distributions of the fluid velocity and the pressure in the water domain, the resultant numerical equations do not involve the special FD scheme. Using time integration approaches, we can numerically solve this simple nonlinear FSI problem. We may give further explanations on the approximate process to solve this simple example. The assumed functions for the velocity and the pressure of the water are linear in the full fluid domain, which allows us to derive the handsolvable numerical equations. Therefore readers, especially students, can clearly understand the mathematical process in completing the integration calculations in the integrated conservation laws based on the updated ALE system, which may emphasize the following: G

G

G

G

G

G

The free surface is defined as a material boundary with its moving velocity v~i equaling the fluid velocity vi , in order to trace the free surface motions, so that v~i 2 vi  0 on integrations over the free surface. On the wet interface, that is, the end of the rod, to guarantee no separation between the rod and the water, the mesh vertical velocity v~1ð3AÞ 5 u_2 is defined, while its horizontal velocity v~2ð3AÞ 5 0 is chosen to allow different horizontal water flows v2ð3AÞ along the line 3A. In calculations of the boundary integrations, since the outside normal vectors of the boundaries are always along the horizontal or vertical directions, its components have the value 0, 1, or 21, from which the boundary integrations are easily completed. The boundary integrations are defined on the boundary, so that generally the differential element dΓ f is a positive area and, in this example, a positive line element. Therefore take care about how to express this positive line element by dxi . For example, in the line integration from points 6 to 2, dΓ f 5 2 dx1 , and the normal vector ν 1 5 0; ν 2 5 1: However, in the line integration from points 1 to 5, dΓ f 5 dx1 , and the normal vector ν 1 5 0; ν 2 5 2 1: The main aim of this “by hand” example is to explain how to use the conservation laws to establish the numerical equations; therefore we do not wish to solve it for the numerical results, which would be too approximate. Although the errors of this simple approximation may be large, it still could be used for initial approximation estimations in the first stage of practical designs in engineering. For practical complex FSI problems, computer calculations are needed to obtain the solutions, and we will discuss the mixed FEFDM for nonlinear rigid bodywater interactions with some application examples in Section 11.4 as follows.

11.4

Numerical examples by the mixed finite elementfinite difference method

In previous subsections, we presented the fundamental theory and methods to deal with nonlinear FSI problems by using the mixed FECFD method. The simple examples have been given for hand-mathematical work to illustrate how to derive the numerical equations based on the described theory. In this subsection, we will discuss some of practical numerical examples of nonlinear structurewater interactions based on the FEFDM to show the applications of the

458

FluidSolid Interaction Dynamics

proposed numerical approaches. In these examples, for simplicity, the structures are treated as rigid bodies undergoing large motions, so that the ideas and concepts are clarified before considering a fully flexible solid. The following examples were reported in previous publications by Price and Xing (2000) and Xing et al. (2002, 2003, 2005).

11.4.1

Description of nonlinear rigid bodywater fluidsolid interaction systems

11.4.1.1 Reference frames Fig. 11.14 illustrates an arbitrarily shaped 3D rigid floating body traveling on the calm surface of a fluid. In this figure, o 2 y1 y2 y3 represents a Cartesian coordinate-axis system fixed in a 3D space; O 2 x1 x2 x3 denotes a moving equilibrium coordinate system with its origin at the mass center of the rigid body, located at a point Oðy01 ; y02 ; y03 Þ at time t 5 0, and this remains parallel to o 2 y1 y2 y3 at all times. The axis system O 2 X1 X2 X3 is fixed in the rigid body at O such that it coincides with O 2 x1 x2 x3 in the absence of any disturbance at time t 5 0. When disturbed, the mass center O is allowed to translate with displacement u0i , velocity v0i , and acceleration w0i , and the body is allowed to rotate with angular displacement θi around the mass center O. Therefore the moving coordinate system O 2 x1 x2 x3 is a noninertial system with, in general, acceleration w0i 6¼ 0. The coordinate system O 2 X1 X2 X3 is the material frame of reference describing the motion of the rigid body. Between the moving frame O 2 x1 x2 x3 and the fixed coordinate system o 2 y1 y2 y3 , there exist the following relations: yi 5 y0i 1 xi ;

y_i 5 v0i 1 vi ;

y€i 5 w0i 1 wi ;

(11.114)

where v0i and w0i denote the velocity and acceleration at the origin of the moving system, while y_i , vi , and y€i wi represent the velocities and accelerations of the particles relative to the fixed and moving frames, respectively.

11.4.1.2 Governing equations The governing equations describing this dynamic systems are presented in Cartesian tensor form as follows. Fluid domain It is assumed that the fluid is incompressible. Based on the equations presented in Section 3.6.2, in the moving coordinate-axis system O 2 x1 x2 x3 by introducing the inertial force, 2ρf w0i , the equations describing fluid flows are as follows. Conservation equation of mass vj;j 5 0;

ðxi ; tÞAΩf 3 ðt1 ; t2 Þ

x2

x3 p0

Ti

Ωs O

Γv vi

ηi

ST

vi

(11.115a)

Fi

ζ(y1, y2, t)

Γf

x1

Σ

Γp p∞

fi

Ωf Γb y3 y2

Reference plane y3 = 0

o FIGURE 11.14 Rigid floating body traveling on the surface of a fluid.

y1

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

459

Conservation equation of momentum ρf vi;t 1 ρf ðvi vj Þ;j 5 ρf f^i 2 ρf w0i 1 σij;j ;

ðxi ; tÞAΩf 3 ðt1 ; t2 Þ

(11.115b)

Constitutive equation σij 5 2 pδij 1 2μVij ;

ðxi ; tÞAΩf 3 ðt1 ; t2 Þ

(11.115c)

Strain-ratevelocity relation Vij 5

 1 vi;j 1 vj;i ; 2

ðxi ; tÞAΩf 3 ðt1 ; t2 Þ

(11.115d)

Nondimensional forms Assuming a characteristic body length L and velocity v0 giving a characteristic time L/v0 and pressure ρf v20 , we find that the nondimensional forms of the governing equations in Eqs. (11.115a)(11.115d) are expressed as 0 1 2 v;t 5 2 @φδij 1vi vj 2 Vij A ; Re ;j (11.116a) ðx ; tÞAΩ 3 ðt ; t Þ; vj;j 5 0; or

ðvi;j 1 vj;i Þ ; Vij 5 2

i

f

1

2

  1 v;t 5 2 φδij 1vi vj 2 vi;j ; Re ;j

(11.116b)

where Re denotes the Reynolds number which with Fr the Froude number are as follows ρ f v0 L ; μ

v0 Fr 5 pffiffiffiffiffiffi; gL x δ l l3 φ 5 p 1 w0i xi 1 2 : Fr Re 5

(11.116c)

Boundary conditions On the free surface Γ f , the unknown wave height disturbance function ζ 5 y03 1 ζ~ relative to the fixed axis o 2 y1 y2 y3 describes the surface elevation measured from a chosen base level, such as the horizontal reference plane y3 5 0 in Fig. 11.14. Here, ζ~ represents the height of the free surface relative to the origin O of the moving system. This disturbance function satisfies the relation ~ 1 ; x2 ; tÞ 5 0; y3 2 ζðy1 ; y2 ; tÞ 5 x3 2 ζðx

(11.117a)

with the material derivative Dðy3 2 ζÞ=Dt 5 0 because the free surface is a material surface. This produces the kinematic condition of the free surface: v3 5 ζ~ ;t 1 v1 ζ~ ;1 1 v2 ζ~ ;2 ;

ðxi ; tÞAΓ f 3 ðt1 ; t2 Þ:

(11.117b)

The traction force on the free surface is the prescribed atmospheric pressure p^0 , which is assumed to be perpendicular to the free surface if the viscosity of the air and the surface tension properties of the fluid are neglected. The dynamic condition on the free surface is expressible in the form σij ν j 5 2 p^0 ν i ;

ðxi ; tÞAΓ f 3 ½t1 ; t2 :

(11.117c)

On the inlet boundary Γ v in Fig. 11.13, the fluid velocity is prescribed, and the boundary condition is vj 5 v^j 2 v0j ;

ðxi ; tÞAΓ v 3 ½t1 ; t2 :

(11.118)

On the outlet boundary Γ p , the outflow is assumed to be approximately unidirectional, and the surface pressure is of the known value p^N 5 p^0 to be in the negative direction of vi. Therefore, on this boundary, the stress continuity condition requires that

460

FluidSolid Interaction Dynamics

σij ν j 5 2 p^0 ν i ;

ðxi ; tÞAΓ p 3 ½t1 ; t2 :

(11.119)

On the fixed wall boundary Γ b , the velocity of viscous fluid is zero, and the appropriate boundary condition is vj 1 v0j 5 0;

ðxi ; tÞAΓ b 3 ½t1 ; t2 :

(11.120a)

For an ideal fluid flow, only the normal velocity on the fixed boundary vanishes, but the tangential velocity is not prescribed and is therefore free, that is,   vj 1 v0j ν j 5 0: (11.120b) Solid domain Rotation tensor and its time derivatives Based on the rotation tensor Rij and its time derivatives given by Eqs. (10.156) and (10.157), we can obtain the current position xi, relative displacement uRi , velocity vRi , and acceleration wRi at time t of a particle Xi in the rigid body under the moving coordinate system: xi 5 Rij Xj ; uRi vRi wRi

Xi 5 Rji xj ;     5 xi 2 Xi 5 Rij 2 δij Xj 5 δij 2 Rji xj ; 5 x_i 5 R_ij Xj 5 eilk θ_ l Rkj Xj 5 eilk θ_ l xk ; 5 R€ij Xj 5 eilk θ€ l Rkj Xj 1 eilk ekrs θ_ l θ_ r Rsj Xj 5 eilk θ€ l xk 1 eilk ekrs θ_ l θ_ r xs

(11.121a)

5 eilk θ€ l Rkj Xj 2 θ_ r θ_ r Rij Xj 1 θ_ l θ_ i Rlj Xj 5 eilk θ€ l xk 2 θ_ r θ_ r xi 1 θ_ l θ_ i xl : In the acceleration formulation of a particle, the first, second, and third terms on the right-hand side denote the tangential, centripetal, and Coriolis accelerations, respectively. In the material coordinate system, the displacement Ui, velocity Vi, acceleration Wi of the particle given in Eq. (11.121a) are obtained by an inverse rotation as follows:   Ui 5 Rji uRj 5 Rji Rjr 2 δjr Xr 5 ðδir 2 Rri ÞXr ; Vi 5 Rji x_j 5 Rji R_jr Xr 5 Rji ejlk Rkr θ_ l Xr ; (11.121b) Wi 5 Rji R€jr Xr 5 Rji ejlk θ€ l Rkr Xr 1 Rji ejlk ekms θ_ l θ_ m Rsr Xr 5 Rji ejlk θ€ l Rkr Xr 2 Rji θ_ m θ_ m Rjr Xr 1 Rji θ_ l θ_ j Rlr Xr : The absolute motion of each particle Xi in the rigid body can be represented by a summation of the translation u0i of the mass center O and its relative motion, that is, ui 5 u0i 1 uRi ;

vi 5 v0i 1 vRi ;

wi 5 w0i 1 wRi :

(11.121c)

The equations describing the motion of the rigid body can be derived by considering the equilibriums of all forces (including the inertial one) and their moments about the mass center at the current position in the system O 2 x1 x2 x3 . When we establish these two equations, we have to remember that the external forces applied in and on the solid body are defined in the material coordinate system O 2 X1 X2 X3 , while the motion variables in Eq. (11.121c) are defined in the system O 2 x1 x2 x3 . Therefore a suitable transformation is needed in order to keep all involving variables in the same coordinate system. Equation of translation motion of the mass center ð ð ð mu€0i 5 ρS Rij F^ j dΩ 1 Rij T^ j dS 1 Rij Tj dS; (11.122a) ΩS

ST

Σ

where m denotes the total mass of the body, and the rotation tensor Rij is used to transform the solid forces defined in the system O 2 X1 X2 X3 into the system O 2 x1 x2 x3 in which the force balance equation is derived. The result of the inertial force of the rotation acceleration wRi in Eq. (11.121a) vanishes, since its expression includes the coordinate Xr , so that the volume integration over the body of mass center at the origin O vanishes.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

461

Equation of rotation motion about the mass center The resultant of inertial force caused by the translation acceleration acts at the origin O, which produces any moment about the mass center, and only the rotation acceleration wRi makes a contribution for the moment balance equation. Considering this equation in the system O 2 x1 x2 x3 , we obtain ð ð ð ð R ^ ^ ρS eijk xj wk dΩ 5 ρS eijk xj Rkl F l dΩ 1 eijk xj Rkl T l dS 1 eijk xj Rkl Tl dS; (11.122b) ΩS

ΩS

Σ

ST

from which, when the acceleration in Eq. (11.121a) is substituted, it follows that ð ð ð I il θ€ l 1 Iiln θ_ l θ_ n 5 ρS eijk xj Rkl F^ l dΩ 1 eijk xj Rkl T^ l dS 1 eijk xj Rkl Tl dS; ΩS

where

ð I il 5

Σ

ST

(11.122c)

ð ΩS

ρS eijk eklm xj xm dΩ 5

ΩS

ρS eijk eklm Rjr Rms Xr Xs dΩ

5 eijk eklm Rjr Rms I~rs 5 δil I~rr 2 Ris Rlr I~rs ; ð ~I rs 5 ρS Xr Xs dΩ; I~rs 5 I~sr ; I il 5 I li ; ΩS ð   Iiln 5 ρS eijk eklt etnm xm xj dΩ 5 δil ejnm 2 δjl einm Rmr Rjs I~rs ;

(11.122d)

ΩS

Iiln 5 Iiln 1 I^iln ; ðIiln 2 Iinl Þ ; I^iln 5 2

ðIiln 1 Iinl Þ ; Iiln 5 2

Iiln 5 Iinl ;

I^iln 5 2 I^inl ;

and the relationship based on the e 2 δ identity in Eq. (2.14)   eijk eklt etnm 5 δil δjt 2 δit δjl etnm 5 δil ejnm 2 δjl einm

(11.122e)

and the identity J iln θ_ l θ_ n 5 0

(11.122f)

have been used in deriving the resultant equation. The I~sr represents the second-order inertial moment matrix of the body about its mass center in the material coordinate system, which is independent of the rotation. If the material coordinate system is chosen as the principal inertial system of the body, then I~sr 5 0; r 6¼ s. For small rotation cases, we do not need to distinguish the coordinates xi and Xi, so that Rij in Eq. (11.122d) can be replaced by δij to obtain its reduced form. Fluidstructure interaction interface Let us assume that there exists no discontinuity on the fluidstructure interaction interface Σ during motions. This implies that both the velocity vi of the viscous fluid and the velocity Vi of the solid are the same at each point, that is, v i 5 Vi ;

ðxi ; tÞAΣ 3 ½t1 ; t2 :

(11.122g)

The dynamic equilibrium condition on the interaction interface is represented as σijf ν j 1 σijS ηj 5 σijf ν j 1 Ti 5 0;

ðxi ; tÞAΣ 3 ½t1 ; t2 :

(11.122h)

If the fluid is assumed to be nonviscous, condition Eq. (11.122g) reduces to v i ν i 5 Vi ν i ;

ðxi ; tÞAΣ 3 ½t1 ; t2 

(11.122i)

with condition Eq. (11.122h) changing to p 1 ηi σSij ηj 5 p 1 ηi Ti 5 0; τ i σSij ηj 5 τ i Ti 5 0:

ðxi ; tÞAΣ 3 ½t1 ; t2 ;

(11.122j)

462

FluidSolid Interaction Dynamics

Here, superindices f and S are used to identify the stress in the fluid and the solid, respectively, and τ i in Eq. (11.122j) denotes the unit tangential vector of the wet interface. Initial conditions At time t1 5 0, the initial configuration and velocity distributions of the system must be given. Therefore the initial conditions are assumed as follows: u0i 5 u^0i 5 y^0i ;

11.4.2

v0i 5 v^0i ;

0 θi 5 θ^ i ;

0 θ_ i 5 θ_^ i ;

vi 5 v^i0 2 v^0i ;

ðxi ÞAΩf ;

ζ 5 ζ 0;

ζ~ 5 ζ~ 0 5 ζ 0 2 y^03 :

(11.123)

Arbitrary LagrangianEulerian description of the water motion

An ALE mesh coordinate system ξ i is introduced to transform the physical domain Ωf under the moving reference frame O 2 x1 x2 x3 into the computational or reference domain Ωz 5 Ωξ . This transformation, as discussed in Section 3.1.3, is represented as ξi 5 ξ i ðxj ; tÞ;

τ 5 t;

(11.124a)

xi 5 xi ðξ j ; τÞ;

t 5 τ;

(11.124b)

with the assumed inverse transformation

and Jacobian determinant and the transformation relation for the corresponding volume elements are, respectively, given by



@xi

~ z 5 JdΩ ~ ξ: (11.124c) J~ 5

j

. 0; dΩf 5 JdΩ @ξ The mesh velocity is given by v~i 5

@xi ; @τ

(11.124d)

which describes the mesh relative velocity in the moving frame O 2 x1 x2 x3 . Its divergence links with the time rate of the Jacobian determinant, given by Truesdell (1991), in the form @J~ ~ ~ v: 5 J v~i;i 5 Jdiv~ @τ

(11.124e)

The application of Eqs. (3.19d)(3.19j) transforms the nondimensional governing equations in Eqs. (11.116a) and (11.116b) to the forms defined in the mesh coordinate system ðξj ; τÞ. That is, from Eq. (11.116a), we have the mass conservation equation   ~ i @ξ j =@xi @ Jv 5 0; (11.125) ~ j J@ξ and, from Eq. (11.116b), the momentum conservation equation 0 1   j ~ i ~ @vi @ξk @ξj @ Jv J @ @ ~ @ξj @ξ A: ~ i ðvm 2 v~m Þ 52 Jφ 1 Jv 2 ~ ~ j @xi @xm Re @ξk @xm @xm J@τ J@ξ

(11.126)

Similarly, the kinematic condition in Eq. (11.117b) on the free surface transforms to @ζ~ @ζ~ @ζ~ 5 v3 2 ðv1 2 v~1 Þ 2 ðv2 2 v~2 Þ @τ @x1 @x2     ~ j =@x1 ~ j =@x2 @ J~ζ@ξ @ J~ζ@ξ 5 v3 2 ðv1 2 v~1 Þ 2 ðv2 2 v~2 Þ : ~ j ~ j J@ξ J@ξ

(11.127)

The rest of the boundary conditions for the fluid, as well as the boundary conditions on the fluidstructure interaction interface, remains unchanged.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

11.4.3

463

Numerical formulations

11.4.3.1 Artificial compressibility Eqs. (11.125) and (11.126) represent the conservation equations of mass and momentum, respectively. For an incompressible fluid, there is no coupling between these equations. The fluid velocity obtained by solving Eq. (11.126) has to satisfy Eq. (11.125), which is considered a constraint condition on Eq. (11.126). To overcome this restriction, as assumed in Eq. (11.101), the artificial compressibility relation (see, e.g., Hirsch, 1995) is introduced:   ~ m @ξj =@xm @ Jv @φ 52β ; (11.128) ~ j @~τ J@ξ where β is an artificial compressibility constant, and τ~ a pseudo-time variable. A relatively large value of Δτ~ is chosen to approximate an incompressible fluid.

11.4.3.2 Discretization of momentum equation To help numerical manipulations, Eq. (11.126) is represented in matrix form:   @ e j 2 eνj @d 52 ; @τ @ξj where 2

3 v1 6 7 d 5 J~4 v2 5;

2

v1 ðv1 2 v~1 Þ 1 Sj1 φ

(11.129a) 3

6 7 j 7 ej 5 6 4 v2 ðv2 2 v~2 Þ 1 S2 φ 5;

v3 2 j k 3 Sm Sm @v1 =@ξk 1 6 j k 7 ejν 5 4 S Sm @v2 =@ξk 5; ~ JRe m Sjm Skm @v3 =@ξk

(11.129b)

v3 ðv3 2 v~3 Þ 1 Sj3 φ ~ j =@xm Sjm 5 J@ξ vj 5 vm @ξj =@xm : v~j 5 v~m @ξj =@xm

(11.129c)

In these equations, vj and v~j are the contravariant components of velocity v and v~ under the mesh coordinate system ðξ ; τÞ, but vi and v~i represent the components of the velocity under the Cartesian coordinate system (xi, t). In this study, only laminar flow is assumed, and a constant viscosity value is chosen to derive viscous fluxes. j

Time derivatives The application of a backward-difference implicit formula of second-order accuracy (Hirsch, 1995) to the time derivative of Eq. (11.126) gives 1:5dn11 2 2dn 1 0:5dn21 5 2 rn11 ; Δτ

(11.130)

where r represents a residual vector, and superscripts n 1 1, n, and n 2 1 denote the physical time points of the calculation. The pseudo-time derivative in Eq. (11.128) is represented by an implicit Euler back-forward difference formula, which gives   ~ j Þ n11;m11 φn11;m11 2 φn11;m β @ðJv 5 2 n11 : (11.131) Δτ~ @ξj J~ This allows Eq. (11.130) to be rewritten as 1:5dn11;m11 2 1:5dn21;m 1:5dn11;m 2 2dn 1 0:5dn21 5 2 rn11;m11 2 ; Δτ Δτ with the second superscript m denoting quantities at the mth pseudo-time iteration level to be introduced. The combination of Eqs. (11.131) and (11.132) gives the matrix numerical iteration equation

(11.132)

464

FluidSolid Interaction Dynamics

2 4Iτ J~

n11



@R 1 @D

n11;m

3

  ~ n11;m 2 2D 5ΔDn11;m11 5 2 Rn11:m 2 Im 1:5D ~ n 1 0:5D ~ n21 ; Δτ

(11.133a)

where  T  T  T D 5 φ v 1 v 2 v 3 5 φ vT ; v 5 v 1 v 2 v 3 ;  T  T ~ 5 JD; ~ Ejν 5 0 ejT ; v~ 5 v~1 v~2 v~3 ; D ν  j  j  @ E 2 Eν ~ j ejT T ; R5 ; Ej 5 Jβv j @ξ   1 1:5 1:5 1:5 ; ; ; Iτ 5 diag ; Im 5 diagð0; 1; 1; 1Þ; Δτ~ Δτ Δτ Δτ

(11.133b)

(11.133c)

and the linearization relation 

R

n11;m11

5R

n11;m

 n11;m n11;m  n11;m11  @R @R n11;m n11;m 5R 1 D 2D 1 ΔDn11;m11 @D @D

(11.133d)

is used.

Spatial derivatives A well-established implicit numerical scheme of study developed by Kwak et al. (1986) and by Rogers and Kwak (1991) has been used for the numerical solution of Eq. (11.133a). The derivatives of the convective fluxes Ej in vectors R and @R=@D are discretized using a first-order upwind difference scheme based on the flux-difference splitting method (see, e.g., Shen et al., 1996; Chen et al., 1998, 1999). These terms are therefore expressed in the forms  j  j @E @Eν Rα;β;γ 5 1 ; (11.134a) j @ξ α;β;γ @ξ j α;β;γ  j Ejα1δ1L ;β1δ2L ;γ1δ3L 2 Ejα;β;γ @E 5 ; (11.134b) @ξj α;β;γ Δξj  j Ejνðα1δ1L ;β1δ2L ;γ1δ3L Þ 2 Ejνðα2δ1L ;β2δ2L ;γ2δ3L Þ @Eν 5 ; (11.134c) @ξj α;β;γ 2Δξj where index L 5 j 5 1, 2, 3. The numerical flux Ejα;β;γ is calculated by adopting the approximate Riemann solver by Roe (1981) as follows: 8 j   9 E Dα;β;γ 1 Ej Dα2δ1L ;β2δ2L ;γ2δ3L > > = <



1

@Ej  L  L ; (11.134d) Ejα;β;γ 5



D 2

D > 2> ; : @D α;β;γ α;β;γ DLα;β;γ 5 Dα;β;γ 2 Dα2δ1L ;β2δ2L ;γ2δ3L ;   L Dα;β;γ 5 0:5 Dα;β;γ 1 Dα2δ1L ;β2δ2L ;γ2δ3L :

(11.134e)

From Eqs. (11.134a)(11.134e), it is observed that Rα;β;γ is represented by the values of D at the 7 points ðα; β; γÞ, ðα 6 1; β; γÞ, ðα; β 6 1; γÞ, and ðα; β; γ 6 1Þ. Similarly, the derivative ð@R=@DÞα;β;γ can be represented by variable values at these 7 points. The final difference representation of Eq. (11.133a) is given by 2 3  n11;m   ~ n11;m 22D 4Iτ J~n11 1 @R 5ΔDn11;m11 5 2 Rn11;m 2 Im 1:5D ~ n 10:5D ~ n21 : (11.134f) α;β;γ α;β;γ α;β;γ α;β;γ @D α;β;γ Δτ

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

465

11.4.3.3 Discretization of kinematic condition on free surface Time derivatives The wave elevation of the free surface is updated at each iterative step by means of the unsteady kinematic boundary condition expressed in Eq. (11.127). The discretization of the time derivative in this equation is the same as the one used in the momentum Eq. (11.126), and therefore Eq. (11.127) is written in the discretization form: 1:5ζ~

n11;m n11;m n n21 2 1:5ζ~ 1:5ζ~ 2 2ζ~ 1 0:5ζ~ n11;m11 52γ ; 2 Δτ Δτ

n11;m11

(11.135)

where the residual term γ n11;m11 is calculated by an approximation similar to Eq. (11.133d), that is, γ n11;m11 5 γ n11;m   ~ j =@x1 @ J~ζ@ξ 1 ~ j J@ξ

  ~ j =@x2 @ J~ζ@ξ ~ j J@ξ

3n11;m11 !n11;m 2 v1 2 v~1 21 Δ4 v2 2 v~2 5 v3

(11.136a)

5 γ n11;m 1 An11;m Δðv1vÞn11;m11 ; where 3n11;m @ξ1 @ξ1 07 " #n11;m 6 7 6 @x1 @x2 ~ @ζ~ @ ζ 7 6 n11;m 2 2 7 6 21 A 5 ; 7 6 @ξ @ξ @ξ 1 @ξ2 07 6 5 4 @x1 @x2 0 0 1 2 3n11;m11 v1 2 v~1 Δðv1vÞn11;m11 5 Δ4 v2 2 v~2 5 : v3 2

(11.136b)

(11.136c)

Spatial derivatives The spatial derivatives of Eq. (11.126) are computed by a second-order central difference scheme of the form ! ζ~ α1δ1J ;β1δ2J 2 ζ~ α2δ1J ;β2δ2J @ζ~ 5 ; J 5 j: (11.137) j 2Δξj @ξ α;β

A substitution of Eq. (11.137) into Eq. (11.136a) gives n11;m11 γ n11;m11 5 γ n11;m 1 An11;m ; α;β α;β α;β Δðv1vÞα;β

(11.138)

and the difference equation describing the kinematic boundary condition on the free surface is expressed by 1:5Δζ~ α;β Δτ

n11;m

n11;m Δζ~ α;β

1:5ζ~ α;β

n11;m

5

2 γ n11;m α;β

n11;m11 5 ζ~ α;β

n11;m11 2 An11;m α;β Δðv1vÞα;β

2

n n21 2 2ζ~ α;β 1 0:5ζ~ α;β ; Δτ

(11.139)

n11;m 2 1:5ζ~ α;β :

In this equation, the increments of the wave elevation ζ~ and the fluid velocity vi are unknowns. Therefore Eqs. (11.134f) and (11.139) are solved by an iteration method at the same iteration time level to obtain all incremental values. That is, we continue the procedure from the previous wave height ζ~ to calculate the incremental ΔD from Eq. (11.134f) from which the corresponding velocity is substituted into Eq. (11.139) to obtain a new approximate wave  height ζ~ . Based on this new free-surface position, the incremental ΔD is calculated in iterations until it is negligibly small. The dynamic pressure condition in Eq. (11.117c) on the free surface is satisfied by introducing the prescribed pressure value p^0 at the mesh nodes on the free surface.

466

FluidSolid Interaction Dynamics

11.4.3.4 Dynamic equation of solid motion The dynamic equations with the boundary conditions of the solid given in Eqs. (11.122a)(11.122j) can be represented in matrix form: € 1 F ðθi ; θ_ i Þ 5 F ðθi ; F^ r ; T^ s Þ 1 f Σ ðθi ; σrs Þ; MU R S

(11.140a)

where  U 5 uoT

θT

T

 ; uoT 5 u01

u03 ;

u02

 θT 5 θ1

θ2

θ3 ;

M 5 diag½m; IðθÞ; m 5 diagðm; m; mÞ; ðIÞil 5 I il ðθÞ; h iT FR 5 0 0 0 θ_ T I 1 θ_ θ_ T I 2 θ_ θ_ T I 3 θ_ ; ðI i Þln 5 Iiln ðθÞ:

(11.140b)

The column matrix FS represents the corresponding structural force vector produced by the body forces F^ i and the traction force T^ i applied to the rigid body being rotated at θi . In general, the force vector FS is a function of the dis_ U; € tÞ, and includes stiffness, damping, and inertial placement, velocity, acceleration, and time, that is, FS 5 FS ðU; U; forces dependent on the motion of the rigid body. In most engineering applications, FS is independent of the acceleration of the solid, although it is included for completeness in the mathematical model. The column matrix FΣ represents the interaction force vector on the FSI interface Σ. This vector depends on the fluid stress tensor σij for viscous fluids or the dynamic pressure p for nonviscous fluids, the current position of the interface involving the rotation θi , and the current displacement of the solid. The column matrix FR represents the forces due to the rigid rotation θi and the rotation velocity θ_ i as given in Eq. (11.140b). From these considerations, Eq. (11.140a) can be rewritten to a general nonlinear equation: _ tÞ 1 f Σ ðU; σÞ; € 5 FS , R ðU; U; MðUÞU

F S , R 5 FS 2 FR :

(11.141)

This nonlinear equation is solved by using an incremental step-by-step approach. Let us assume that the solution at the discrete time t 5 nΔt, t is known and that the solution for the discrete time t 1 Δt 5 ðn 1 1ÞΔt is required, where Δt is a suitably chosen time increment. The application of the NewtonRaphson iteration method (see, e.g., Bathe, 1982, 1996; Xing, 2015b) to Eq. (11.141) gives the iterative equations: € 1 CΔU _ 1 KΔU 5 FS , R ðUt ; U _ ; t 1 ΔtÞ 1 f Σ ðUt ; σt1Δt Þ 2 MðUt ÞΔU € ; MðUt ÞΔU t

t

t

(11.142)

where 0

1t @F S , R A; Ct 5 @ _ @U

1t 1t 0 1t 0 € @F @f @ðM UÞ S , R Σ A; A1 @ A1 @ Kt 5 @ @U @U @U 0

(11.143)

Ut1Δt 5 Ut 1 ΔUt : A direct integration method (see Bathe, 1982, 1996) can be used to solve this equation to obtain the displacement, _ ΔU: € velocity, and acceleration increments, ΔU; ΔU;

11.4.3.5 Numerical solution _ and acceleration U € of the structure at time t 1 Δt are In the previous equation, the displacement U, velocity U, unknowns. The fluid stress σðt 1 ΔtÞ or pressure pðt 1 ΔtÞ is also unknown. Eq. (11.142), describing the motion of the structure, and Eqs. (11.134f) and (11.139), describing the fluid dynamics, are coupled. To solve these coupled equations describing the nonlinear dynamics of the fluidstructure interactions, the partitioned iterative approach is used based on the computation flowchart shown in Fig. 11.10. The main steps of the iteration have been given in Section 11.3.2. Now we have derived all the numerical equations for simulating nonlinear rigid FSI problems in term of computers. To illustrate the developed nonlinear mathematical model and numerical iterative scheme, a selection of examples is chosen to highlight the breadth of application of the developed approach in the following subsections. In many ways, these simplistic examples illustrate how these fluidstructure interaction mechanisms couple the various components of the posed problem and the modeling of large excursion motions.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

467

x

Free surface

Fluid

Yf

yf

M

Y0

y0 K

o FIGURE 11.15 Fluidmassspring interaction system with a free surface (Xing et al., 2003).

11.4.4

A fluidmassspring interaction system

Fig. 11.15 illustrates the fluidmassspring interaction system under examination. It consists of a nonviscous fluid column of height h 5 Yf 2 Y0 with a free surface, a piston of mass M, and sectional area A with a spring of stiffness K. The fluid is assumed to be incompressible and of density ρf . In the static equilibrium condition, y 5 Y0 and y 5 Yf denote the bottom and top positions of the fluid column which, at any time t during motion, move to positions y0 and yf, respectively. Y^ represents the position of the top face of the piston when the spring is in its unstretched state. For this problem, it is assumed for simplicity that the rigid body translates along the y axis with no rotational motion. Its motion can therefore be fully determined by displacement y0. The dynamic equation of the rigid body is given by M

d 2 y0 ^ 1 Mg 5 2 pA: 1 Kðy0 2 YÞ dt2

(11.144a)

In the static equilibrium state, d2 y0 =dt2 5 0; p 5 ρf hg, and y0 5 Y0 , which, when substituted into Eq. (11.144a), gives ^ 1 Mg 5 2 ρf Ahg: KðY0 2 YÞ

(11.144b)

The subtraction of Eq. (11.144b) from Eq. (11.144a) yields the dynamic equation of motion of the mass relative to its static equilibrium position, that is,

468

FluidSolid Interaction Dynamics

M

d2 y~0 ~ 1 K y~0 5 2 Aðp 2 ρf hgÞ 5 2 Ap: dt2

(11.144c)

where y~0 5 y0 2 Y0 and p~ 5 p 2 ρf hg represent the dynamic displacement of the mass and the dynamic pressure of the fluid excluding the hydrostatic pressure ρf hg, respectively. A moving coordinate axis, x, fixed on M, is chosen to describe the motion of fluid of height ζ~ 5 yf 2 y0 relative to the mass. The governing equations describing the fluid motion now take the following forms: ρf v;t 5 2 p~;x 2 ρf y€~0 ;

~ 0 , x , ζ;

v;x 5 0;

ζ~ ;t 5 v;

p~ 5 0;

~ x 5 ζ;

v 5 0;

x 5 0:

(11.144d)

~ 5 h; ζð0Þ

(11.144e)

The initial conditions of the system are given by y~0 ð0Þ 5 0;

y_~0 ð0Þ 5 V0 ;

~ 0Þ; vðx; 0Þ 5 0 5 pðx;

where V0 is the velocity of the mass M at time t 5 0. Let us assume that h and V0 represent the characteristic length and velocity, respectively, from which the characteristic time h=V0 and pressure ρf V02 are derived. The nondimensional forms of Eqs. (11.144a)(11.144e), using the same notation, are, respectively, ~ ~ tÞ; y~0 ð0Þ 5 1; v;t 5 2 p~;x 2 y€~0 ; 0 , x , ζ; (11.145) y€~0 1 ω2 y~0 5 2 ψpð0; pffiffiffiffiffiffiffiffiffiffiffi where ω 5 ðh=V0 Þ K=M is the nondimensional natural frequency of the massspring system, and ψ 5 ρf Ah=M represents the mass ratio of the fluid and solid. The nondimensional forms of the other equations in Eqs. (11.144a) (11.144e) are the same as their original derivations. On solving Eq. (11.145) taking into account the stated boundary and initial conditions, it follows that v 5 0;

p~ 5 ð1 2 xÞy€~0 ;

ζ~ 5 1;

~ y€~0 1 ω2 y~0 5 0: ð1 1 mÞ

(11.146a)

Here, m~ represents an additional mass coefficient of the fluid to the mass M (i.e., the added mass). The corresponding theoretical solution is y~0 5

1 sinðωsf tÞ; ωsf

(11.146b)

where ωsf is the natural frequency of the FSI system, that is, ω ωsf 5 pffiffiffiffiffiffiffiffiffiffiffiffi: 1 1 m~

(11.146c)

Figs. 11.16 and 11.17 show a comparison of the theoretical solution of the nondimensional displacement y~0 in ~ tÞ in Eq. (11.146a) with the numerical solutions obtained by the FSI iteration approach Eq. (11.146b) and pressure pð0; described here. It can be seen that there is very good agreement between the solutions. The numerical results for the case ofpffiffiffiffiffi a viscous fluidmassspring system with Reynolds number Re 5 hV0 =ν 5 100 and Froude number Fr 5 V0 = hg 5 0:32 are also shown in these two figures. As observed, the introduction of viscosity into the fluid slightly increases the period of vibration of the system. It is noted that in this mathematical model, the stress σij replaces the pressure in the numerical scheme of study as represented in Eq. (11.142). Therefore, for this viscous case, the fluid is contained in a finite thickness, 2D domain having a grid density idealization

FIGURE 11.16 Comparison of the nondimensional dynamic displacement of the mass for viscous and nonviscous fluids (Xing et al., 2003).

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

469

FIGURE 11.17 Comparison of the nondimensional dynamic fluid force acting on the mass for viscous and nonviscous fluids (Xing et al., 2003).

FIGURE 11.18 Wave height ζ~ 2 1, relative to the free-surface position in the nonviscous case, of the free surface of the fluid in the viscous case. x1 represents a nondimensional coordinate in the horizontal direction (Xing et al., 2003).

of 41 3 71 3 3 elements in width, height, and thickness directions. In the latter direction perpendicular to the plane of the paper, a symmetry condition is used. In the calculation, a time step Δt 5 0:01 is chosen and the properties of the system are represented by a spring stiffness K 5 4π2 , fluid density ρf 5 1, and mass M 5 1. Fig. 11.18 shows the wave height of the free surface of the fluid in the viscous case. This height is calculated relative to the free-surface position in the nonviscous case, which is the same as the displacement of the mass as shown in Fig. 11.16, since the fluid is assumed incompressible. It is observed in Fig. 11.17 that there exists a jump at t 5 0 in the curve of the dynamic force for the viscous case. This arises because the exact pressure and velocity distributions in the water corresponding to the initial position and velocity in Eq. (11.144e) of the mass cannot be absolutely defined, and therefore the initial conditions of fluid pressure and velocity are assumed in the calculations. This implies that the initial conditions of the mass position and fluid pressure in Eq. (11.144e) do not satisfy the dynamic equilibrium equations of the system at time t 5 0. This phenomenon is often found in the literature (see, e.g., Yeung and Ananthakrishnan, 1992).

11.4.5

A spring-supported fluidrigid-dam interaction system

^ acting Fig. 11.19 illustrates a spring-supported fluidrigid-dam interaction system excited by a prescribed torque MðtÞ on a rigid dam located at a distance OA 5 XA . This system also models a simple flap or paddle wave maker used in towing tank experiments (see, e.g., Gentaz et al., 2000). The mass center of the dam of total mass M is located at point B at a distance H in the direction S from point A. A fluid of density ρf occupies the movable domain Ωf with free surface Γ f and fixed boundaries Γ 1v and Γ 2v . The angle θ0 denotes the initial rotation of the spring from its nondeformed state. In static equilibrium, the depth of fluid is ζ 0 , whereas in the dynamic case, the fluid depth of the free surface is represented by the unknown function ζðx1 ; tÞ. The equation of dam motion is given by ðζΣ σij ν i ν j x2 dx2 ^ 1 ; (11.147a) I θ€ 1 Kðθ0 1 θÞ 2 MgH sin θ 5 MðtÞ cos2 θ 0

470

FluidSolid Interaction Dynamics

X2

vi

S

θ Γf

^

^ p=0

M ΩS B

Γv1 Ω

ζ(X ,t)

ζ0

1

Σ

f

Mg

ζ

Σ

^ fi Γv 2

o

A

x

1

FIGURE 11.19 Spring-supported fluidrigid-dam interaction system or simple wave maker (Xing et al., 2003).

where I represents the moment of inertia of the rigid dam around its simply supported base point A, and K denotes the spring stiffness constant. For the nonviscous fluid case, σij ν i ν j 5 p; in which vi denotes the unit normal vector pointing out from fluid inside to outside, as shown in Fig. 11.19. It is assumed that at the initial time t1, the system is in its static equilibrium state with θ 5 0. The relation between the initial fluid depth ζ 0 and angle θ30 0 in the static equilibrium state of the system is derived from Eq. (11.147a) with M^ 5 0; θ€ 5 0, and p 5 ρf gðζ 0 2 x2 Þ, giving Kθ0 5 ρf g

ðζ0

ðζ 0 2 x2 Þdx2 5

0

ρf gζ 30 : 6

(11.147b)

The characteristic p length of fluid ζ 0 and natural vibraffiffiffiffiffiffiffiffi and time chosen as nondimensional parameters are the depth 2 tion p period T 5 2π I=K , respectively, defining the Reynolds number Re 5 ζ =ðνTÞ and Froude number 0 ffiffiffiffiffiffiffiffiffiffi Fr 5 ζ 0 =g=T: The nondimensional form of Eq. (11.147a) is ^ MðtÞ MgH sin θ 5 1 Mf ; θ€ 1 4π2 θ 2 I=T 2 I=T 2

(11.148a)

where Mf 5

1 I=T 2

ðζΣ 0

px2 dx2 2 4π2 θ0 cos2 θ

(11.148b)

represents the dynamic moment about the point A produced by the fluid pressure. In the calculations, the following parameters are used: Fr 5 0.178, Re 5 200, H 5 1.0, θ0 5 0:1; ζ 0 5 1:0; M 5 1.0, I 5 4MH 2 =3; ρf 5 1:0; g 5 9.8, Δt 5 0:02:

11.4.5.1 Free vibration The dam is moved to an angle θ 5 θ1 5 0:1 rad by an external force. When this force is removed, the system is assumed to remain in its free-vibration state. Figs. 11.2011.22 show the dynamic moment Mf on the dam produced by the fluid pressure, the dynamic rotation angle θ, and the wave height ζ 2 1 of the free surface.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

471

FIGURE 11.20 Dynamic moment Mf on the dam produced by the fluid pressure in the free-vibration case (Xing et al., 2003).

FIGURE 11.21 Dynamic rotation angle θ of the dam in the free-vibration case (Xing et al., 2003).

FIGURE 11.22 The wave height ζ 2 1 of the free surface in the free-vibration case.

11.4.5.2 Forced vibration The dam or wave maker is excited by an external nondimensional moment   ^ T 2 MðtÞ T2 7t 5 sin ; I 3 2I

(11.149)

from its static equilibrium state, θ 5 0, where t denotes a nondimensional time. Figs. 11.2311.25 show the dynamic moment Mf on the dam produced by the fluid force, the dynamic rotation angle θ, and the wave height ζ 2 1 of the free surface, respectively.

11.4.6

Two-dimensional rigid body floating in a free-surface fluid

11.4.6.1 Free vibration Fig. 11.26 shows a rigid uniform square cylinder of infinite length, total mass M per unit length, and breadth L and floating in an infinite unbounded fluid with a free surface. Let us assume that, in the static equilibrium state, the center

472

FluidSolid Interaction Dynamics

FIGURE 11.23 Dynamic moment Mf on the dam produced by the fluid pressure in the forced vibration case (Xing et al., 2003).

FIGURE 11.24 Dynamic rotation angle θ of the dam in the forced vibration case (Xing et al., 2003).

FIGURE 11.25 Wave height ζ 2 1 of the free surface in the forced vibration case (Xing et al., 2003).

C of the cylinder is located at point O, where an absolute coordinate-axis system O 2 y2 y3 is fixed. Therefore the mass per unit length of the cylinder is M 5 ρf L2 =2 equaling the weight of the water displaced. At time t1 5 0, the center C of the cylinder is statically displaced to the position y3 5 0:15L by an external force. After this external force is released, the cylinder undergoes a freevibration interaction with the fluid. A moving axis-coordinate system Cp2ffiffiffiffiffiffiffiffiffiffi x2 x3 fixed at the center C of the cylinder is used in this calculation. The characteristic time T 5 2π=ω , ω 5 2g=L , 0 0 pffiffiffiffiffiffiffiffi Reynolds number Re 5 L2 =ðTνÞ 5 200, Froude number Fr 5 L=g=T 5 0:225, and time increment Δt 5 0:02 are chosen. Figs. 11.2711.29 show the dynamic force applied to the cylinder by the fluid, the displacement of the center of the cylinder, and the wave height disturbance of the free surface. In Fig. 11.29, free surface wave heights vary very quickly near the cylinder and seem to be very large. These effects are due to the imposed no-slip condition for the viscous fluid

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

473

FIGURE 11.26 Free vibration of a 2D rigid square cylinder floating in a free-surface fluid (Xing et al., 2003).

23 22 21 CL 20 19 18 17 0

0.5

1.0

1.5 t

2.0

2.5

3.0

FIGURE 11.27 Dynamic force coefficient CL 5 F=ð0:5ρf L4 =T 2 Þ on the cylinder produced by the fluid pressure. F represents the total fluid force acting on the cylinder (Xing et al., 2003).

FIGURE 11.28 Dynamic displacement of the center of the cylinder (Xing et al., 2003).

on the surface of the cylinder. Therefore fluid particles within a very thin boundary layer near this interface are moving at the speed of the cylinder.

11.4.6.2 Forced vibration Now we investigate the forced motion of this system excited by the nondimensional sinusoidal force f 5 4 sinð2παtÞ applied vertically at the mass center of the cylinder in its static equilibrium position at time t 5 0. The values α 5 ω=ω0 5 0:5; 0:99; 2:0, where ω denotes the forced frequency. Fig. 11.30 shows the time histories of the

474

FluidSolid Interaction Dynamics

FIGURE 11.29 Wave height ζ of the free surface of the viscous fluid (Xing et al., 2003).

– ω = 0.5ω0 × ω = 0.99ω0 º ω = 2.0ω0

0.0

0.5

1.0

1.5

2.0

2.5

t

y3

0.0

3.0

–0.25 FIGURE 11.30 Time histories of the dynamic response of the mass center of the cylinder excited by a nondimensional sinusoidal force f 5 4 sinð2παtÞ (Re 5 1000) (Xing et al., 2005).

nondimensional displacement at the mass center. It is observed from this figure that the amplitudes of the displacements of the mass center at frequency ratio α 5 0.5, 0.99 are nearly the same. As was demonstrated by Xing et al. (2003), the frequency ratio of the damping system is 1/1.35D0.74, which is the midpoint between 0.5 and 0.99. Because the frequency ratio 2.0 corresponds to the case of a much higher force frequency than the natural frequency of the floating cylinder, the amplitude of displacement response is very small. Similarly, the time histories of the coefficient of total fluid force acting on the cylinder can be obtained, but they are neglected here.

11.4.7

Prescribed motion of a rigid cylinder floating in the water

The dynamic behavior of a ship traveling in an irregular seaway is commonly characterized by the symmetric (heave, pitch) and antisymmetric (sway, yaw, roll) motions of the vessel (see, e.g., Bishop and Price, 1979). To aid in understanding the fluidstructure interaction, studies have been undertaken relating to the prescribed heave or roll motion of a cylinder floating in a free surface (see, e.g., Ursell, 1949a,b; Tasai and Koterayama, 1976; Chung, 1977; Yeung and Ananthakrishnan, 1992). Because the motion of the cylinder and therefore the boundary conditions on the wetted interfaces are prescribed, the relevant problem in fluid dynamics is defined, although the flow around the cylinder and description of the surface waves generated on the free surface remain unknown. From this consideration, the problem of a prescribed heave motion of a 2D rigid body floating in a free surface fluid, as shown in Fig. 11.31, is not an FSI problem. However, the simulation of this example was required by the referees of the journal, to which the paper by Xing et al. (2003) was submitted, in order to compare the results with the available publication at that time. Therefore the problem studied by Yeung and Ananthakrishnan (1992), as well as by Yeung and Liao (1999) as described in Fig. 11.31, is reexamined to validate the proposed mathematical model, as well as to illustrate its engineering application.

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

475

FIGURE 11.31 Cylinder heaving in a viscous fluid (Xing et al., 2003).

As shown in Fig. 11.31, to model this system, a fixed coordinate system is chosen with the x axis lying in the mean water level and the y axis pointing upward. The equilibrium draft of the rectangular body is denoted by d and its beam by B. The free surface is denoted by Γ f , the body contour by Σ, and the open boundary by Γ N . The vertical oscillation of the cylinder is given by yðtÞ 5 a sinðωtÞ, where a denotes the amplitude of the prescribed heave motion and ω its frequency. In the calculation, the time is scaled by 1=ω, length by B, velocity by ωB, and force by ρf ω2 B3 . The parameters chosen to illustrate the application of the proposed mathematical model are ω 5 1, Reynolds number Re 5 1000.0, Froude number Fr 5 1:0, a=B 5 0:2; d=B 5 0:5, and time increment Δt 5 0:02. Fig. 11.32 shows a selection of velocity-vector plots at nondimensional times t 5 14.2 and t 5 15.8, at which the locations of the cylinder bottom and its velocities are y 5 2 0:300; 2 0:518, and y_ 5 2 0:013; 2 0:200, respectively. In this figure, the motion of the free surface and the vortices generated by the motion of the cylinder are clearly illustrated. Fig. 11.32 clarifies the pressure and viscous-stress contributions to the heave force FðtÞB=a. Here, Fp and Fs denote the components of the heave force produced by the dynamic pressure and viscous shear stress, respectively, and y(t) is the prescribed motion of the cylinder. A calculated amplitude of Fp 5 0:51 is comparable with values 0:47 and 0:60 deduced from graphical data presented in Yeung (1975) and Yeung and Ananthakrishnan (1992), respectively. These investigators also produced harmonic time histories similar to the ones displayed in Fig. 11.33 predicated by the present study. This evidence provides further validation of the generalized numerical method developed.

11.4.8

Flows around a bluff body

11.4.8.1 Overset grids This example involves the multizone overlap grid technique, which provides a conceptually simple method for domain decomposition. For instance, the background grid is constructed to cover the main body element, and the overset grid overlays the background grid to improve geometric flexibility and so effectively captures the interesting features of the flow field as described by Kao et al. (1994), Koomullil and Soni (1999), and Xing et al. (2002, 2003, 2005). The overset grid may be composed of several grid blocks. Each of these blocks may be treated in exactly the same way as those grid blocks of the background grids, or, alternatively, different solution methods in the different subdomains may be employed. In addition to the boundaries of the computational domain, subgrids have intergrid boundaries with their neighboring donor subgrids but may also contain holes. Therefore proper boundary conditions are needed at the intergrid boundaries, and the hole points are blanked or excluded from the solutions of the NavierStokes equations. More details regarding the overset grid technique are described by Baysal et al. (1991) and Kao et al. (1994). The information transfer between the overlapping grids is of primary importance for improved efficiency of the overset technique. Here, a trilinear interpolation function, which has proved more effective than a Taylor series expansion method, is used to obtain the boundary values of the flow field variables φ in the overset grid derived from the background grid. This means φ 5 a1 1 a2 ξ 1 1 a3 ξ 2 1 a4 ξ 3 1 a5 ξ 1 ξ 2 1 a6 ξ 1 ξ 3 1 a7 ξ 2 ξ 3 1 a8 ξ 1 ξ 2 ξ 3 ;

(11.149)

476

FluidSolid Interaction Dynamics

FIGURE 11.32 Velocity-vector plots and free-surface disturbance (y 5 0, calm water surface) at nondimensional times (A) t 5 14.2 and (B) t 5 15.8 (Xing et al., 2003).

FIGURE 11.33 Pressure and viscous-stress contributions to the heave force FðtÞB=a. Fp and Fs denote the components of the heave force produced by the dynamic pressure and viscous shear stress, respectively, and y(t) is the prescribed motion of the cylinder (Xing et al., 2003).

where 0 # ξj # 1; ðj 5 1; 2; 3Þ are mesh coordinates, and the coefficients ai for i 5 1 to 8 are determined from the known values of variable φ at the eight base vertices of the hexahedron surrounding the target point. Since trilinear interpolation can be used only on cubes and the hexahedron formed around a target point is generally warped due to the curvilinearity of a body-fitted grid, it is therefore necessary to map a warped hexahedron to a unit cube by means of an

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

477

isoperimetric mapping. The relations between the coordinates ðx1 ; x2 ; x3 Þ of the target point and those ðξ1 ; ξ2 ; ξ3 Þ of the mesh point can be expressed in the following manner: xi 5 a1i 1 a2i ξ 1 1 a3i ξ2 1 a4i ξ 3 1 a5i ξ 1 ξ 2 1 a6i ξ1 ξ3 1 a7i ξ 2 ξ3 1 a8i ξ1 ξ 2 ξ 3 ;

(11.150)

aji

where the coefficient is determined by the coordinates of the eight vertices of the hexahedron similar to the calculation of the coefficients ai in Eq. (11.149). A Newton iteration method is used to derive the coordinates ðξ1 ; ξ2 ; ξ3 Þ corresponding to any point Pðx1 ; x2 ; x3 Þ in this cubic volume, and then the value of the field variable φ at point Pðx1 ; x2 ; x3 Þ can be obtained by Eq. (11.149).

11.4.8.2 The definition of the problem As a further validation example of the numerical method developed, numerical simulations were carried out to investigate the nature of flow around a 2D bluff body, as shown in Fig. 11.34, which is a benchmark problem in CFD. Durao et al. (1988) experimentally measured and Freitas and Runnels (1999) and Caughey (2001) numerically predicted flows past a fixed cylinder of square cross section. The experiments undertaken by Durao et al. (1988) were performed in a horizontal rectangular water tunnel, 0.120 m 3 0.156 m, with a 0.02 m 3 0.02 m obstacle set across the narrow dimension. The duct was extended to 1.56 m upstream and 0.44 m downstream of the obstacle. The square cross-sectional cylinder was centered between the two walls of the duct in a uniform flow velocity U0 5 0:68 m=s, giving a Reynolds number Re 5 LU0 =ν 5 14; 000, where L is the length of the cylinder side. In the numerical simulations undertaken, a square cylinder of the same size as used in the experiment by Durao et al. (1988) is chosen. However, in order to include the fluidbody interaction cases, the cylinder is tethered by a rigid cable of zero mass to a frictionless hinge on the axial centerline of the duct, and the other end of the rigid cable is fixed to the mass center of the cylinder. The computational configuration is shown in Fig. 11.34. The length of the cable is the same as the side length of the cylinder. The motion of the cylinder is described by two variables γ, θ, as indicated in Fig. 11.34, and the cable constrains the motion of the cylinder’s mass center to an arc of a circle of unit radius. At the same time, the cylinder is allowed to rotate around its mass center through angle θ. The characteristic length and velocity chosen in the calculations (and experiment) are cylinder length L(50.02 m) and uniform flow velocity U0 ð 5 0:68 m=sÞ. The Reynolds number is the same as the one used in the experiment, that is, Re 5 LU0 =ν 5 14; 000: The characteristic time, that is, the dimensionless unit of time in the numerical simulation, is L=U0 5 1=34 5 0:0294 second. The lift force and moment coefficients are given by Cf 5 F=ð0:5ρU02 L2 Þ and Cm 5 M=ðρU02 L3 Þ, respectively. In all computations, the mass per unit length of the cylinder is 1.5 kg/m. For the numerical predictions, the computational domain extends from five dimensionless lengths upstream and 20 downstream of the cylinder, so the total computational domain is of length 26 and width 7.8. This nondimensional width corresponds to the real width of the channel in the experiment. The system is simulated as a dynamical fluidcylinder interaction problem defined in similarly chosen experimental and numerical domains. The boundary conditions in the numerical calculations are assumed as follows. On the inflow boundary, the uniform flow velocity is given and the pressure extrapolated from inner to outer points of the domain. On the outflow boundary, all variables are set by a zero-normal gradient condition. On the fixed wall boundary, a no-slip condition is implemented and a zeronormal gradient for the pressure imposed. On the FSI interface, the velocity at every point on the wet surface of the rigid body equals the velocity at the same point in the viscous fluid field, and a zero-normal gradient for the pressure is imposed.

FIGURE 11.34 Definition of variables and sketch of cylinder and numerical simulation domain (Xing et al., 2005).

478

FluidSolid Interaction Dynamics

TABLE 11.1 Three cases to be calculated in the simulation.

Initial conditions

Geometrical conditions of the cylinder motion

Case 1

Case 2

Case 3

γ505θ

γ50

γ 5 π=4

θ 5 π=6

θ50

γ50

γ Variable

θ Variable

θ Variable

γ505θ

TABLE 11.2 Comparison of numerical results obtained by using different time increments, grid resolutions, and experimental results (γ 5 0 5 θ) (Xing et al., 2005). Time increment Δt

Background grid resolution

Overlap grid resolution

Drag coefficient time-average/ amplitude

Lift coefficient time-average/ amplitude

Frequency of vortex shedding (Hz)

0.1 and 0.2

81 3 41

73 3 73

2.06/0.16

0.00/1.96

4.53

0.2

81 3 41

99 3 99

2.23/0.19

0.00/2.07

4.72

0.2

121 3 71

99 3 99

2.18/0.20

0.00/2.11

4.72

0.01

81 3 41

73 3 73

2.03/0.10

0.00/1.78

4.66

0.05

81 3 41

73 3 73

2.06/0.12

0.00/1.82

4.62

0.05

121 3 71

73 3 73

2.01/0.14

0.00/1.62

4.63

0.05

81 3 41

99 3 99

2.15/0.11

0.00/1.98

4.87

Numerical Koobus et al. (2000)

1.97/Not available

Experimental Koobus et al. (2000)

2.052.23/Not available

Experimental

0.00/Not available

4.7

Durao et al. (1988)

In the numerical simulations undertaken in that study, the three cases listed in Table 11.1 are examined in the flowing manner. Case 1 investigates the solution for the wake flow past a fixed cylinder of square cross section. To examine the sensitivity of the proposed numerical scheme of study, Table 11.2 displays a comparison of findings using different time increments and grid resolutions to assess convergence of the numerical solution. This is undertaken at Reynolds number Re 5 14,000 in order to compare with the numerical and experimental results presented by Koobus et al. (2000) and the experimental findings of Durao et al. (1988). In general, the predicted results followed the trends observed in the experiments. The majority of the calculated drag coefficient time-average values lie within the range of experimental data reported by Koobus et al. (2000), although the predicted value of Koobus et al. (2000) is smaller than the currently predicted and experimental results. The calculated value of the frequency of vortex shedding lies within a 6 3% bandwidth of the experimental value of Durao et al. (1988), which is deemed acceptable. The time history traces of drag and lift coefficients (as shown in Figs. 11.3511.37) are similar to the regular curve patterns presented by Koobus et al. (2000). As observed, changes to time and mesh values cause relatively small changes in the overall results (though they

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

479

FIGURE 11.35 Time histories of the drag and lift coefficients for Case 1, Re 5 14,000, (Xing et al., 2005).

Fine grid Coarse grid 2.2

CD

2.1 2.0 1.9 1.8 100.0

102.5

105.0 Time

107.5

110.0

FIGURE 11.36 Time histories of lift coefficient (Case 1, Re 5 14,000) using a small time increment Δt 5 0:05 and different grid resolutions: dashed line for background grid 121 3 71 with overlap grid 73 3 73; solid line for background grid 81 3 41 with overlap grid 73 3 73 (Xing et al., 2005).

Fine grid Coarse grid 2

CL

1 0 –1 –2 90

95

100 Time

105

110

FIGURE 11.37 Time histories of drag coefficient (Case 1, Re 5 14,000) using a small time increment Δt 5 0:05 and different grid resolutions: dashed line for background grid 121 3 71 with overlap grid 73 3 73; solid line for background grid 81 3 41 with overlap grid 73 3 73 (Xing et al., 2005).

may influence detailed aspects of time histories, i.e., amplitudes), but the computational effort increases significantly as the time increment decreases in size and mesh resolution becomes finer. Cases 2 and 3 allow the cylinder to move with rotation angles γ and θ, as shown in Fig. 11.34 to investigate the fluidrigid cylinder interactions. Initially, the cylinder is held at prescribed deflections and then released in order to study its dynamic interaction with the flow. The time histories of the cylinder motion, the forces (i.e., drag and lift), and moments produced by the dynamical interactions between cylinder and viscous fluid flow are calculated. A more detailed computation of individual findings is now discussed.

480

FluidSolid Interaction Dynamics

FIGURE 11.38 Predicted velocity vectors of the fluid flow past the fixed square cylinder at different dimensionless time units 102, 104, 106, and 108 (Case 1, Re 5 14,000) (Xing et al., 2005).

11.4.8.3 Case 1: Fixed cylinder at γ 5 0 5 θ In this study, the cylinder is fixed on the centerline of the water channel (i.e., γ 5 0 5 θ) and not allowed to move. Fig. 11.35 displays a typical time evolution of drag and lift coefficients. In the first numerical experiment of the sensitivity exercise, for computational efficiency, we used a dimensionless time step of Δt 5 0:2 or 0.1, a background grid defined by 81 3 41 cells in the fluid flow and transverse directions, respectively, together with an overlap computational domain surrounding the square cylinder of grid resolution 73 3 73. The computational time increment for the first 80 steps is Δt 5 0:2 (nondimensional time), which was reduced to Δt 5 0:1 for the next 30 steps. The predicted oscillation period of the lift force is 7.5/34 5 0.22 second, giving an oscillation frequency of 4.53 Hz corresponding to the frequency of vortex shedding and lying within 3% of the experimental frequency 4.7 Hz reported by Durao et al. (1988). Because of the symmetry of the problem, the oscillation period of the drag force is half the oscillation period of the lift force. Fig. 11.38 displays the flow vectors calculated at different times in a period of vortex shedding. In these diagrams, the overlap computational domain of the base grid with the overset grid in the region of the obstacle is removed. Figs. 11.36 and 37 show the predicted time histories of the lift and drag coefficients, respectively, for time increment Δt 5 0:05 and a fine background grid resolution 121 3 71. For all the different time increments and grid resolutions considered in Table 7.2, the time-averaged drag coefficients and the frequencies of vortex shedding are not greatly changed, and they display reasonable agreement with the available experimental results. As expected for this symmetric flow, the time-averaged lift coefficient is zero, although minor amplitude differences are observed. Figs. 11.3511.37 show time histories of drag and lift coefficients displaying relatively regular patterns confirming the findings of Koobus et al. (2000) and are similar to those derived by Caughey (2001) at Re 5 120. They differ in form to the irregular patterns presented by Tamura and Kuwahara (1990), Kondo and Yamada (1995) at Re 5 10,000, and Okajima et al. (1992) at Re 5 1000. All time histories presented in the sensitivity analysis discussed in Table 7.2 produced relatively regular patterns, indicating that the sampling rate per cycle is not the cause of the difference. This leaves the mathematical model chosen and imposed boundary conditions as possible explanations. For example, in the FE model of Kondo and Yamada (1995), a circular fluid domain is used to represent the computational region, which had a diameter of 60 times the breadth of the cylinder. Although they give no clear definition of the boundary conditions adopted, they are very likely to be different from those used here, that is, the fixed wall boundary conditions on the channel walls. Franke and Rodi (1991) demonstrated that the influence of the mathematical model in deriving solutions is significant. To illustrate more clearly the convergence criterion adopted, time histories of the local fluid pressure and residual norms for Case 1 were calculated using a background grid 121 3 71, overlap grid 7 3 3 73, and time increment 0.05. Fig. 11.39 presents the time history of the local dynamic pressure at the left-top corner of the square cylinder, and Fig. 11.40 illustrates the pressure distributions along the top and bottom surfaces of the cylinder at times 25, 120, 140. Fig. 11.41 shows the time histories of L2 norms of the residuals of the continuity equation and momentum equations for the overlapping zone. These results further confirm that the numerical convergence is achieved not only for the global lift and drag coefficients, as

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

481

FIGURE 11.39 Time history of pressure calculated at the left-top corner of the square cylinder (Case 1, Re 5 14,000, time increment Δt 5 0:05, background grid 121 3 121 and overlap grid 73 3 73) (Xing et al., 2005).

FIGURE 11.40 Pressure distributions along the bottom and top surfaces of the square cylinder at different times (A) 25, (B) 120, and (C) 140 (Xing et al., 2005).

shown in Fig. 11.35, but also for the local dynamic variables. Although initially the L2 norms show higher values than in the remaining time, nevertheless, they satisfy the chosen convergence criterion.

11.4.8.4 Case 2: Cylinder rotating about its mass center at γ 5 0 line In this investigation, the mass center of the cylinder is held at γ 5 0, and the cylinder is not allowed to move from the centerline of the water channel. The cylinder rotates only about its mass center, and its motion is determined by evolution of the cylinder-flow interaction. The initial angle of rotation between the central lines of the cylinder and the duct is θ 5 π=6. The time increment and grid resolution used in this example is the coarse model described in Case 1. Figs. 11.42 and 11.43, respectively, display the time histories of the dynamic moment and angular displacement of the oscillating cylinder produced by the flowcylinder interaction. The oscillations in the moment are caused by variations

482

FluidSolid Interaction Dynamics

FIGURE 11.41 Time histories of the L2 norms of the residuals of the continuity equation (A), momentum equations in x direction (B), and momentum equations in y direction (C) for the overlapping zone (Xing et al., 2005).

Moment M

0.1

0.0

–0.1

–0.2 0

50

100 Time

150

200

FIGURE 11.42 Moment acting on the cylinder induced by the fluid flow field for Case 2, Re 5 14,000 (Xing et al., 2005)

in the flow field around the bluff body and are manifested in the angular displacement. Fig. 11.44 shows the velocity vectors at different times.

11.4.8.5 Case 3: Cylinder with 2 degrees of freedom of motion, γ and θ variables In this case, at time t 5 0, the cylinder is held off center of the water channel at a position defined by γ 5 π=4 and with zero rotation (i.e., θ 5 0). After release, the cylinder moves with 2 degrees of freedom represented by γ and θ. If the cylinder’s rotation is fixed (i.e., θ 5 0) and not allowed freedom to rotate, the case reduces to one similarly investigated by Freitas and Runnels (1999), but the length of cable may differ because this dimension is not defined in our study. Fig. 11.45 illustrates the time histories of the variables γ and θ after release of the cable-cylinder system calculated using the coarse model described in Case 1. Fig. 11.46 shows the corresponding drag and lift coefficients experienced

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

483

FIGURE 11.43 Angular position as a function of dimensionless time for the cylinder about its mass center for Case 2, Re 5 14,000 (Xing et al., 2005).

FIGURE 11.44 Predicted velocity vectors of the fluid flow past the rotating square cylinder at different dimensionless time units (A) 20, (B) 24, (C) 28, and (D) 32 for Case 2, Re 5 14,000 (Xing et al., 2005).

FIGURE 11.45 Predicted time histories of the angular variables angles of θ and γ as defined in Fig. 11.1 for Case 3, Re 5 14,000 (Xing et al., 2005).

by the cylinder created by interaction influences. Fig. 11.47 illustrates the predicted moments Mγ and Mθ about the attachment point of the cable on the centerline and the other at the center of the cylinder, respectively. A periodic wake is developed in a very short time, and the lift coefficient varies periodically. The period of oscillation of the lift force is approximately 6.6/34 5 0.194 second, which corresponds to a frequency of vortex shedding of

484

FluidSolid Interaction Dynamics

FIGURE 11.46 Predicted drag and lift forces acting on the cylinder induced by the fluid flow field for Case 3, Re 5 14; 000, Cd: drag coefficient; Cl: lift coefficient (Xing et al., 2005).

FIGURE 11.47 Predicted moments acting on the cylinder induced by the fluid flow field for Case 3, Re 5 14,000; Mγ and Mθ represent the moment about the fixed end of the cable and the center of the cylinder, respectively (Xing et al., 2005).

5.15 Hz, indicating that the dynamical interaction causes a higher frequency of vortex shedding. The high-frequency but low-amplitude oscillation in the drag is caused by variations of the periodic shedding of vortices around the cylinder and the oscillating rotation of the cylinder. The oscillation periods of the two moments are the same as experienced in the lift force predictions, and they correspond to changes of the angular displacements. The angular displacement γ decreases quickly due to the effect of the moment caused by the drag force. Predicted flow vectors, at various stages of the fluidcylinder interaction, are plotted in Fig. 11.48. Here the overlap part of the base grid with the overset grid in the region of the cylinder is removed to aid observation. The first diagram shows the flow at dimensionless time t 5 20 when the lift is near its maximum point at time t 5 20:4, and the moment Mγ is near zero at time t 5 20:2. It is seen that a clockwise vortex is formed behind the cylinder. The last diagram at time t 5 44 is nearly identical to the first one, demonstrating a degree of repetition in the interaction cycle. The four diagrams displayed in Fig. 11.48 were not chosen in a period of shedding of vortices because our intention here is to show the change of the cylinder’s position during the interaction process. It is observed from Figs. 11.46 and 47 that both angles γ and θ, as well as the moment coefficients Mγ and Mθ , are approximately out of phase. This is because an increase of angle θ, as shown in Fig. 11.34, provides a decrease in the angle of attack of the flow on the cylinder, and therefore the lift and moment coefficients decrease.

11.4.8.6 Remark As discussed previously, the method adopted in this investigation is based on a partitioned or iterative solution procedure in which the equations describing the motions of the fluid and structure are solved separately using dedicated solvers, and the data are exchanged at every time step until convergence of the iteration is reached. This information is transferred forward to the next time step (see Xing et al., 2003) and the numerical scheme of study is repeated. Due to the generalized and varying relation between the computational grid and the fluid itself, as represented by the fluid variables, a rigorous stability/convergence analysis cannot be performed (see, e.g., Hirt et al., 1974, 1997). In the calculations presented, the

Mixed finite elementcomputational fluid dynamics method for nonlinear fluidsolid interactions Chapter | 11

485

FIGURE 11.48 Predicted velocity-vector fields of the flow past the translating and rotating square cylinder having 2 degrees of freedom at four nondimensional times of (A) 20, (B) 28, (C) 36, and (D) 44 (Xing et al., 2005).

time step chosen in the numerical scheme of study provides a sufficient number of points to describe the oscillation behavior of the predicted dynamic variables defining the interaction process. For example, the vortex shedding oscillation frequency 4.53 Hz predicted in this simulation corresponds to an oscillation period of 7.5 nondimensional time units. The maximum time step Δt 5 0:02 corresponds to sampling the captured vortex-shedding period in 7.5/0.2 5 37.5 time steps. This has been found acceptable to meet the condition for stability and convergence in computations of unsteady viscous flows as discussed, for example, by Koobus et al. (2000). This is because the Jacobian of the residual is diagonally dominant due to the upwind difference process adopted in our implicit scheme. This greatly enhances the subiterative procedure and provides an improved performance which, alternatively, can be obtained using a small time step increment. In the water channel, the cylinder imposes a 13% blockage (Freitas and Runnels, 1999) to the flow, which results an increase in the mean velocity in the vicinity between the channel walls and the cylinder. The modeling of this practical channel wall (i.e., wind tunnel, circulating water channel, etc.) is achieved by implementing, as described previously, a no-slip condition for the velocity and a zero-normal gradient condition for the pressure on the fixed wall boundary and the FSI interface. The velocity at every point on the wet surface of the rigid body is required to equal the velocity at the same point in the viscous fluid field. A blockage correction will be needed for these results to be useful to someone interested in a cylinder in an infinite stream, as discussed by Freitas and Runnels (1999). The first mesh point in the boundary layer of the cylinder surface is taken at distant dimensionless lengths 1:19 3 1024 , which is of the order of the viscous sublayer thickness and satisfies the requirement as defined by Hirsch (1990) through the expression rffiffiffiffiffi y τ0 y1 5 D1B10; (11.151) ν ρ where τ 0 is the wall shear stress, and y denotes the distance to the wall. This study focuses on the motion of the cylinder, and therefore the grid spacing near the channel walls does not strictly satisfy Eq. (11.151) but could nevertheless be made to do so with increased computational effort. During the iterative process developed within the numerical scheme of study, the convergence criterion adopted relates to a small residual value (i.e., 10221023, say) derived in the mass/momentum conservation calculations. For the three cases under examination, it has been numerically demonstrated that iterative convergence is achieved by adopting such a criterion. To rigorously prove the convergence mathematically is beyond the scope of this investigation since here we demonstrate the applicability and versatility of the proposed theoretical model (see, e.g., Xing et al., 2003) and numerical scheme of study. Further validations and verifications are necessary to develop a generic algorithm to solve fluidstructure interaction problems, but in this study, where possible, comparisons with limited available experimental data are examined.