Mixed finite element—smoothed particle methods for nonlinear fluid–solid interactions

Chapter 12

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter Outline 12.1 Introduction 12.1.1 Limitations of grid-based methods for violent flows 12.1.2 Key ideas of meshfree particle methods 12.1.3 History and developments with applications 12.2 Fundamentals of smoothed particle hydrodynamics 12.2.1 Smoothed particle hydrodynamics interpolations 12.2.2 Particle approximation 12.2.3 Construction of smoothing functions 12.2.4 Smoothed particle hydrodynamics formulation for Navier
Stokes equations 12.2.5 Numerical techniques for fluid flows 12.2.6 Improved methods based on smoothed particle hydrodynamics

488 488 488 489 492 492 499 503 516 520 527

12.3 Meshfree Galerkin methods 12.3.1 Moving least square representing kernel interpolant 12.3.2 Shepard interpolant 12.3.3 Orthogonal basis for local approximations 12.3.4 Applications of the moving least square reproducing kernels 12.4 Mixed finite element—smoothed particle method for fluid
solid interaction problems 12.4.1 Generalized solution procedure 12.4.2 Modeling of fluid
solid interaction involving large rigid motions with small elastic deformation 12.4.3 Application examples

531 531 534 535 535 536 536

537 547

The proposed finite element (FE)
computational fluid dynamic (CFD) method for nonlinear fluid
solid interaction (FSI) problems in Chapter 11, Mixed finite element
computational fluid dynamics method for nonlinear fluid
solid interactions, is one of grid-based numerical methods in which various types of fluid solvers have been used for different situations such as potential flow models and mesh-based Navier
Stokes (NS) solvers. The potential flow model is an efficient and stable approach for many wave
structure interaction problems. However, it cannot handle the situations where breaking waves occur. The general CFD solvers based on the NS equation are capable of dealing with these violent-free surface motions, but the large free surface deformation has caused extreme mesh distortion, which can either terminate the computation or result in drastic deterioration in accuracy (Swegle et al., 1995). However, in marine engineering, these types of violent-free surface flow
structure interactions are very important and common, for example, slamming impacts on ships, sloshing in liquefied natural gas tanks, green water loading onto the upper structure on the deck, etc. Therefore with the rapid development of computational science and technology, scientists, and engineers have been seeking alternative numerical approaches to overcome the shortcomings of grid-based numerical methods. Since the 1970s, meshless numerical methods have been developed, in which the smoothed particle method (SPM) has been well formulated and widely practiced in simulating violent fluid flows. In this chapter, we present a mixed FE (MFE) method
SPM, in which the solid structure is modeled by a powerful nonlinear FE scheme, while the fluid motions are modeled by the SPM, to deal with nonlinear FSI problems involving violent fluid flows. To fully understand the involved SPM in this mixed scheme, there is a summary briefing on the history, characteristics, developments with applications, and the fundamental theory of SPM in Sections 12.1
12.3 before a discussion of the proposed MFE method
SPM in Section 12.4. Several publications present more details on the meshfree particle methods (MPMs), such as Liu and Liu (2003), Li and Liu (2004), and Atlury (2004); readers may refer to these books for further detailed information.

Fluid
Solid Interaction Dynamics. DOI: https://doi.org/10.1016/B978-0-12-819352-5.00012-4 © 2019 Higher Education Press. Published by Elsevier Inc. All rights reserved.

487

488

Fluid
Solid Interaction Dynamics

12.1 12.1.1

Introduction Limitations of grid-based methods for violent flows

As previously discussed for linear FSI problems in Chapter 7, Finite element models for linear fluid
structure interaction problems, and for nonlinear FSI problems in Chapter 11, Mixed finite element
computational fluid dynamics method for nonlinear fluid
solid interactions, the proposed numerical methods require the continuum system to be divided into small connected elements or volumes, as in FE method (FEM), finite volume method (FVM), or finite difference method (FDM). These elements, called mesh or grid, allow the governing partial differential equations (PDEs) to be converted into algebraic equations. Today, these grid-based FEMs in solids and CFDs in fluids in computational mechanics—a far-reaching field extending from basic science to applied research to applications—play a prominent role in technological advancement. Despite great success, grid-based numerical methods suffer from inherent difficulties in many aspects that limit their application to a number of problems, especially nonlinear FSI problems with violent flows. In grid-based numerical methods, mesh generation in the problem domain is a prerequisite for numerical simulations. For the traditional mesh-based CFD approaches, constructing a regular grid for irregular or complex geometry has never been an easy task, and it usually requires additional complex mathematical transformation that can be even more expensive than solving the problem itself. Determining the precise locations of the inhomogeneities, free surfaces, large deformable boundaries, and moving FSI interfaces within the fixed Eulerian grid is also a formidable task. The grids could be generated or rezoned to fit these dynamic shapes of moving boundaries or interfaces by using most popular approaches, such as marker and cell (McKee et al., 2008), volume of fluid (Hirt and Nichols, 1981), level set (Osher and Fedkiw, 2001), and constrained interpolation profile (Hu and Kashiwagi, 2004; Zhu, 2006). These types of moving surface
capturing approaches can be very complicated, tedious, and time-consuming, which may introduce additional inaccuracy into the solutions. The difficulties and limitations of the grid-based methods are especially evident in simulating hydrodynamic problems with breaking waves, such as explosions, slamming, green water, and high-velocity impacts; tracing the millions of flying water particles caused by an explosion, for instance, is impossible on the fixed Eulerian grid. The motions of the FSI system dominated by physical governing equations are described in either Lagrangian or Eulerian terms. The Lagrangian description is traditionally used in solid mechanics, and the corresponding typical numerical scheme has been well-known as the FEM. The grid or mesh is fixed on the material and moves/deforms with the material. The physical variables defined at each material point and developing over time can be calculated and recorded. The updated Lagrangian grids and mesh rezoning techniques have been effectively adopted by the FEM codes to simulate the nonlinear problems of solid structures, such as penetrations, in which tracking the large elastic/plastic deformation is not as difficult as in the case of breaking waves of fluids. However, in fluid mechanics, the Eulerian description is often used, and the related numerical grid or mesh is fixed in the space, making it difficult to fit with the moving Lagrangian grid in solids on the FSI interfaces for nonlinear problems. To overcome the drawbacks of these two descriptions, the arbitrary Lagrangian
Eulerian (ALE) description was developed, as discussed in Chapter 11, Mixed finite element
computational fluid dynamics method for nonlinear fluid
solid interactions. In ALE formulations, the fluid grid is allowed to move in an arbitrary velocity, so that on the FSI interface, the fluid grid velocity can be updated to fit the material solid grid motion to satisfy the coupling conditions. Unfortunately, this updated ALE scheme still suffers from inherent difficulties in simulating nonlinear FSI problems with violent fluid motions, such as breaking waves and separations between solids and fluids. Numerical experiments have revealed two main reasons for these difficulties: one is due to different motion descriptions for fluids and solids, and the other is the grid-based numerical scheme. Therefore using the same motion description for fluids and solids and giving up the grid-based numerical methods could help in overcoming this inherent difficulty. As we have seen, the Lagrangian description has been very successful in simulations of solid motions, and there are many very powerful FE analysis (FEA) programs, so it would be desirable if the Lagrangian method could be used in both solid and fluid motion descriptions. This implies that we need to develop a numerical scheme to model fluid motions without grid requirements and based on the Lagrangian description. The meshfree, Lagrangian particle methods that were developed in this background required numerical simulations of complex systems in science and engineering.

12.1.2

Key ideas of meshfree particle methods

The key idea of MPMs is to employ a set of finite numbers of discrete particles to represent the state of a system and to trace the motion of the system. Each particle can be either directly associated with one discrete physical particle or

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

489

generated to represent a part of the continuum domain of the problem. The sizes of particles can range from very small nano or microscale to macroscale and even to astronomical scale. For continuum dynamics, each particle has a set of field variables such as mass, momentum, energy, etc. The evolution of the dynamic system is governed by the conservation laws. Depending on the mathematical models, MPMs can be deterministic or probabilistic. The particle evolution in the deterministic MPMs is totally governed by physical laws with regard to the initial and boundary conditions, and the solution of the problem is theoretically precisely determined. For probabilistic MPMs, the corresponding solutions have a probabilistic nature governed by statistical principles. MPMs fall into the category of Lagrangian methods, in which the particles describing the physical system move in the Lagrangian frame according to the conservation laws and the internal interactions, as well as external force derivations, and thus they evolve into its new state in time. The inherent Lagrangian nature of MPMs provides an essential basis for its modeling in fluids in combination with FEMs in solids to deal with nonlinear FSI problems with violent flows. The important and essential characteristics of MPMs lie in using discretized particles with no fixed connectivity to represent the problem domain and the related physical variables, as well as approximating the functions, derivatives, and integrals, etc. over the particles under the Lagrangian description rather than over the mesh in the Eulerian description. Since these characteristics are essential in nature, MPMs can overcome the previously mentioned difficulties and provide the following advantages over conventional grid-based numerical methods: (1) they make it relatively simple to solve problems with complex geometry; (2) particle refinement is expected to be easier than the mesh refinement; (3) large deformation is easier to deal with; (4) it is more convenient to trace particle motions and the time history of material points and to obtain the entire system features under the Lagrangian description and therefore to identify and trace moving boundaries with breaking waves. Similarly to the numerical process used in the grid-based approaches, a typical simulation for a defined problem by governing equations and initial and boundary conditions based on MPMs involves creating particles by domain discretization techniques, generating numerical equations
based particle techniques, and solving the resultant algebraic or ordinary differential equations (ODEs).

12.1.3

History and developments with applications

The origin of meshfree methods can be traced back a few decades, but it was not until after the early 1990s that substantial and significant advances were made in the field. The initial interest in developing meshfree methods was to develop meshfree interpolant schemes for FEM, such as the techniques for adaptive refinement and moving objects in fluids. The development of meshfree interpolants by Nayroles et al. (1992), Belytschko et al. (1994a,b, 1995a,b), and Liu et al. (1993a,b, 1995a,b,c) has led to a comprehensive understanding of meshfree discretization, associated variational formulations, and the birth of a new class of meshfree methods and partition of unity methods. During the past decade, meshfree technology has achieved many successes and attracted much attention due to its applications to computations for failure mechanics and multiscale nanotechnology. Some special issues, devoted to MPMs, have been published in the major international journals in computational mechanics, such as Hughes et al. (1996) and Chen and Liu (2000). Since then, there have been many review papers as well as important books in this field, such as the earlier ones by Monaghan (1982, 1985, 1988, 1992), Benz (1989, 1990), Li and Liu (2002, 2004), Liu (2002), Babuska et al. (2003), and Liu and Liu (2003). In general, there are two main computational modeling techniques in MPMs: one is the meshfree Galerkin methods (MGMs) and the other a direct approximation of the strong form of a PDE.

12.1.3.1 Meshfree Galerkin methods The MGMs are based on various Galerkin weak formulations of PDE to generate the numerical equations for problems. In these approaches, a field variable u for a particle located at x within the problem domain can be approximated using the information on the particles within the domain of the particle at x in the following interpolation form: uðxÞ 5

N X

φI ðxÞuI ;

(12.1)

I51

where N is the number of particles within the support domain of the particle at x, uI denotes the field variable at particle I, and φI is the shape function at particle I. The shape functions can then be used to generate a set of discretized system equations based on the weak form of the problem. The publications for the MGMs are too many to be listed here

490

Fluid
Solid Interaction Dynamics

because various PDE weak forms can be used. A typical representative may be the diffuse element method (Breitkopf et al., 1998, 2000; Nayroles et al., 1992). Many other different examples with publications, such as the element-free Galerkin method, the reproducing kernel particle method (Jun et al., 1998), the h
p cloud method, the partition of unity method, the meshless local Petrov
Galerkin method (Atluru and Zhu, 1998, 2000a,b; Atluri et al., 1998a,b, 2000; Lin and Atluri, 2001), the finite sphere method, the free mesh method, the moving particle FEM, the natural element method and reproducing kernel element method, and others may be consulted in Liu and Liu (2003) and Li and Liu (2004). MGMs have been successfully used in solving many difficult engineering problems. One of the important application examples is using its flexibility in interpolation to simulate crack growth and propagation problems (Belytschko et al., 1994a,b, 1995a,b; Daux et al., 2000). The visibility criterion was developed in modeling discontinuous fields (Belytschko et al., 1995b; Krysl and Belytschko, 1996; Dolbow et al., 2000). A series of publications by Li and Liu (1996a,b, 1997, 1999a,b, 2002, 2004) have shown that the MGM has the strong ability to sustain large mesh distortion.

12.1.3.2 Smoothed particle hydrodynamics Smoothed particle hydrodynamics (SPH) directly approximates the strong form of PDE in the discretization equation by using a special collocation technique. In SPH, the field variable f for a particle located at x is represented by the integral form ð f ðxÞ 5 f ðx0 ÞWðx 2 x0 ; hÞdx0 ; (12.2) Ω

where W is the smoothing kernel function with its parameter h or kernel function, and Ω is the volume of integral containing x. The derivatives of the field variable are also represented by the corresponding integrals to be used to obtain the numerical particle equations transformed from the PDEs of the problem. SPH is a meshfree, Lagrangian particle method for modeling fluid flows that has been largely developed in the last three decades since its first invention to solve astrophysical problems in three-dimensional (3D) open space (Gingold and Monaghan, 1977; Lucy, 1977). Several review papers and comprehensive books on SPH have been published, for example, Benz (1989, 1990), Monaghan (1992, 2012), and Liu and Liu (2003), and Li and Liu (2004). The elaboration of SPH, including general issues such as derivations of the particle form of continuity equation, momentum equation/ thermal equation, construction of artificial viscosity, a nearest particle searching algorithm, and improved approaches, as well as various applications, were presented in both the early and the recent works (for example, Schussler and Schmitt, 1981; Monaghan, 1982, 1985, 1987, 1988, 1989, 1992, 1994, 1996, 2001, 2002, 2012; Monaghan and Gingold, 1983; Monaghan and Poinracic, 1985; Monaghan and Lattanzio, 1985; Monaghan and Kocharyan, 1995; Stellingwerf and Wingate, 1994; Monaghan and Kos, 1999; Capuzzo-Dolcetta and Robert, 2000; Pozorski and Wawrenczuk, 2002; Liu, 2003; Zhang and Batra, 2004; Liu and Liu, 2003; Quinlan et al., 2006; Oger et al., 2006, 2007; Crespo et al., 2007; Jiang et al., 2007; Basa et al., 2009; Dalrymple and Knio, 2010). Since the collective movement of those particles is similar to the movement of a fluid or gas flow, it is identified by the name “smoothed particle hydrodynamics.” Applications of the SPH method are shortly summarized as follows. Application on incompressible fluid SPH was first applied on incompressible fluids with the free surface problem by Monaghan (1994) and then for incompressible flows of low Reynolds number by Morris et al. (1997) and further for incompressible viscous flows by Pozorski and Wawrenczuk (2002). In those early works, fluid such as water was assumed to be slightly compressible. Therefore the quasiincompressible equation of state was applied to calculate pressure. Time-stepping size depended on sound speed, which was adjusted to confine the fluid density variation. This weakly compressible SPH (WCSPH) method has proven able to successfully simulate Poiseuille flow by comparing the results obtained from the FVM (Lobovsky´ and Vimmr, 2007). However, WCSPH requires a very small time step, and a small density error may cause significant nonphysical pressure fluctuation (Lee et al., 2008, 2010). An approximated pressure projection method was developed (Cummins and Rudman, 1999) to enforce fluid incompressibility by solving a pressure Poisson’s equation. Afterward, a truly incompressible SPH (ISPH) method was proposed (Shao and Lo, 2003; Shao, 2009) using prediction
correction fractional steps to upgrade the related physical properties. The temporal velocity field was integrated forward in time without considering the pressure effect in the first step. The temporal density obtained from the first step was then implicitly projected onto a velocity divergence
free space to satisfy incompressibility in the second

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

491

step. The pressure values were calculated through a Poisson’s equation matrix. This ISPH was widely used since then (Hosseini et al., 2007). However, a more straightforward way to ensure incompressibility is to use constant density (Lee et al., 2010). Application on multiphase flow For more realistic simulations of fluid motions, multiphase flow needs to be considered. The early SPH application on multiphase flow is about compressible fluid like dusty gas. The mixed fluid can be treated as a new type of fluid, and void fraction is used (Johnson and Beissel, 1996; Monaghan and Kocharyan, 1995). The mass density of this new fluid can be updated based on the continuity equation. Since the density of dust is known, the gas density can be obtained through the void fraction. It is more difficult to consider the situation when air is mixed into water because the large density ratio may cause instability in the algorithm. One needs to be careful to choose an appropriate particle evolution form in order to ensure that the density is continuous for each phase. If WCSPH is used, the densities may need to be reinitialized during the computation process (Colagrossi and Landrini, 2003; Colagrossi et al., 2009). More methods have been developed to deal with multiphase flow. In the latest work about compressible multiphase flows, the pressure of gas was assumed to be constant, and the density was calculated from the equation of state (Ritchie and Thomas, 2001). Different forms of particle evolution were developed to get rid of the influence of a large density ratio (Grenier and Touze, 2008; Grenier et al., 2009; Hu and Adams, 2006, 2007). However, most of these works are based on the WCSPH method; improvement in using ISPH for multiphase flows related to incompressible flow is still necessary. Application on solids SPH is also applied on solids. When a solid continuum is in a state of isotropic tension, the solid is stretched, and the SPH particles attract one another to resist the stretching, causing particle clumping. The use of the standard SPH equation in conjunction with an explicit time-stepping scheme leads to an unstable time integration for any time step size, and this instability cannot be eliminated by introducing artificial viscosity (Bonet and Lok, 1999; Bonet et al., 2004). This phenomenon is called tensile instability, which is a major deficiency of traditional SPH. Besides, SPH cannot produce correct estimations on the boundaries because not enough particles are nearby. Various methods were developed to improve the accuracy, stability, and consistency of SPH by adopting some of the techniques used in MGMs. Application on fluid
solid interaction problems SPH is effective in simulating large fluid motions, and it has been applied to FSI problems by combining FEM for structure simulations. An approximate master
slave algorithm has been proposed (Attawy et al., 1994; Johnson, 1994; Johnson and Beissel, 1996, Johnson et al., 1996a,b). In this algorithm, the contact constraint was satisfied by applying a contact force to both the slave node and the master surface to remove any penetration in the next time step. The force is always normal to the corresponding element surface, and sliding in the tangential direction on the interface between particles and elements is allowed in order to be consistent with velocity changes. When a slave node overlaps the master segment, the normal velocities of the involved three nodes are adjusted to conserve linear momentum, angular momentum, and a velocity matching the slave node on master segment. This approximate coupling method has an advantage in that the particle sums remain symmetric, so that the deformation gradient is computed correctly (Rabczuk et al., 2000, 2006), and the approximation is independent of the relation between the particle distance and the FE length. Another possibility to simplify coupling is to apply the SPH method on solids as well. The shear stress, deviator stress, and pressure formalism can be derived by applying SPH directly to the strain rate tensor and the rotation rate tensor. This makes the transfer of information between the fluid and structure domains easier because similar algorithms are used on both parts. It also makes the simulation more efficient in the case of large deformations of the solid. There are two coupling models for using SPH for both solids and fluids. One is to treat all the particles in the same way regardless of their nature, and a correction is applied to avoid particle penetration (Rafiee and Thiagarajan, 2008, 2009); the other is to treat the fluid and solid separately and to find the right position of FSI interface and its normal direction (Antoci et al., 2007). The force acting on the solid from the fluid is computed through the evaluation of an approximated surface integral of fluid pressure. This force on the interface follows the action
reaction principle. In some cases when the tensile instability becomes serious, an artificial viscosity and artificial stress term can be used to deal with the problem (Antoci et al., 2007; Bui et al., 2007).

492

Fluid
Solid Interaction Dynamics

The MFE
SPM provides an effective approach to simulate nonlinear FSI problems involving breaking waves and fluid
solid separation (FSS) (Sun, 2013, 2016; Sun et al., 2011, 2012, 2013, 2014, 2015a,b, 2016, 2017; Javed et al., 2013a,b, 2014a,b, 2016; Javed, 2015; Wu et al., 2014). In this approach, the solid motion includes large rigid motion with small elastic deformation represented by mode summation, while the fluid is formulated by SPM. On the FSI interface, the velocity consistency and force balance conditions must be satisfied, and at each time step the solid current position is the boundary of SPH particles, and large free surface fluid motion is traced (Monaghan, 1994). The new developments are as follows: (1) to coordinate with the projection method for incompressible fluids, the force balance condition is expressed in pressure gradient form using the acceleration at the FE nodes of the solid, which can avoid the calculation of the solid stress; (2) a mixed source term in the pressure Poisson equation with no artificial term in the formulation; (3) a new version of “cell-link” neighbor particle searching strategy, which reduces about 6.5/9 (B72%) of the searching area compared with traditional the “cell-linked” algorithm (Sun, 2016; Sun et al., 2014, 2015a,b, 2016, 2017); and (4) the hybrid SP
Eulerian grid is used to deal with FSI problems involving viscous flows in large domains, in which the active domain around the FSI interface adopts an SP scheme while the outside domain uses the Eulerian grid to improve the calculation efficiency (Javed, 2015; Javed et al., 2013a,b, 2014a,b). The developed method has been validated by experiments or other publications. In this chapter, we will present the detailed theory, numerical process, and applications of this method.

12.2

Fundamentals of smoothed particle hydrodynamics

For readers to more clearly to understand MFE
SPM for nonlinear FSI problems, in this subsection, we briefly present the fundamental formulations and techniques of SPH, which has been successful in FSI analysis. For detailed information on this method, see Liu and Liu (2003) and Li and Liu (2004).

12.2.1

Smoothed particle hydrodynamics interpolations

12.2.1.1 Integral representation of a function The concept of integral representation of a function used in the SPH method starts with the following identity: ð f ðxÞ 5 f ðx0 Þδðx 2 x0 Þdx0 ;

(12.3)

Ω

where x is a 3D position vector in an integration domain Ω in 3D space, and δðx 2 x0 Þ is the Dirac delta function defined by  ð N; x0 -x; δðx 2 x0 Þ 5 δðx 2 x0 Þdx0 5 1; xAΩ: (12.4) 0; x0 6¼ x; Ω Eq. (12.3) implies that a function can be represented in an integral form. Since the Dirac delta function is used, the integral representation in Eq. (12.1) is exact and rigorous, as long as f(x) is defined and continuous in Ω. If the Delta function kernel δðx 2 x0 Þ is replaced by a smoothing function Wðx 2 x0 ; hÞ, the integral representation of f(x) is given by Eq. (12.2), that is, ð f ðxÞ  f ðx0 ÞWðx 2 x0 ; hÞdx0 ; (12.5) Ω

0

where Wðx 2 x ; hÞ is called the smoothing kernel function, or kernel function. In the smoothing function Wðx 2 x0 ; hÞ, h is the smoothing length defining the influence area of the smoothing function. Since the function Wðx 2 x0 ; hÞ is not the Dirac function, the integral representation in Eq. (12.5) can only be an approximation. This is the origination of the term of kernel approximation. The kernel approximation operator is marked by angle brackets (h i), and therefore Eq. (12.5) is rewritten as ð   f ðxÞ 5 f ðx0 ÞWðx 2 x0 ; hÞdx0 : (12.6) Ω

The smoothing function Wðx 2 x0 ; hÞ is usually chosen to satisfy the following conditions: 1. Unity: The normalization or unity condition states ð Wðx 2 x0 ; hÞdx0 5 1: Ω

This condition implies that the chosen kernel function can exactly represent the unity function f ðxÞ  1:

(12.7)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

493

2. Delta function property: Requires that when the smoothing length approaches zero lim Wðx 2 x0 ; hÞ 5 δðx 2 x0 Þ:

h-0

(12.8)

3. Compact support, positivity, and decay: Wðx 2 x0 ; hÞ 5 0; Wðx 2 x0 ; hÞ $ 0;

when jx 2 x0 j . κh; when jx 2 x0 j # κh;

(12.9)

where κ is a constant related to the smoothing function for point at x, and it defines the effective (nonzero) area of the smoothing function. This effective area is called the support domain for the smoothing function of f(x) at x. Based on this compact condition, integration over the entire problem domain is localized to integration over the support domain. Therefore the integration domain Ω can be the same as the support domain. The value of the smoothing function for a particle should be monotonically decreasing with the increase of the distance away from the particle, decay. 4. Symmetric property and smoothness: The smoothing function should be an even function and sufficiently smooth, that is, its derivatives should be sufficiently high with respect to the spatial coordinates. To estimate the accuracy of the SPH integration approximation, we use the Taylor series expansion of f(x0 ), assumed differentiable, around x. From Eq. (12.6), it follows that 1 f vðxÞðx0 2xÞ2 1 ? 2   k  n X ð21Þk hk f ðkÞ ðxÞ x2x0 x 2 x0 1 rn 5 ; k! h h k50

f ðx0 Þ 5 f ðxÞ 1 f 0 ðxÞðx0 2 xÞ 1

from which, when Eq. (12.6) is substituted, it follows   f ðxÞ  f ðxÞ   k  ð X n ð21Þk hk f ðkÞ ðxÞ x2x0 x 2 x0 0 0 Wðx 2 x ; hÞdx 1 rn ; 5 k! h h Ω k50   f ðxÞ  f ðxÞ    ð  n X ð21Þk hk f ðkÞ ðxÞ x2x0 k x 2 x0 Wðx 2 x0 ; hÞdx0 1 rn 5 k! h h Ω k50   n X x 2 x0 Ak f ðkÞ ðxÞ 1 rn 5 ; h k50 k ð  ð21Þk hk x2x0 Wðx 2 x0 ; hÞdx0 : Ak 5 k! h Ω

(12.10a)

(12.10b)

(12.10c)

Comparing the left- and right-hand sides of Eqs. (12.10a)
(12.10c), for the function f(x), in order to be approximated to the nth order, the coefficients Ak must be equal to the counterparts for f (k)(x) on the left-hand side. Therefore the following conditions for the smoothing function W are obtained: ð A0 5 M 0 ; M0 5 Wðx 2 x0 ; hÞdx0 5 1; ðΩ M1 5 ðx 2 x0 ÞWðx 2 x0 ; hÞdx0 5 0; A1 5 2 M 1 ; ðΩ M2 (12.11) ; M2 5 ðx2x0 Þ2 Wðx 2 x0 ; hÞdx0 5 0; A2 5 2! Ω ^ ð21Þn Mn An 5 ; n!

ð Mn 5

Ω

ðx2x0 Þn Wðx 2 x0 ; hÞdx0 5 0:

Here Mk is called the kth moments of the smoothing function. The first equation in Eq. (12.11) is the unity condition expressed in Eq. (12.7), and the second represents the symmetric property. Satisfaction of these two conditions ensures

494

Fluid
Solid Interaction Dynamics

first-order consistency for the SPH kernel approximation of a function. This implies that the kernel approximation is of second-order accuracy if the normalization and even function conditions are valid.

12.2.1.2 Integral representation of the derivatives of a function First derivative The approximation of the spatial derivative f;i ðxÞ; i 5 1; 2; 3; of the function f(x) is obtained simply by replacing f(x) in Eq. (12.6) by f,i(x), that is, ð   (12.12) f;i ðxÞ 5 f;i ðx0 ÞWðx 2 x0 ; hÞdx0 : Ω

Here, the derivative is operated with respect to the primed coordinate. Using the rule for the derivative of a product of two functions, we obtain f;i ðx0 ÞWðx 2 x0 ; hÞ 5 ½ f ðx0 ÞWðx2x0 ; hÞ;i 2 f ðx0 ÞW;i ðx 2 x0 ; hÞ;

(12.13)

from which Eq. (12.12) is reduced to ð ð   f;i ðxÞ 5 ½ f ðx0 ÞWðx2x0 ; hÞ;i dx0 2 f ðx0 ÞW;i ðx 2 x0 ; hÞdx0 :

(12.14)

Ω

Ω

The first integration on the right-hand side of Eq. (12.14) can be converted into an integral over the surface S of the integration domain Ω by using the Green theorem, so that we derive ð   f;i ðxÞ 5 IS 2 f ðx0 ÞW;i ðx 2 x0 ; hÞdx0 ; Ω ð (12.15) 0 0 0 IS 5 f ðx ÞWðx 2 x ; hÞν i dx ; S

where vi denotes the unit vector along the outer normal to the surface S. Substituting the Taylor expansion given in Eqs. (12.10a)
(12.10c) into the second integral of Eq. (12.15), we obtain   f;i ðxÞ  f;i ðxÞ   k )  ð (X n ð21Þk hk f ðkÞ ðxÞ x2x0 x 2 x0 0 0 5 IS 2 W;i ðx 2 x ; hÞdx 1 rn k! h h Ω k50 8 9  k  ð  n X ð21Þk hk f ðkÞ ðxÞ < x2x0 = x 2 x0 W;i ðx 2 x0 ; hÞdx0 1 rn 5 IS 2 ; k! h h (12.16) Ω: k50   n X x 2 x0 ðkÞ Að1Þ f ðxÞ 1 r 5 IS 1 ; n k h k50  ð  ð21Þk11 hk x2x0 k 5 W;i ðx 2 x0 ; hÞdx0 : Að1Þ k k! h Ω It is clear that the derivative f,i(x) can be approximated to nth order if the following conditions are satisfied: ð ð1Þ M0 5 W;i ðx 2 x0 ; hÞdx0 5 0; ðΩ ð1Þ M1 5 ðx 2 x0 ÞW;i ðx 2 x0 ; hÞdx0 5 1; ðΩ (12.17) ð1Þ M2 5 ðx2x0 Þ2 W;i ðx 2 x0 ; hÞdx0 5 0; Ω

^ Mnð1Þ and

5

ð Ω

ðx2x0 Þn W;i ðx 2 x0 ; hÞdx0 5 0;

 Wðx2x0 ; hÞS  0:

(12.18)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

495

The condition in Eq. (12.18) is in fact the compact support condition. As shown in Fig. 12.1, the supporting domain of the smoothing function W is located within the problem domain, for which the surface integration IS in Eq. (12.15) vanishes. However, as shown in Fig. 12.2, when the supporting domain of the smoothing function W intersects with the problem domain, there exists a truncated boundary on which the value of the smoothing function W is not zero. Therefore the surface integration IS in Eq. (12.15) is no longer zero. For convenience in numerical simulation, this surface integration can be treated as zero, so that Eq. (12.15) is simplified as ð   f;i ðxÞ 5 2 f ðx0 ÞW;i ðx 2 x0 ; hÞdx0 : (12.19) Ω

Some modifications should also be made to remedy the boundary effects. Eq. (12.19) reveals a very important characteristic of the SPH method, which allows the spatial derivatives of a field function to be determined by the derivatives of the smoothing function W rather than by the field function itself. Second derivative The approximation of the second derivative can be obtained by substituting f(x) in Eq. (12.6) with its second derivative f;ij ðxÞ, which gives ð   f;ij ðxÞ  f;ij ðxÞ 5 f ;ij ðx0 ÞWðx 2 x0 ; hÞdx0 : (12.20) Ω

Using integration by parts and the Green theorem, this equation is derived as follows:   f;ij ðxÞ  f;ij ðxÞ ð ð 0 0 5 f ;i ðx ÞWðx 2 x ; hÞν j dS 2 f ;i ðx0 ÞW; j ðx 2 x0 ; hÞdx0 ðS ðΩ ð 0 0 5 f ;i ðx ÞWðx 2 x ; hÞν j dS 2 f ðx0 ÞW; j ðx 2 x0 ; hÞν i dS 1 f ðx0 ÞW;ij ðx 2 x0 ; hÞdx0 : S

Ω

S

Problem domain

W

κh Support domain FIGURE 12.1 Supporting domain of the smoothing function W is located within the problem domain.

Problem domain

W Boundary

Interior

κh Support domain FIGURE 12.2 Supporting domain of the smoothing function W intersects with the problem domain, so that it is truncated by the boundary.

(12.21)

496

Fluid
Solid Interaction Dynamics

Replacing the function f(x0 ) in the third integration of Eq. (12.21) by the Taylor expansion in Eqs. (12.10a)
(12.10c), we have   k  ð ð ð X n ð21Þk hk f ðkÞ ðxÞ x2x0 x 2 x0 0 0 0 0 0 0 f;ij ðxÞ 5 f ;i ðx Þ Wðx 2 x ; hÞν j dS 2 f ðx ÞW;j ðx 2 x ; hÞν i dS 1 W;ij ðx 2 x ; hÞdx 1 rn k! h h S S Ω k50     ð ð ð k n X ð21Þk hk f ðkÞ ðxÞ x2x0 x 2 x0 5 f ;i ðx0 ÞWðx 2 x0 ; hÞν j dS 2 f ðx0 ÞW;j ðx 2 x0 ; hÞν i dS 1 W;ij ðx 2 x0 ; hÞdx0 1 rn k! h h S S Ω k50   ð ð n 0 X ð2Þ x2x 5 f ;i ðx0 ÞWðx 2 x0 ; hÞν j dS 2 f ðx0 ÞW;j ðx 2 x0 ; hÞν i dS 1 Ak f ðkÞ ðxÞ 1 rn ; h S S k50 (12.22) where Að2Þ k 5

ð21Þk hk k!

 ð  x2x0 k W;ij ðx 2 x0 ; hÞdx0 : h Ω

(12.23)

It can be realized that f;ij ðxÞ is able to be approximated to nth-order accuracy if the following conditions are satisfied: ð ð2Þ M0 5 W;ij ðx 2 x0 ; hÞdx0 5 0; ðΩ ð2Þ M1 5 ðx 2 x0 ÞW;ij ðx 2 x0 ; hÞdx0 5 0; ðΩ (12.24) M2ð2Þ 5 ðx2x0 Þ2 W;ij ðx 2 x0 ; hÞdx0 5 2; Ω

^ Mnð2Þ 5

and

ð Ω

ðx2x0 Þn W;ij ðx 2 x0 ; hÞdx0 5 0;

 Wðx2x0 ; hÞ S  0;

(12.25)

 W;i ðx2x0 ; hÞ S  0:

(12.26)

The conditions in Eqs. (12.25) and (12.26) confirm the surface integral terms in Eq. (12.22) for an arbitrary function f(x0 ) to vanish.

12.2.1.3 Rules for nth-order accuracy For SPH approximations with nth-order accuracy, we have learned that the conditions in Eqs. (12.11), (12.17), (12.18), and (12.24)
(12.26) should be satisfied for a function, its first and second derivatives, respectively. Now we further demonstrate that these conditions can be synthesized as follows. Based on the Green theorem, we have the following equations: ð ð W;i ðx 2 x0 ; hÞdx0 5 Wðx 2 x0 ; hÞν i dS 5 0; ðΩ ðS 0 0 W;ij ðx 2 x ; hÞdx 5 W;i ðx 2 x0 ; hÞν j dS 5 0; Ω S ð ð ð (12.27) 0 0 0 0 0 0 0 ðx 2 x ÞW;ij ðx 2 x ; hÞdx 5 ðx 2 x ÞW;i ðx 2 x ; hÞdS 1 W;i ðx 2 x ; hÞdx Ω Ω ðS ð 0 0 5 ðx 2 x ÞW;i ðx 2 x ; hÞν j dS 1 Wðx 2 x0 ; hÞν i dx0 : S

S

The first equation in Eq. (12.27) confirms that the first condition in Eq. (12.17) is actually the result of Eq. (12.18) or Eq. (12.25), and the second and third equations in Eq. (12.27) reveal that the first and second conditions in Eq. (12.24) are actually determined by Eqs. (12.25) and (12.26).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

497

Now we can prove that the rest of conditions in Eqs. (12.17) and (12.24) can be derived from Eqs. (12.11), (12.25), and (12.26). Actually, by using integration by parts and the Green theorem, we have ð ðx2x0 Þk Wðx 2 x0 ; hÞdx0 Ω ð 1 52 ½ðx2x0 Þk11 ;i Wðx 2 x0 ; hÞdx0 k11 Ω ð ð 21 1 0 k11 0 5 ðx2x Þ Wðx 2 x ; hÞν i dS 1 ðx2x0 Þk11 W;i ðx 2 x0 ; hÞdx0 k11 S k11 Ω ð 1 5 ðx2x0 Þk11 W;i ðx 2 x0 ; hÞdx0 (12.28) k11 Ω ð 21 5 ðx2x0 Þk12 W;i ðx 2 x0 ; hÞdx0 ðk 1 1Þðk 1 2Þ Ω ð ð 21 1 5 ðx2x0 Þk12 W;i ðx 2 x0 ; hÞν j dx0 1 ðx2x0 Þk12 W;ij ðx 2 x0 ; hÞdx0 ðk 1 1Þðk 1 2Þ Ω ðk 1 1Þðk 1 2Þ Ω ð 1 5 ðx2x0 Þk12 W;ij ðx 2 x0 ; hÞdx0 : ðk 1 1Þðk 1 2Þ Ω Therefore for each value of k 5 0; 1; 2; . . .; n, the condition in Eq. (12.11) gives, respectively, the corresponding condition with k 1 1 in Eq. (12.17), as well as the corresponding one with k 1 2 in Eq. (12.24). To summarize this demonstration, we confirm that if a function and its first two derivatives are to be represented to nth-order accuracy, then the smoothing function should satisfy the rules in Eqs. (12.11), (12.25), and (12.26), which are rewritten as follows: ð M0 5 Wðx 2 x0 ; hÞdx0 5 1; ðΩ M1 5 ðx 2 x0 ÞWðx 2 x0 ; hÞdx0 5 0; ðΩ (12.29) M2 5 ðx2x0 Þ2 Wðx 2 x0 ; hÞdx0 5 0; Ω

^

ð

Mn 5 and

Ω

ðx2x0 Þn Wðx 2 x0 ; hÞdx0 5 0;

 Wðx2x0 ; hÞ S  0;  W;i ðx2x0 ; hÞ   0:

(12.30)

S

These rules can be used to construct the smoothing functions. The rules in Eq. (12.29) imply that the smoothing function is able to reproduce polynomials, and the conditions in Eq. (12.30) require the property of compact support for the smoothing function and its first derivative. Satisfying these two conditions, the smoothing function can represent a function and its first two-order derivatives to the nth-order accuracy. Following the same procedure, the rules for producing higher-order derivatives of a function can also be derived. Except for the rule in Eq. (12.29), in general, the compact support conditions for the mth-order derivative to be approximated to nth-order accuracy are @ðαÞ Wðx2x0 ; hÞ   0; @x0i α1 @x0j α2 @x0k α3 S

α 5 α1 1 α2 1 α3 ;

α 5 0; 1; 2; . . .; m 2 1:

(12.31)

12.2.1.4 Consistency of kernel approximation The conditions in Eqs. (12.29) and (12.30) are derived using the Taylor series expansion, which is directly related to the consistency concept for the traditional FDM. In FEM analysis, to ensure convergence, the shape function must

498

Fluid
Solid Interaction Dynamics

satisfy a certain degree of consistency. For a displacement model, its shape function must, at least, be able to exactly represent a rigid displacement and a constant strain, that is, a linear displacement field of the body. For a higher degree of consistency, the FEM shape function is required to model the corresponding polynomial with the higher order. Referring to the consistency concept from FEM, we investigate whether the kernel approximation can represent polynomials as follows. For a constant f(x) 5 c0, from the first condition in Eq. (12.29), it follows that ð c0 Wðx 2 x0 ; hÞdx0 5 c0 ; (12.32) Ω

implying the zeroth-order consistency of the kernel approximation.For a linear function f ðxÞ 5 c0 1 c1 x, from the first and second conditions in Eq. (12.29), it follows that ð ð xWðx 2 x0 ; hÞdx0 5 x 5 x0 Wðx 2 x0 ; hÞdx0 ; (12.33) Ω

so that

Ω

ð Ω

ðc0 1 c1 x0 ÞWðx 2 x0 ; hÞdx0 5 c0 1 c1 x;

which confirms the first-order consistency of the kernel approximation. Furthermore, from the first three conditions in Eq. (12.29) and the result in Eq. (12.33), we have ð ðx2x0 Þ2 Wðx 2 x0 ; hÞdx0 Ω ð 5 ðx2 2 2xx0 1 x02 ÞWðx 2 x0 ; hÞdx0 Ω ð ð ð 2 5 x Wðx 2 x0 ; hÞdx0 2 2x x0 Wðx 2 x0 ; hÞdx0 1 x02 Wðx 2 x0 ; hÞdx0 Ω Ω Ω ð 2 2 02 0 0 5 x 2 2x 1 x Wðx 2 x ; hÞdx 5 0;

(12.34)

(12.35)

Ω

which implies

ð x2 5

Ω

x02 Wðx 2 x0 ; hÞdx0 :

Therefore for a quadratic function f ðxÞ 5 c0 1 c1 x 1 c2 x2 , we can obtain ð ðc0 1 c1 x0 1 c2 x02 ÞWðx 2 x0 ; hÞdx0 5 c0 1 c1 x 1 c2 x2 ;

(12.36)

(12.37)

Ω

which confirms the second-order consistency of the kernel approximation. Generally, using the binomial theorem and the following related results: ðx2x0 Þk 5

k X ð21Þm Ckm xk2m x0m ; m50

0 5 ð121Þk 5 k21 X m50

k X

ð21Þm Ckm ;

m50

ð21Þm Ckm 5 ð21Þk11 Ckk 5 ð21Þk11 ;

(12.38)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

499

we can calculate the integration of the kth moment in Eq. (12.29): ð ðx2x0 Þk Wðx 2 x0 ; hÞdx0 Ω ð X k 5 ð21Þm Ckm xk2m x0m Wðx 2 x0 ; hÞdx0 Ω m50

k X

ð

x0m Wðx 2 x0 ; hÞdx0 m50 ð ð ð k 0 0 1 k21 01 0 0 2 k22 5 x Wðx 2 x ; hÞdx 2 Ck x x Wðx 2 x ; hÞdx 1 Ck x x02 Wðx 2 x0 ; hÞdx0 1 ? Ω Ω Ω ð 1 ð21Þk Ckk x0k Wðx 2 x0 ; hÞdx0 5 0; 5

ð21Þm Ckm xk2m

(12.39a)

Ω

Ω

where the combination  Ckm 5

k m

 5

k! kðk 2 1Þ?ðk 2 m 1 1Þ 5 ; ðk 2 mÞ!m! m!

m # k:

(12.39b)

For each value of k 5 0, 1, 2, Eq. (12.33) gives equations in Eqs. (12.32), (12.33), and (12.36), respectively. Following the same argument, if the conditions in Eq. (12.29) are satisfied, we obtain xn

n21 X

ð ð21Þm Ckm 5 ð21Þn11 xn 5 ð21Þn11

m50

Ω

x0n Wðx 2 x0 ; hÞdx0 ;

k 5 n:

(12.40)

As a result of this, an nth-order polynomial can be exactly represented by the kernel function, and the nth-order consistency is reached.

12.2.2

Particle approximation

12.2.2.1 Approximation of a function In the SPH method, the investigated system is represented by a finite number of particles, each of which carries an individual mass and occupies an individual space, which is achieved by the particle approximation process as follows. Fig. 12.3 shows the supporting domain, a circle of radius of κh, of the smoothing function W for particle I. This circle is filled by N small particles, of which each particle, such as particle J located at position x0 J , carries its mass mJ and occupies an infinitesimal volume dx0 5 ΔVJ . Introducing the mass density ρJ , (J 5 1, 2, . . ., N), we can denote the particle mass in the form mJ 5 ΔVJ ρJ :

(12.41)

The integral representation of function f(x) can be written in the following summation: ð   f ðxÞ 5 f ðx0 ÞWðx 2 x0 ; hÞdx0 Ω

5

N X

f ðxJ ÞWðx 2 xJ ; hÞΔVJ

J51

5

N X mJ J51

ρJ

f ðxJ ÞWðx 2 xJ ; hÞ:

(12.42)

500

Fluid
Solid Interaction Dynamics

S J

W

rIJ

I κh Ω

FIGURE 12.3 Particle approximations based on the particles within the supporting domain, the circle of radius of κh, of the smoothing function W for particle I.

The particle approximation of a function at particle I can finally be represented as 

N  X mJ f ðxJ ÞWIJ ; f ðxI Þ 5 ρ J51 J

(12.43)

where WIJ 5 WðxI 2 xJ ; hÞ:

(12.44)

Physically, Eq. (12.43) represents that the value of a function at particle I is approximated by the average of those values of the function at all the particles in the supporting domain of particle I weighted by the smoothing function. Following the same derivation, the particle approximation for the spatial derivative of the function is N X   mJ f;i ðxÞ 5 2 f ðxJ ÞW;i ðx 2 xJ ; hÞ; ρ J51 J

(12.45)

where W;i ðx 2 xJ ; hÞ 5

@W @W 52 : @xJi @xi

(12.46)

The particle approximation for the spatial derivative of a function at particle I can finally be written in the form * + N X   @f ðxI Þ mJ @WðxI 2 xJ ; hÞ f ðxJ Þ f;i ðxI Þ 5 5 @xIi @xIi ρ J51 J (12.47) N X mJ @WIJ 5 f ðxJ Þ ; ρ @xIi J51 J where @WIJ xIJ @WIJ 5 i ; @xIi rIJ @rIJ rIJ 5 ½ðxIi 2xJi ÞðxIi 2xJi Þ

(12.48) 1=2

;

xIJ i

5 xIi 2 xJi :

Eq. (12.48) represents that the value of the spatial derivative of a function at particle I is approximated by the average of those values of the function at all the particles in the supporting domain of particle I weighted by the derivative of the smoothing function. It has been seen that the particle approximation in Eqs. (12.43) and (12.47) actually converts the continuous integral representations of a function and its derivatives to the discretized summations based on an arbitrary set of particles. This is a key approximation, making SPH method very simple without using a background mesh for numerical integration.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

Substituting the function f(x) with the density function ρ in Eq. (12.43), we obtain N X ρI 5 mJ WIJ ;

501

(12.49)

J51

which is generally referred to as the summation density approach. Eq. (12.49) states that the density of a particle is a weighted average of all the particles in its support domain.

12.2.2.2 Techniques in deriving smoothed particle hydrodynamics formulations Based on the preceding procedure of kernel and particle approximations, we can always convert PDEs into their corresponding SPH formulations. A number of ways are used. One approach is to multiply each term of the PDEs with the smoothing function and then to integrate over the volume by means of integration by parts, the Green theorem, and Taylor expansions. A more straightforward approach is directly using Eqs. (12.43) and (12.47). In this approach, the following two identifies are often employed to place the density inside the spatial derivative operator:

1 ðρf ðxÞÞ;i 2 ρ;i f ðxÞ ; ρ   f ðxÞ f ðxÞ f;i ðxÞ 5 ρ 1 2 ρ;i : ρ ;i ρ

f;i ðxÞ 5

Using Eq. (12.47), we obtain that

(12.50)

N   X @WIJ ðρf ðxÞÞ;i 5 mJ f ðxJ Þ ; @xIi J51 N   X @WIJ ρ;i 5 mJ ; @xIi J51  X  N f ðxÞ f ðxJ Þ @WIJ 5 mJ 2 ; ρ ;i J51 ρJ @xIi

(12.51)

from which, when substituted into Eq. (12.50), it follows that the derivative of function f(x) at particle I is obtained as 8 9 N 21
12.2.2.3 Consistency of the particle approximation The consistency concept for the kernel approximation is based on the continuous form of the integral representation, which cannot ensure consistency for the discrete form of the particle approximation. This phenomenon is called particle inconsistency. Consider the first and second conditions in Eq. (12.29); we have their discrete forms N X J51

Wðx 2 x0J ; hÞΔx0J 5 1;

N X ðx 2 x0J ÞWðx 2 x0J ; hÞΔx0J 5 0; J51

(12.54)

502

Fluid
Solid Interaction Dynamics

(A)

(B)

W

W

Particle Boundar x

x

FIGURE 12.4 SPH approximation for 1D case: (A) particles in the support domain are of irregular distribution; (B) support domain is truncated by the boundary. SPH, Smoothed particle hydrodynamic; 1D, one-dimensional.

where N is the total number of particles in the support domain for the given particle located at x. Since the particle distribution in the support domain is irregular, as shown in Fig. 12.4A, being nonsymmetrical about point x, or truncated by the boundary, as shown in Fig. 12.4B, the summations in Eq. (12.54) are not necessarily satisfied. One of proposed approaches to obtain a kth-order consistency in the discrete form is as follows. We may assume the smoothing function in the form k

x2x I X J Wðx 2 xJ ; hÞ 5 AI ðx; hÞ ; (12.55) h I50 from which, the discretized form of Eq. (12.29) can be written as  I k X x2xJ AI ðx; hÞ Wðx 2 xJ ; hÞ 5 h I50 2 3  I I N k k N  X X X X x2xJ 4 AI ðx; hÞ x2xJ 5ΔxJ 5 AI ðx; hÞ ΔxJ 5 1; h h J51 I50 I50 J51 2 3   X I I11 N  k k N  X X X x 2 xJ 4 x2xJ 5 x2xJ ΔxJ 5 AI ðx; hÞ AI ðx; hÞ ΔxJ 5 0; h h h J51 I50 I50 J51

(12.56)

^ 2 3   k k I I1k N k N  X X X x2xJ 4X x2xJ 5 x2xJ ΔxJ 5 AI ðx; hÞ AI ðx; hÞ ΔxJ 5 0: h h h J51 I50 I50 J51 Eq. (12.56) is a matrix equation to solve k 1 1 coefficients AI ðx; hÞ; that is, 2

H0 6H 6 1 6 4 ^ Hk

H1 H2 ^ Hk11

32 3 2 3 A0 1 Hk 6A 7 607 H11k 7 76 1 7 6 7 76 7 5 6 7; & ^ 54 ^ 5 4 ^ 5 0 Ak ? Hk1k ? ?

where HI ðx; hÞ 5

N X x2xJ I J51

h

ΔxJ :

(12.57)

(12.58)

In practice, this approach shows that restoring the consistency in discrete form leads to some simulation problems. For example, the resultant smoothing function is negative in some parts of the region, leading to the unphysical representation of field variables, or it may not be symmetric, nonmonotonically decreasing with the increase of the particle distance. Therefore in solving practical problems, special care must be taken (Liu and Liu 2003).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

12.2.3

503

Construction of smoothing functions

12.2.3.1 Polynomial form We assume a polynomial form as Wðx 2 x0 ; hÞ 5 WðRÞ 5

n X

AI R I ;

I50 0

0 1=2

R 5 ½ðx2x ÞUðx2x Þ 5 ½ðxi 2x0i Þðxi 2x0i Þ1=2 ; @R 2 ðxi 2 x0i Þ ; R;i 5 0 5 @xi R

(12.59a)

δij R2 2 ðxi 2 x0i Þðxj 2 x0j Þ @ R R;ij 5 0 0 5 ; @xi @xj R3 2

R;ii 5

3R2 2 R2 2 5 ; R R3

of which the first and second derivatives with respect to x0i are given by W;i ðRÞ 5

dW R;i ; dR

W;ij ðRÞ 5

d2 W dW R;ij ; R;i R;j 1 dR2 dR

d2 W ðxi 2 x0i Þðxi 2 x0i Þ 2dW d2 W 2dW 5 : 1 1 W;ii 5 dR2 R2 RdR dR2 RdR

(12.59b)

The condition in Eq. (12.30) requires the first derivative of the smoothing function, for it exists at R 5 0. We should have  dW 2 A1 ; ðxi 2 x0i Þ-01 ; ð0Þ lim R 5 W;i ð0Þ 5 (12.59c) ;i 0 A1 ; ðxi 2 x0i Þ-02 ; xi -xi dR so that A1 must be zero, and the polynomial in Eq. (12.59a) takes the form WðRÞ 5 A0 1 A2 R2 1 ? 1 An Rn ;

R5

r : h

Example I: Quadratic smoothing function ðκ 5 1Þ Taking n 5 2 in Eq. (12.59a), we have a generalized quadratic smoothing function

r 2 WðR; hÞ 5 A0 1 A2 R2 5 A0 1 A2 : h

(12.59d)

(12.60a)

For one-dimensional (1D), two-dimensional (2D), and 3D cases, the unity condition in Eq. (12.29), respectively, yields   ð1   A 0 1 A2 2 h A0 1 A2 R2 dR 5 2h 5 1; (12.60b) 3 0   ðhh ð1  2 i

A0 A2 1 A0 1 A2 hr 2πrdr 5 h2 A0 1 A2 R2 2πRdR 5 2πh2 5 1; (12.60c) 2 4 0 0 ðh 0

A0 1 A2

 r 2 h

4πr 2 dr 5

 

A0 A2 1 h3 A0 1 A2 R2 4πR2 dR 5 4πh3 5 1: 3 5 0

ð1

(12.60d)

From the compact condition in Eq. (12.30), we require WðR 5 1; hÞ 5 A0 1 A2 5 0;

(12.60e)

504

Fluid
Solid Interaction Dynamics

2 dW/dR 1 W (R)

R

–1

1 –1 –2

FIGURE 12.5 Smoothing function WðRÞ 5 1 2 R2 and its derivative dW=dR 5 2 2R.

which, when respectively combined with Eqs. (12.60b)
(12.60d), gives the smoothing function in the quadratic form 8 1D; > < 3=ð4hÞ; 2 2 (12.60f) WðR; hÞ 5 A0 ð1 2 R Þ; A0 5 2=ðπh Þ; 2D; > : 3 15=ð8πh Þ; 3D: Fig. 12.5 shows the curves of the values of the smoothing function W and its derivative dW/dR affected by R. Example II: Quartic smoothing function ðκ 5 1Þ This function has the form WðRÞ 5 A0 1 A2 R2 1 A3 R3 1 A4 R4 ;

R5

r : h

(12.61a)

Using the condition in Eq. (12.31), if the function and its first two derivatives have compact supports, we have Wð1Þ 5 A0 1 A2 1 A3 1 A4 5 0; W 0 ð1Þ 5 2A2 1 3A3 1 4A4 5 0;

(12.61b)

Wvð1Þ 5 2A2 1 6A3 1 12A4 5 0: Using the unity condition in Eq. (12.29), we obtain   ð1 A0 1 A2 A3 A4 2 3 4 1 1 5 1; 1D: 2 hðA0 1 A2 R 1 A3 R 1 A4 R ÞdR 5 2h 3 4 5 0  

A2 A3 A4 2 A0 1 1 1 h A0 1 A2 R 1 A3 R 1 A4 R 2πRdR 5 2πh 5 1; 2 4 5 6 0

ð1 2D:



2



2

2

ð1 0

4





h A0 1 A2 R 1 A3 R 1 A4 R 4πR dR 5 4πh 3

3D:

3

3

4

2

3

A0 A2 A3 A4 1 1 1 3 5 6 7

(12.61c)

(12.61d)

 5 1:

Solving Eq. (12.61b) with one of Eqs. (12.61c)
(12.61e), we obtain the quartic function in the form 8 > < 5=ð4hÞ; 1D WðRÞ 5 A0 ð1 2 6R2 1 8R3 2 3R4 Þ; A0 5 5=ðπh2 Þ; 2D > : 105=ð16h3 Þ; 3D:

(12.61e)

(12.61f)

The performance of SPH, in terms of accuracy, stability, and computational efficiency, is influenced by the choice of kernel functions. Aluru (1999) presented a reproducing kernel particle method for meshless analysis to improve the performance of kernel functions. In the PhD thesis, Sun (2013) presented the numerical investigations of the nine commonly used kernel functions (Liu et al., 2003) in a 1D case. These nine kernel functions are as follows.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

505

12.2.3.2 Kernel functions Quadratic (Hicks and Liebrock, 2000) 

r 2  3 12 W ðr; hÞ 5 ; 4h h

0#

r # 1: h

(12.62a)

This quadratic smoothing function was used in the grid-free finite integration method. The main advantage of this kernel function is the simplicity and ease of computation, whereas the drawback is that the first derivative is not zero on the boundary of the support domain, which means that it does not have compact support for its first derivative.

Quartic (Lucy, 1977) 

r 2

r 3

r 4  5 W ðr; hÞ 5 126 18 23 ; 4h h h h

0#

r #1 h

(12.62b)

This quartic smoothing function and its first two derivatives satisfy the compact support condition.

Johnson’s quadratic (Johnson and Beissel, 1996) W ðr; hÞ 5

  1 3 r 2 3 r 3 1 2 ; h 16 h 4h 4

0#

r #2 h

(12.62c)

The specialty of this kernel is that the first derivative increases as the particles move closer, and it decreases as they move apart. This is advantageous for adjusting the position of particles to maintain stability. However, the derivative of this kernel function is not smooth at r 5 0.

Gaussian (Gingold and Monaghan, 1977) W ðr; hÞ 5

1 ðπh2 Þn=2

e2ðr=hÞ ;

2 1 W ðr; hÞ 5 pffiffiffi e2ðr=hÞ ; πh

2

n: dimension of space; n 5 1;

0#

r #2 h

(12.62d)

Gaussian kernel was used to simulate nonspherical stars originally. It is sufficiently smooth even for the secondorder derivative. However, it is not really compact as theoretically it never goes to zero. This can result in a large support domain with an inclusion of many particles for the particle approximation.

Super Gaussian (Monaghan and Poinracic, 1985a)  3 r 2 2ðr=hÞ2 W ðr; hÞ 5 2 ; n: dimension of space; e h ðπh2 Þn=2 2   1 3 r 2 2ðr=hÞ2 r 2 W ðr; hÞ 5 pffiffiffi ; n 5 1; 0 # # 2 e h h πh 2 1



(12.62e)

This is one of the higher-order smoothing functions that are devised from lower-order forms. Its main disadvantage is that the kernel is negative in some region of its support domain, which may lead to unphysical results for hydrodynamics problems.

506

Fluid
Solid Interaction Dynamics

Cubic spline (Monaghan and Poinracic, 1985a) 8 3 r 2 3 r 3 > > 1 2 1 > > 2 h 4 h > > C< r 3 W ðr; hÞ 5 n 1 22 h > > > 4 h > > > : 0

0,r,h h # r , 2h

;

(12.62f)

r $ 2h

where n denotes the dimension of space, and C is the normalization factor, of which the values 2/3, 10/7π, and 1/π are, respectively, for 1D, 2D (circular), and 3D (spherical) cases. The cubic spline function is the most widely used smoothing functions since it resembles a Gaussian function while having compact support. However, the second derivative of the cubic spline is a piecewise linear function, so the stability properties can be inferior to those of smoother kernels. Quartic spline 0 0 11 8   

> 4 4 4 > 1 @ r r r > > 12:5 2 @5 11:5 1 10 10:5 AA; > > 24h h h h > > > > > 0 1 > > >

4 4 < 1 r r @ 2:52 A 2 5 1:52 W ðr; hÞ 5 24h h h > > > > > > 1 r 4 > > > 2:52 > > 24h h > > > : 0

0 , r , 0:5h

0:5h # r , 1:5h 1:5h #

(12.62g)

r , 2:5h h

r $ 2:5h

Higher-order splines were introduced because they are better approximations of the Gaussian smoothing kernel and more stable. Quintic spline

0 1 8

r 5

r 5

r 5 > > 1 > @ 32 A > 2 6 22 1 15 12 > > 120h h h h > > > > > 0 1 > > > 1

r 5

r 5 < @ 32 A 2 6 22 W ðr; hÞ 5 120h h h > > > > > > 1 r 5 > > > 32 > > 120h h > > > : 0

0,r,h

h # r , 2h 2h #

(12.62h)

r , 3h h

r $ 3h

New quartic (Liu et al., 2003)   1 2 9 r 2 19 r 3 5 r 4 W ðr; hÞ 5 2 1 2 ; h 3 8 h 24 h 32 h

0#

r #2 h

(12.62i)

This quadratic smoothing function satisfies the compact support condition for the first derivative, and it has a smoother second derivative than the piecewise linear second derivative of the cubic function. Therefore the stability properties should be superior to those of the cubic function. However, the second derivative is not a monotonic function of r. This may lead to an incorrect approximation.

12.2.3.3 Numerical tested functions and error estimation The preceding kernel functions were used to approximate the following five common functions in a 1D case by Sun et al. (2011) and Sun (2013).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

F1 ðxÞ 5 x; F2 ðxÞ 5 x3 ; F3 ðxÞ 5 e2x ; F4 ðxÞ 5 sin x; F5 ðxÞ 5 tan x:

507

(12.63)

In the numerical tests, these functions first were approximated in integral form and then transformed into particle approximations. In each approximation, the effects of different smoothing lengths and different particle numbers were considered separately. Integration approximation As shown in Fig. 12.6, the 1D domain is uniformly divided into a set of points with distance δr, and the smoothing domain of each point is further divided into a set of points with spacing of δr0 . The integral approximation is formulated by ð N X   F ðxa Þ 5 F ðxb ÞW ðxa 2 xb Þdx0 5 F ðxb ÞWab δx0 ; (12.64) Ω b51 Wab 5 W ðxa 2 xb Þ: Here, δx0 is the spacing of points within the smoothing domain, that is, δx0 5 δr0 . The analytical result is calculated and compared with the integral approximation results, and then the error is obtained. Particle approximation In particle approximation, δx0 represents a particle volume, which is related to material properties. In order to stay consistent with the integral approximation, the volume of a particle should satisfy the following condition: mb δx0 5 ; (12.65a) ρb so the mass density can be expressed as ρ5

mb : δx0

(12.65b)

The mass of the system is assumed to be a unit as no specific material is considered. Since a uniform distribution of mass is preferred, the discretization in particle approximation is slightly different from integral approximation. The problem domain is just uniformly divided into a number of particles with same mass mb of spacing δr0 as shown in Fig. 12.7. From the integral approximation in Eqs. (12.64) and (12.65b), it follows that the particle approximation of the mass density N N X   X mb 0 ρ 5 W δx 5 mb Wab : ab δx0 b51 b51

Hence, the particle approximation of the function FI(x), (I 5 1, 2, 3, . . ., 9) at each particle is expressed as δr

δr′ a–1

b a

a+1

Smoothing domain of point a FIGURE 12.6 Discretization of 1D domain for an integral approximation. 1D, One-dimensional.

δr′ a–1

a

Smoothing domain of point a FIGURE 12.7 Discretization in particle approximation.

a+1

(12.65c)

508

Fluid
Solid Interaction Dynamics



N N N X  X mb X Wab F ðxa Þ 5 F ðxb ÞWab δx0 5 F ðxb ÞWab 5 F ð x b Þ PN : ρb b51 b51 b51 b51 Wab

(12.66)

Here subindex I is neglected, and we have used mb being a constant for the uniform division of the domain. The main difference between the integral approximation in Eq. (12.64) and the particle one in Eq. (12.66) is the introduction of material property in Eq. (12.66). Error parameter for accuracy estimation Smoothing length and particle numbers are the two key factors affecting the accuracy of a kernel approximation. Hence, the accuracy of each kernel function is tested with different smoothing lengths and different particle numbers. The error parameter for accuracy is defined by a mean square root (MSR): vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X   ε5t ε2 ; εb 5 Fðxb Þ 2 Fðxb Þ ; (12.67) N b51 b where ε is the MSR error, εb is the error between the analytical results and SPH results at the bth particle, and N is the number of particles used for error analysis. In order to further analyze the relationship between the error and the number of neighboring particles N, as well as the value of smoothing length h, we respectively assume that εN 5 AN α ;

log εN 5 log A 1 α log N;

(12.68)

ln εh 5 ln B 1 β ln h;

(12.69)

and εh 5 Bhβ ;

where α and β are variables to be determined, and A and B are constant coefficients. Based on these assumptions, when logarithmic scales are used, the relationship between ε and N will appear as a straight line, as shown in Eq. (1.68), and so is the case for ε and h, as shown in Eq. (12.69).

12.2.3.4 Numerical test results Integration approximation Fig. 12.8 shows that for the five functions in Eq. (12.63), the results obtained from the integral approximation in Eq. (12.64) agree well with the analytical data except for function F 5 tan(x) due to its singularity. In this case, the problem domain should be divided into two sections to avoid incorrect computation through the singular point. The results of the integral approximations with the fixed number of neighboring particles N 5 100 for the five functions in Eq. (12.63), influenced by the smoothing length h, are shown in Fig. 12.9, where the logarithmic scales based on Eq. (12.69) are adopted. From Fig. 12.9, it can be observed that in most cases decreasing the smoothing length h reduces the error. However, the function F 5 tan(x) is a special case, for which the error is high and fluctuates. It is difficult to improve the results by decreasing the smoothing length. The reason for this is that this function is singular in the domain. Overall, the super-Gaussian kernel function provides smaller errors, and the results are not sensitive to the change of smoothing length. Using the logarithmic scales, Fig. 12.10 shows the results of the errors for five functions affected by the neighboring particles when the smoothing length is fixed at h 5 0.1. From this figure, it is noted that, generally, the five neighboring particles can produce an approximation with an error of less than 5% for all kernel functions except for function F 5 tan(x). This means that the ratio of about k 5 2 between the spacing and the smoothing domain radius (kh) can provide accurate enough results for engineering applications. Increasing the number of neighboring particles can slightly improve the accuracy of the approximation for most kernels except cubic spline. It seems that the quartic kernel provides the smallest error in most cases and that cubic spline kernel is not sensitive to the decreasing smoothing length. Particle approximations The results of particle approximations for the five functions in Eq. (12.63) are shown in Fig. 12.11 and are compared with the corresponding theoretical results, which also agree well except for function F 5 tan(x) due to its singularity

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

F=x integral approxiamtion vs analytical

F=x3 integral approxiamtion vs analytical

2.5

10

Analytical results

1.5 1 0.5 0

6 4 2

0

0.5

1 1.5 Problem domian legth

0

2

0.5

1 1.5 Problem domian legth

2

1

Analytical results Integral approximation N=100,h=0.1

1.2

0

F=sin(x) integral approxiamtion vs analytical

F=exp(–x) integral approxiamtion vs analytical

1.4

0.8

1

Function vavlue

Function vavlue

Integral approximation N=100,h=0.1

8 Function vavlue

Function vavlue

Analytical results

Integral approximation N=100,h=0.1

2

509

0.8 0.6 0.4

0.6

0.4 0.2 Analytical results

0.2 0

0 0

0.5

1

1.5

2

Integral approximation N=100,h=0.1

0

0.5

Problem domian legth

1

1.5

2

Problem domian legth F=tan(x) integral approxiamtion vs analytical

60 Analytical results

Function vavlue

40

Integral approximation N=100,h=0.1

20 0 –20 –40 –60

0

0.5

1

1.5

2

Problem domian legth FIGURE 12.8 Results of integral approximations by new quartic for the five functions in Eq. (12.63), compared with the corresponding theoretical results (Sun, 2013).

510

Fluid
Solid Interaction Dynamics

0

Error vs smoothing length in integral approximation F=x

In(average error)

–5 –10 –15 –20 –25 –30 –35 –40 –4

–3

–3.5

–3 –2.5 ln(smoothing length)

–2

–1.5

Error vs smoothing length in integral approximation F=x3

In(average error)

–4 –5 –6 –7 –8 –9 –10 –4

–5

–3.5

–3 –2.5 ln(smoothing length)

–2

–1.5

Error vs smoothing length in integral approximation F=exp(-x)

Quadratic

In(average error)

–6 –7

Quartic

–8

New-quartic

–9

Johnson-quadratic

–10 –11

Super-Gaussian

–12

Cubic-spline

–13 –4

–4

–3.5

–3 –2.5 ln(smoothing length)

–2

–1.5

Gaussian

Error vs smoothing length in integral approximation F=sin(x)

In(average error)

–5 –6 –7 –8 –9 –10 –11 –12 –4

7

–3.5

–3 –2.5 ln(smoothing length)

–2

–1.5

Error vs smoothing length in integral approximation F=tan(x)

In(average error)

6 5 4 3 2 1 0 –1 –4

–3.5

–3 –2.5 ln(smoothing length)

–2

–1.5

FIGURE 12.9 The errors ln εh in Eq. (12.69) of integral approximations of five functions affected by the smoothing length ln h, with N 5 100 (Sun, 2013).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

Error vs the number of neighbouring particles in integral approximation F=x

Error vs the number of neighbouring particles in integral approximation F=x3

0

0 –1 Log (average error)

Log (average error)

511

–10

–20

–30

–2 –3 –4 –5

–40

2

4

6

8

10

16

14

12

18

–6

20

2

4

Number of neighbouring particles

–4

–4 Log (average error)

Log (average error)

–2

–6 –8 –10 –12 6

8

10

12

14

16

10

12

14

16

18

20

Error vs the number of neighbouring particles in integral approximation F=sin(x)

–2

4

8

Number of neighbouring particles

Error vs the number of neighbouring particles in integral approximation F=exp(-x)

2

6

18

Number of neighbouring particles

20

–6 –8 –10 –12 2

4

6 8 10 12 14 16 Number of neighbouring particles

18

20

Error vs the number of neighbouring particles in integral approximation F=tan(x)

3.5

Log (average error)

3 2.5 2 1.5 1 0.5

2

4

6

8 10 12 14 16 Number of neighbouring particles

18

20

FIGURE 12.10 The errors log εN of integral approximations of five functions affected by the neighboring particle number N, with fixed h 5 0.1 and legend the same as shown in Fig. 12.9 (Sun, 2013).

point. The influence of smoothing length and the number of particles on errors in particle approximation are also studied, and the results are shown in Figs. 12.12 and 12.13, respectively. From Fig. 12.12, it is clear that a smaller smoothing length provides smaller error. In the case of linear function approximation, the error is very small, although the curves are not showing a clear trend, which implies that the particle approximation of the liner function is normally very accurate even with a large smoothing length. Therefore decreasing the smoothing length is not meaningful in this case. It seems that the super-Gaussian gives the best accuracy in this particle approximation and quartic is the second best. Comparing the integral and particle approximations, we have found that the error is reduced with smaller smoothing length. The results can also be improved by increasing the number of particles, but this is not as efficient as using smaller smoothing length. This indicates that a proper smoothing length is more important to obtain accurate approximations. For function F 5 tan(x), as the value of the function approaches infinity, its singular point within the problem

512

Fluid
Solid Interaction Dynamics

F=x particle approxiamtion vs analytical

2.5

F=x3 particle approxiamtion vs analytical

10 Analytical

Function value

Analytical

Particle approximation N=100,h=0.1

2 1.5

6

1

4

0.5

2

0

0

Particle approximation N=100,h=0.1

8

0.5 1 1.5 Problem domain length

0

2

0

F=exp(-x) particle approxiamtion vs analytical

F=sin(x) particle approxiamtion vs analytical

1.4

1 Analytical Particle approximation N=100,h=0.1

1.2

0.8

1

Function value

Function value

2

1.5 1 0.5 Problem domain length

0.8 0.6

0.6 0.4

0.4

0

Analytical Particle approximation N=100,h=0.1

0.2

0.2 0

0.5 1 1.5 Problem domain length

2

0

0

0.5 1 1.5 Problem domain length

2

F=tan(x) particle approxiamtion vs analytical 100

Function value

50 0 –50 Analytical

–100

–150

Particle approximation N=100,h=0.1

0

0.5 1 1.5 Problem domain length

2

FIGURE 12.11 Results of particle approximations by new quartic for five functions in Eq. (12.63), compared with the corresponding theoretical results (Sun, 2013).

domain, it is difficult to obtain an accurate approximation. In this case, the problem domain can be divided into two sections, and error can be assessed in each section. In both the integral and particle approximations, all the kernels produce similar results; the choice of kernel is not very important in most cases. As shown in Figs. 12.9 and 12.12, the derivatives of error with respect to smoothing length are approximately constant for the functions x3 ; e2x , and sin x in both integral and particle approximations, which implies that the second derivative of error with respect to smoothing length approximates zero and that SPH approximation has a second order

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

Error vs smoothing length in particle approximation F=x

Error vs smoothing length in particle approximation F=x3

–35

–2 –3 In(average error)

–35.5 In(average error)

513

–36 –36.5

–4 –5 –6 –7

–37 –8 –37.5 –3.4 –3.2

–3

–2.8 –2.6 –2.4 –2.2 –2 ln(smoothing length)

–9 –3.4 –3.2

–1.8 –1.6 –1.4

–6

–5

–7

–6

–9

–1.8 –1.6 –1.4

–7 –8 –9 –10

–11 –12 –3.4 –3.2

–2

Error vs smoothing length in particle approximation F=sin(x) –4

In(average error)

In(average error)

Error vs smoothing length in particle approximation F=exp(-x)

–10

–2.8 –2.6 –2.4 –2.2

ln(smoothing lenghth)

–5

–8

–3

–3

–2.8 –2.6 –2.4 –2.2

–2

ln(smoothing length)

–1.8 –1.6 –1.4

–11 –3.4 –3.2

–3

–2.8 –2.6 –2.4 –2.2 –2 ln(smoothing length)

–1.8 –1.6 –1.4

Error vs smoothing length in particle approximation F=tan(x) 2.8

In(average error)

2.7 2.6 2.5 2.4 2.3 2.2 2.1 2 –3.4 –3.2 –3 –2.8 –2.6 –2.4 –2.2 –2 –1.8 –1.6 –1.4 ln(smoothing length)

FIGURE 12.12 The errors ln εh in Eq. (12.69) of particle approximations of five functions affected by the smoothing length ln h, with N 5 100 and legend the same as shown in Fig. 12.9 (Sun, 2013).

of accuracy with smoothing length. From Figs. 12.10 and 12.13, we have found that the influence of the total number of particles on the errors of the particle approximations is less sensitive than the influence of the number of neighboring particles on the error of integral approximations. This indicates that increasing the neighboring particles is more efficient than increasing the total number of particles, and therefore the number of neighboring particles is an important factor affecting the accuracy. To summarize, integral approximation does not differ from particle approximation in most cases, which means that the particle approximation is consistent with integral approximation. Different kernel functions give similar approximations; especially for new quartic, cubic spline, quartic spline, and quintic spline functions, they provide close results in most cases. It seems that the quartic kernel function shows the best performance generally. SPH approximation has a second order of accuracy with smoothing length. Decreasing the smoothing length or increasing the number of neighboring particles is normally useful in improving accuracy.

514

Fluid
Solid Interaction Dynamics

Error vs the number of particles in particle approximation F=x3 –4

–35

–4.5 log(average error)

log(average error)

Error vs the number of particles in particle approximation F=x –34.5

–35.5 –36 –36.5 –37 –37.5 30

–6 –6.5 –7 –7.5

40

50

60

70 80 90 Number of particles

100

110

120

–7

–6.5 log(average error)

–6.5

–6

–7.5 –8 –8.5 –9

40

50

60 70 80 90 Number of particles

100

110

120

–7 –7.5 –8 –8.5 –9

–9.5 –10 30

30

Error vs the number of particles in particle approximation F=sinx

Error vs the number of particles in particle approximation F=exp(-x)

log(average error)

–5 –5.5

40

50

60 70 80 90 Numer of particles

100

110

–9.5

120

30

40

50

60 70 80 90 Number of particles

100

110

120

Error vs the number of particles in particle approximation F=tanx

log(average error)

5 4 3 2 1 30

40

50

60

70

80

90

100

110

120

Number of particles

FIGURE 12.13 The errors log εN of particle approximations of five functions affected by the particle number N, with fixed h 5 0.1 and legend the same as shown in Fig. 12.9 (Sun, 2013).

12.2.3.5 Choosing smoothing length and symmetrization of particle interaction As the numerical test results for nine kernel functions presented in Sections 12.2.3.3 and 12.2.3.4, the smoothing length h is very important in the SPH method, which directly influences the efficiency of the computation and the accuracy of the solution. Physically, in practical simulations, if the smoothing length is too small, there may not be enough particles in the support domain κh to exert forces on a given particle, which results in low accuracy. Also, if the smoothing length is too large, all the details of the particle or local properties may be smoothed out, and the accuracy suffers too. For a given accuracy, the number N of the particle inside the effective support usually is fixed. Practice shows that for 1D, 2D, and 3D problems, the number of neighboring particles and the particle itself should be about 5, 21, 57, respectively, if the particles are placed in a lattice with a smoothing length of 1.2 times the particle spacing and κ 5 2. However, in dynamic simulations, most particles are moving, and particle density varies with time, so that for a given

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

515

accuracy and computation efficiency, the smoothing length should vary with time as well. So, for problems where the fluid expands or contracts locally to maintain consistent accuracy throughout the space, a different smoothing length could be used for each particle. Therefore it is crucial to set a criterion to determine the smoothing length to be adjusted during computations in order to maintain numerical accuracy and efficiency. One way to obtain this criterion is to require that each particle has approximately the same amount of mass. For example, under the distribution of the Gaussian kernel function in Eq. (12.62d), from Eq. (12.65c), the mass density at point a is given by N N X   X ρa 5 mb Wab 5 mb b51

and the mass of particle a satisfies the equation

b51

1 ðπh2 Þn=2

e2ððxa 2xb Þ=hÞ ;

    ma 5 ρa ΔVa ~ ρa hna ;

2

(12.70a)

(12.70b)

so that, to keep particle mass ma as a constant for all particles in the supporting domain, the smoothing length ha should satisfy  1=n : (12.70c) ha ~ ρ a This can be implemented based on the following nondimensional equation by considering the averaged density to update the smoothing length, that is,  1=n h ρ 5 0 ; (12.71a) h0 ρ where the subindex 0 indicates the initial values of the smoothing length and the density, respectively, and n is the number of dimensions of the problem. Another proposed method to evolve the smoothing length is by solving the following deferential equation: dh h dρ 52 ; dt dρ dt

(12.71b)

which implies that the smoothing length should be small in the region of the density increasing with time. The updated smoothing length approach leads to its variation in time and space, which results in each particle having its own smoothing length. For two interaction particles I and J with the corresponding different smoothing lengths hI and hJ, the influencing domain of particle I may cover particle J but not necessarily vice versa. Therefore the interaction between these two particles is a violation of Newton’s third law. In order to overcome this problem, it is necessary to find a way to preserve the symmetry of particle interactions. The proposed approach is to choose a common smoothing length hIJ for the interaction of two particles by using the arithmetic or geometric mean: Arithmetic: Geometric:

hI 1 hJ ; 2 2hI hJ hIJ 5 ; hI 1 hJ hIJ 5

(12.72a)

and the choice may be hIJ 5 maxðhI ; hJ Þ or hIJ 5 minðhI ; hJ Þ: The value of the smoothing function can be obtained using the chosen smoothing length WIJ 5 WðRIJ ; hIJ Þ;

(12.72b)

where RIJ is the distance between these two particles. Another approach is to use the average smoothing function value W IJ 5 without using a symmetric smoothing length.

½WðhI Þ 1 WðhJ Þ ; 2

(12.72c)

516

Fluid
Solid Interaction Dynamics

12.2.3.6 Smoothed particle hydrodynamics approximation rules To transform the PDEs of the conservation laws in continuum mechanics into SPH formulations, the following approximation rules are used: G

The average of the product of two functions f1 and f2 can be approximated by the product of the individual ensembles, that is,       (12.73) f1 Uf2 5 f1 U f2 :

G

The gradient of a scalar field function Φ can be approximated by the gradient of its kernel representation, that is, hrΦi 5 rhΦi:

(12.74a)

The approximation in Eqs. (12.74a)
(12.74d) may be an exact relation in an infinite domain, which, by using the Green theorem, is proved as follows: ð ð 0 0 0 hrΦi 5 Φ;i ðx ÞWðx 2 x ; hÞdx 5 f½Φðx0 ÞWðx2x0 ; hÞ;i 2 Φðx0 ÞW;i ðx 2 x0 ; hÞgdx0 ðΩ ðΩ ð (12.74b) @W @W 0 0 0 5 Φðx ÞWðx 2 x ; hÞν i dx 2 Φðx0 Þ 0 ðx 2 x0 ; hÞdx0 5 Φðx0 Þ ðx 2 x0 ; hÞdx0 5 rhΦi: @xi @xi S Ω Ω Here, we have considered the first integration on the surface as vanishing. Eq. (12.74a) provides a very useful relationship, hrΦi 5 rhΦi 2 ΦhrW i;

(12.74c)

since the kernel function W is localized, and hrW i 5 0. Actually, based on Eq. (12.74a) and the utility Eq. (12.7), we have ð hrW i 5 rhW i 5 r Wð0ÞWðx 2 x0 ; hÞdx0 5 rWð0Þ 5 0: (12.74d) Ω

12.2.4

Smoothed particle hydrodynamics formulation for Navier
Stokes equations

12.2.4.1 Particle approximation of density The density approximation is very important in the SPH method since it basically determines the particle distribution and the smoothing length evolution. The following two approaches can be used to obtain the particle approximation of density. Summation density This approach directly applies the SPH approximation to the density itself. For a given particle I, the density is obtained by Eq. (12.49), that is, N X ρI 5 mJ WIJ ; (12.75) J51

where N is the number of particles in the support domain of particle I, mJ is the mass of particle J, and WIJ 5 WðxI 2 xJ ; hÞ 5 W ðjxI 2 xJ j; hÞ 5 WðRIJ ; hÞ; jxI 2 xJ j ; RIJ 5 h

(12.76)

representing the smoothing function of particle I evaluated at particle J. Here, WIJ has a unit of the inverse of volume. Eq. (12.75) states that the density of a particle can be approximated by the weighted average of the densities of the particles in the support domain of that particle. The summation density approach conserves the total mass M in the entire unchanged problem domain Ω, which may be proven as follows. Based on the unity condition of the smoothing function, we have ð ð ð 0 0 0 M 5 ρðxÞdVðxÞ 5 ρðx ÞWðx 2 x ÞdVðx Þ dVðxÞ Ω Ω Ω ð (12.77) ð ð 0 0 0 0 0 Wðx 2 x ÞdVðxÞ dVðx Þ 5 ρðx ÞdVðx Þ 5 M; 5 ρðx Þ Ω

Ω

Ω

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

517

implying that the total mass is unchanged. However, when Eq. (12.75) is applied to the particles on the boundary of the problem, given the boundary particle deficiency, the result will be spurious. This edge effect could be remedied by using boundary virtual particles discussed later in the chapter. Another disadvantage for the density summation approach is that it entails more computational cost to evaluate the density before other parameters are evaluated. This method well represents the essence of the SPH approximation, so that it seems more popular in practical applications. For simulating general fluid phenomena, this approach yields better results. To improve the accuracy near both the free boundaries and the material interfaces with a density discontinuity, Eq. (12.75) is normalized with the SPH summation of the smoothing function itself over the neighboring particles as follows (Randles and Libersky, 1996): PN mJ WIJ ρI 5 PN J51 : (12.78) ðm J =ρJ ÞWIJ J51 Continuity density This approach approximates the density based on the continuity equation Dρ 5 2 ρvi;i ; Dt

(12.79)

in association with some transformation. Different transformations lead to different forms of density approximation equations. Form I: Using Eqs. (12.47) and (12.48), we approximate the velocity divergence in Eq. (12.79) so that N N X X DρI mJ J @WIJ mJ J @WIJ 5 ρI vi 5 2 ρ v : I J Dt ρ ρ i @xIi @x i J51 J J51 J

(12.80)

Form II: Consider the following expression of the derivative of the unit 0 5 1;i 5

N X mJ @WIJ J51

ρJ @xIi

;

(12.81)

from which it follows 0 5 ρI vIi

N X mJ @WIJ J51

ρJ @xIi

5 ρI

N X mJ J51

ρJ

vIi

@WIJ : @xIi

(12.82)

The summation of Eqs. (12.80) and (12.82) leads to another form of density approximation N X DρI mJ IJ @WIJ 5 ρI v ; Dt ρ i @xIi J51 J

 I  J vIJ i 5 vi 2 vi :

(12.83)

This particle approximation introduces velocity difference in the SPH formulations, which is preferred due to accounting for the relative velocities of particle pairs in the support domain. Form III: Using the continuity equation in the form Dρ 5 2 ρvi;i 5 2 ½ðρvi Þ;i 2 vi ρ;i ; Dt and particle approximation formulation in Eq. (12.47), we obtain 2 3 N DρI X m @W m @W J IJ J IJ 4 ρ J vI 5 5 2 ρ vJ i Dt ρJ ρJ J i @xIi @xIi J51 N X @WIJ 5 mJ vIJ : i @xIi J51

(12.84)

(12.85)

This equation shows clearly that the time change rate of density of a particle is closely related to the relative velocities, weighted by the gradient of the smoothing function, between this particle and all the other particles in the support domain.

518

Fluid
Solid Interaction Dynamics

The continuity density approach gives an approximation of the material derivative of the density, from which the total mass over the problem domain might be different at each time step. For simulating events with strong discontinuity, such as an explosion or a high-velocity impact, this approach is preferred.

12.2.4.2 Particle approximation of momentum The conservation law of momentum in continuum mechanics is given by Dvi 1 5 σij;j ; ρ Dt

(12.86)

in which we assume that the external body force is neglected, and σij is the stress tensor in fluids. Similar to the continuity density approach, we can obtain different forms of momentum approximation equations. Form I: Using Eqs. (12.47) and (12.48), we approximate the stress divergence in Eq. (12.86) so that N DvIi 1X mJ J @WIJ 5 σ : ρI J51 ρJ ij @xIj Dt

(12.87)

Form II: Adding the following identity derived from Eq. (12.83), 05

N N m σI X σIij X mJ @WIJ J ij @WIJ 5 ; ρI J51 ρJ @xIj ρ ρ @xIj J51 I J

(12.88)

N X σIij 1 σJij @WIJ DvIi 5 mJ : Dt ρI ρJ @xIj J51

(12.89)

we obtain

This formulation is symmetric to particle I and J, which reduces the error arising from the particle inconsistency, and it is therefore frequently used to evolve momentum. Form III: Using the identity

  Dvi 1 σij σij 5 σij;j 5 1 ρ ; ρ Dt ρ ;j ρ2 ;j

(12.90)

and particle approximation formulation in Eq. (12.47), we obtain N N σJij @WIJ σIij X DvIi X mJ @WIJ 5 mJ 1 ρ I Dt ρJ ρJ @xj ρI ρI J51 ρJ J @xIj J51 ! N X σIij σJij @WIJ 5 mJ 1 2 : @xIj ρ2I ρJ J51

(12.91)

This is another symmetric formulation evolving momentum. Introducing the constitutive equation for Stokes fluid σij 5 2 pδij 1 τ ij ;

τ ij 5 μγ ij ;

2 γ ij 5 vi;j 1 vj;i 2 vk;k δij ; 3

(12.92)

where p is the isotropic pressure, τ ij denotes viscous stress tensor, and γ ij is the tensor of shear strain rate. For linear viscous fluids, the viscous stress is proportional to the γ ij through the coefficient of dynamic viscosity μ. Eqs. (12.89) and (12.91) can be rewritten, respectively, as follows: N N X X μI γ Iij 1 μJ γ Jij @WIJ DvIi pI 1 pJ @WIJ 52 mJ 1 mJ ; I Dt ρI ρJ @xi ρI ρJ @xIj J51 J51

and  I  N N X X μI γ Iij μJ γ Jij DvIi p pJ @WIJ 52 mJ 2 1 2 1 m 1 J Dt ρI ρJ @xIi ρ2I ρ2J J51 J51

!

@WIJ : @xIj

(12.93)

(12.94)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

519

In these two formulations, the SPH approximations for the pressure and the viscous stress are separately expressed. From Eq. (12.92), the SPH approximation of shear strain rate tensor is derived as ! N N X mJ J @WIJ 2 X mJ J @WIJ I J @WIJ δ vi 1 v v : (12.95) γ ij 5 2 ij j I I 3 J51 ρJ k @xIk ρ @xj @xi J51 J Again, from Eq. (12.81), it follows 0 5 vIj

N X mJ @WIJ J51

0 5 vIi

ρJ @xIi

N X mJ

5

J51

5

J51

ρJ

@WIJ ; @xIi

N mJ @WIJ X mJ I @WIJ 5 v ; I ρ ρ i @xIj @x j J51 J J51 J

N X mJ @WIJ

ρ @xIk J51 J

5

N X mJ J51

which, when substituted into Eq. (12.93), gives N X mJ

vIj

N X

0 5 vIk

γ Iij

ρJ

vIJ i

@WIJ @WIJ 1 vIJ j @xIj @xIi

!

ρJ

vIk

(12.96)

@WIJ ; @xIk

N 2 X mJ IJ @WIJ v : 2 δij 3 J51 ρJ k @xIk

(12.97)

This equation includes the relative velocity between two particles I and J. After the shear strain rate for two particles are calculated, the acceleration can be calculated by Eq. (12.93) or Eq. (12.94).

12.2.4.3 Particle approximation of energy The conservation law of energy states that the time change rate of energy inside an infinitesimal fluid cell equals the summation of the net heat flux into that cell, the time rate of work done by the body, and the surface forces acting on that fluid cell. If the heat flux and the body forces are neglected, the time change rate of internal energy e consists of the work done by the isotropic pressure in the volumetric strain and the energy dissipation due to the viscous shear forces and is governed by the following equation:  De p μ p μ 4 5 2 vj;j 1 γ ij vi;j 5 2 vj;j 1 γ ij ðvi;j 1 vj;i Þ 1 ðvi;j 2 vj;i Þ 2 vk;k δij Dt ρ ρ ρ 2ρ 3 (12.98) p μ γ ij γ ij ; 5 2 vj;j 1 ρ 2ρ where in the derivation process, the following identities have been used: γ ij ðvi;j 2 vj;i Þ 5 0;

γ ij δij 5 γ jj 5 0:

(12.99)

For the particle approximation of Eq. (12.98), Eq. (12.97) can be used for the second part, and the pressure work can be modeled in different ways as follows. Form I: Considering the continuity Eq. (12.84) and its particle expression in Eq. (12.85), we have p p p Dρ ; 2 vj;j 5 2 ð2 ρvj;j Þ 5 2 ρ ρ ρ Dt 2 Form II: The identity

N p @vIi pX @WIJ 5 2 mJ vIJ : i I ρ @xi ρ J51 @xIi

    p @ p @ p vj 1 vj 2 vj;j 5 2 ; ρ @xj ρ @xj ρ

(12.100)

(12.101)

520

Fluid
Solid Interaction Dynamics

leads to another approximation, 2

N X p @vIi mJ pJ IJ @WIJ 5 v : I ρ @xi ρ2J i @xIi J51

(12.102)

Form III: Combining Eqs. (12.100) and (12.102), we obtain the most popular form of SPH approximation for the pressure work   N p @vIi 1X pI pJ IJ @WIJ 2 5 m 1 : (12.103) v J ρ @xIi 2 J51 ρ2I ρ2J i @xIi In this form, the variables appear symmetrically. Form IV: Using Eq. (12.83), we have 2 Form V: Based on the identity,

N p @vIi pI X mJ IJ @WIJ 5 v : I ρ @xi ρI J51 ρJ i @xIi

 p 1 @ @p 2 vj;j 5 2 ðpvj Þ 2 vj ; ρ ρ @xj @xj

(12.104)

(12.105)

we derive 2

N p @vIi 1X pJ @WIJ 5 mJ vIJ : i I ρ @xi ρI J51 ρJ @xIi

(12.106)

Form VI: Combining Eqs. (12.104) and (12.106), we obtain another symmetrical form of SPH approximation for the pressure work N p @vIi 1X pI 1 pJ IJ @WIJ 5 mJ v : (12.107) 2 I ρ @xi 2 J51 ρI ρJ i @xIi The two most frequently used SPH approximation forms for energy equations are:   N DeI 1X pI pJ IJ @WIJ μ 5 mJ 1 1 I γ Iij γ Iij ; vi I 2 2 2 J51 Dt @xi 2ρI ρI ρJ N DeI 1X pI 1 pJ IJ @WIJ μ 5 mJ vi 1 I γ Iij γ Iij : I 2 J51 Dt ρI ρJ @xi 2ρI

(12.108)

According to reported references, there is no obvious difference between the results obtained from these two approximation forms (Liu and Liu, 2003).

12.2.5

Numerical techniques for fluid flows

12.2.5.1 Artificial viscosity To remove numeric oscillations in fluid simulations involving shock waves, artificial viscosity was introduced. One such simulation is the von Neumann
Richtmyer artificial viscosity (Hirsch, 1988):  aN Δx2 ρðvj;j Þ2 ; vj;j , 0 DN 5 (12.109) 0; vj;j $ 0; where aN is an adjustable nondimensional constant, and another is a linear artificial viscosity  aL Δxcρvj;j ; vj;j , 0 DL 5 0; vj;j $ 0;

(12.110)

where aL is an adjustable nondimensional constant, and c denotes the speed of sound. These two artificial viscosities are widely used in fluid simulations by many numerical methods, such as FDM, FVM, FEM (Hirsch, 1988). For SPH simulations (Liu and Liu, 2003), the introduced artificial viscosity is given by

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

8 2 > < 2 αcIJ ϕIJ 1 βϕIJ ; ρIJ DIJ 5 > : 0;

IJ vIJ j xj , 0

521

(12.111)

IJ vIJ j xj $ 0;

where φIJ 5

IJ hIJ vIJ j xj ; IJ 2 xIJ j xj 1 ϕ

ϕ 5 0:1hIJ ;

cIJ 5

cI 1 cJ ; 2

I J vIJ j 5 vj 2 vj

ρIJ 5

ρI 1 ρJ ; 2

hIJ 5

hI 1 hJ ; 2

(12.112)

I J xIJ j 5 xj 2 xj :

In these equations, α and β are constants set around 1.0. The term with α produces a bulk viscosity, while the term with β suppresses particle interpenetration at a high Mach number. The factor φ prevents numerical divergences when two particles are approaching each other. The artificial viscosity given in Eq. (12.111) has the dimension of the pressure, which will be added in the pressure terms of the SPH equations. Another artificial viscosity was introduced into SPH simulations in the form qI qJ DIJ 5 2 1 2 ; ρI ρJ 8     (12.113) < 2αh ρ c vI  1 βh2 ρ vI 2 ; vI , 0 I I I j;j I I j;j j;j qI 5 : 0; vIj;j $ 0: Here, the artificial viscosity depends on the divergence of the velocity field.

12.2.5.2 Artificial heat Artificial heat is introduced into an energy equation to deal with excessive heating under severe circumstances such as the wall heating in the classic example of a stream of gas being brought to rest against a rigid wall. This artificial heat is added in the form (Liu and Liu, 2003) N X q

@WIJ IJ eI 2 eJ xIJ ; 2 j @xI IJ IJ ρ x x 1 φ j J51 IJ j j    2 qI 1 qJ     ; qI 5 2 αhI ρI cI vIj;j  1 βh2I ρI vIj;j  ; qIJ 5 2

HI 5 2

(12.114)

which will be added in the energy equation. Taking into account the artificial viscosity and heat, a very popular set of SPH formulations for the NS equations is as follows: N DρI X @WIJ 5 mJ vIJ ; i Dt @xIi J51 0 1 N σIij σJij DvIi X @WIJ 5 mJ @ 2 1 2 1 DIJ A I ; Dt @xj ρ ρ I J J51 0 1 N DeI 1X pI pJ @WIJ μ 5 mJ @ 2 1 2 1 DIJ AvIJ 1 I γ Iij γ Iij 1 HI ; i 2 J51 Dt @xIi 2ρI ρI ρJ

(12.115)

DxIj 5 vIj : Dt

12.2.5.3 Spurious zero-energy mode For some types of elements in FEM, due to an insufficient accuracy of numerical integrations, a spurious zero-energy mode with no strain or volume change in the element results in zero stress, leading to no resistance to the element

522

Fluid
Solid Interaction Dynamics

W′

x

FIGURE 12.14 A regular particle distribution for a 1D case, which leads to zero-energy mode with W0 (0) 5 0. 1D, One-dimensional.

deformation (Bathe, 1996). For fluid analysis, the FDMs, such as a second-order central difference scheme, also suffer from this type of zero-energy mode for which the derivative of the variable at a certain grid point is zero, although the variable itself oscillates greatly. The same spurious mode problem also occurs in the SPH method when evaluating the derivatives. As shown in Fig. 12.14 for a 1D problem, the regular particle distribution leads to zero-energy mode for the first derivative calculation. To avoid the spurious zero-energy mode, the proposed method is to use two separate types of particles: one for velocity evaluation and the other for stress evaluation. Also, if the particles are not uniformly distributed, the summed contributions from these particles generally will not lead to zero derivatives.

12.2.5.4 Pressure calculation for incompressible fluids For solving compressible flows, the particle motion is driven by the pressure gradient, while particle pressure is calculated by the local particle density and internal energy by means of the state equation. However, for incompressible flows, the actual equation of the fluid state leads to an extremely small time step. Artificial compressibility To remedy this, some artificial state equation for incompressible fluids is proposed as follows (Hirsch, 1988). One of these equations is  γ ρ p5B 21 ; (12.116) ρ0 where a constant γ 5 7 is used in most circumstances, ρ0 is the reference density, and B is a problem-dependent parameter used to set a limit for the maximum change of the density, such as the initial pressure. The subtraction of 1 can remove the boundary effect for free surface flows. Another artificial compressibility equation is simpler in the form of the multiplication of the square of the speed of sound with the density, that is, p 5 c2 ρ;

(12.117)

which was used to simulate incompressible flows with a low Reynolds number. In the artificial compressibility technique, the sound speed is a key factor that deserves careful consideration, for which many discussions and investigations are available; see, for example, Hirsch (1988). Projection method Another way of computing pressure is to treat the fluid as truly incompressible and to use a Poisson equation to obtain the pressure values. This method was developed by Chorin (1967) and is called the projection method, which has proven very effective in the simulations by MPMs (Cummins and Rudman, 1999; Ataie-Ashtiani et al., 2008; Brown et al., 2001; Sun, 2013, 2016; Sun et al., 2017). The numerical process is as follows. For incompressible fluids, the mass density is a constant, and the corresponding mass conservation law requires the divergence of velocity vanishing, that is, vi;i 5 0:

(12.118)

Dvi 2 p;i 1 5 1 τ ij;j 1 gδ3i ; ρ Dt ρ

(12.119)

The momentum equation governing fluid flows is

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

523

where g is the gravitation acceleration along the negative direction of x3. The proposed time integration process for Eq. (12.119) in a time step from tn to tn11 is split into two parts as follows: n11=2

vi

5 vni 1

1 n τ Δt 1 gδ3i Δt; ρ ij;j

n11=2

vn11 2 vi i

5

21 n11 p Δt; ρ ;i

(12.120)

and n11=2

vi;i

5

1 n11 p Δt; ρ ;ii

n11=2

pn11 ;ii 5

ρvi;i Δt

:

(12.121)

The first part considers the effect from the body force and the viscous stress but not the pressure gradient, so that an n11=2 unknown velocity vi is assumed, which will be determined by the second part considering the pressure gradient. Taking a divergence of the second equation in Eq. (12.120) and considering it satisfied for the time being, that is, vn11 i;i 5 0, we obtain a Poisson equation of the pressure in Eq. (12.121). From Eqs. (12.45) and (12.52), it follows that N

 1X n11=2 IJ @WIJ mJ vi ; ρ J51 @xIi N

J @W 1X IJ pn11 mJ pn11 : ;ii 5 ;i ρ J51 @xIi

(12.122)

@p xIJ @p xIJ pJ 2 pI PIJ xIJ 52 i 52 i 5 2i ; @xJi rIJ @rIJ rIJ rIJ rIJ PIJ 5 pI 2 pJ ;

(12.123)

n11=2

vi;i

52

Using Eq. (12.48), we obtain pJ;i 5

which, when substituted into Eq. (12.122), gives pn11 ;ii 5

N N

IJ @W IJ 1X mJ pn11 1 X IJ n11=2 IJ xi @WIJ 5 2 m v ; J i 2 1 ϕ2 @xI ρ J51 rIJ Δt J51 @xIi i

N X mJ pn11 xIJ @WIJ IJ

J51

i

2 1 ϕ2 @xI rIJ i

N

IJ @W ρ X IJ n11=2 52 mJ vi : Δt J51 @xIi

(12.124)

This equation can be solved efficiently. Finally, the velocity and position of each fluid particle can be renewed for the next time step using the following equation: n1ð1=2Þ

vn11 5 vi i

1 2 Pn11 Δt ρ ;i

n1ð1=2Þ 5 vi

2 Δt

N X J51

mJ



 Pn11 Pn11 @WIJ I J 1 2 1 DIJ : 2 ρ ρ @xIi

(12.125)

Here, we have used Eq. (12.52) to express the pressure gradient and introduced the artificial viscosity. Kinematic constraint method Ellero et al. (2007) presented a SPH model for incompressible fluids. As opposed to solving a pressure Poisson equation in order to get a divergence-free velocity field, the incompressibility was achieved by requiring, as a kinematic constraint, that the volume of the fluid particles was constant. They used Lagrangian multipliers to enforce this restriction. These Lagrange multipliers play the role of nonthermodynamic pressures whose actual values are fixed through kinematic restriction. The authors used the SHAKE methodology familiar in constrained molecular dynamics as an efficient method for finding the nonthermodynamic pressure satisfying the constraints. The model was tested for several flow configurations.

524

Fluid
Solid Interaction Dynamics

12.2.5.5 Boundary treatment Free surface On the free surface, the kinematic and dynamic conditions given in Eqs. (3.131a)
(3.131c) should be stratified. The kinematical condition in Eq. (3.131a) physically implies that the free surface is a material surface that can be automatically satisfied by the moving particle method since the MPM is a Lagrangian method in which the material point motion is traced. The dynamic conditions given in Eqs. (3.131b) and (3.131c), which respectively corresponds to the cases considering and not considering the surface tension. Normally, it can be assumed that the atmospheric pressure is p0 5 0 as the initial reference state of the system; also, the dynamic pressure p derived from Eq. (3.131a) or Eq. (3.131b) should be added to the partials on the free surface identified by the following approaches. 2D case: The method presented by Koh et al. (2012), Sun et al. (2014, 2015a), and Sun (2016) can be adopted. As shown in Fig. 12.15, each particle is taken as the center of a circle discretized by 360 points (dashed line). If each of these points on a circle is completely covered by its neighbors, the particle at the center of the circle is recognized as an inner fluid particle; otherwise it is a particle on the free surface. For example, in Fig. 12.15, particle A is identified as a particle on the free surface, since it is one of the dashed points on the circle not covered by its neighbors, while particle B is an inner fluid particle. 3D case: The efficient multilevel strategy for free surface identification is widely used (Marrone, 2011; Marrone et al., 2011; Tamai and Koshizuka, 2014) and is implemented in two steps. Step 1 aims preliminarily to find the possible particles on the free surface using the neighbor particle number Np based on the condition Np , βN0 ;

(12.126)

where N0 is the neighbor particle number of a typical inner fluid particle, such as 32 for a uniform particle distribution, and β is a tuning parameter of value about 0.96. The particles satisfying Eq. (12.126) may be the ones on the free surface. Step 2 refines the searching by checking the geometry property. As shown in Fig. 12.16, a sphere of radius r 5 1.05r0 (r0 is the initial particle distance) and centered at a particle identified by Step 1 is drawn, where nz denotes a unit vector pointing from the center of the sphere to the direction of the sparsest particles. Around this vector, a circular patch on the surface of the sphere, determined by angle θ 5 π=4, is discretized by evenly distributed points. If these points are all covered by the same sphere of its neighbor particles, this particle is recognized as an inner particle; otherwise it is regarded as a free surface particle. Fixed solid boundary In most of the cases, the solid boundary in the particle method will normally be handled by the ghost (virtual) particle method that was used by Libersky and Petschek (1991) and Monaghan (1994). Liu and Gu (2001a,b) suggested using two types of virtual particles: type I virtual particles are located right on the solid boundary. Type II particles fill the boundary region outside the type I particles, which are constructed in the following way. For a real particle I located within the distance of κhI from the boundary, a vertical particle is placed symmetrically on the outside of the boundary.

FIGURE 12.15 Identification of the particles on the free surface (Sun et al., 2017).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

525

FIGURE 12.16 Identification of free surface particles for 3D cases (Sun et al., 2017). 3D, Three-dimensional.

These virtual particles have the same density and pressure as the corresponding real particles but opposite velocity. Since type II particles are not enough to fully prevent the real particles penetrating the solid boundary, the application of type I particles to produce a sufficient repulsive boundary force fIJx when a particle approaches the boundary is needed. The repulsive force is calculated using a similar approach for the molecular force of Lennar-Jones form (Liu and Liu, 2003; Sun 2013): ( h α  β i   2 ; r0 =rIJ # 1 =r 2 r =r D r xIJ =rIJ 0 IJ 0 IJ x fIJ 5 (12.127) 0; r0 =rIJ . 1; where α and β are usually chosen as 12 and 4, respectively, and D is a parameter related to the problem, which should be in the same scale of the square of the largest velocity; rIJ denotes the distance between the real particle I and the virtual type I particle J, and r0 is selected approximately close to the initial particle spacing. Type I virtual particles take part in the kernel and the particle approximations for the real particles, but the position and physical variables do not evolve in the simulation process. Type II virtual particles are produced symmetrically according to the corresponding real particles, so that its parameters are determined by the related real particles in the simulation process. Enforced essential boundary conditions SPH is based on nonlocal interpolants, which makes it difficult to incorporate prescribed or developed data on an essential boundary into the discrete governing equations (Morris et al., 1997; Randles and Libersky, 1996; Takeda et al., 1994). Extra care has to be taken in order to enforce the essential boundary conditions. An approach to introduce the prescribed boundary condition into the particle equation was presented by Li and Liu (2004) as follows. Suppose particle I is a boundary particle of which the supporting domain consists of three subdomains: an interior fluid domain Ωf , a boundary domain Γ, and an outside solid domain Ωs . On the boundary domain Γ on which particle I is located, the value v^ of a variable v is prescribed. Based on the partition of unity at position I, from Eq. (12.43) with ΔVJ ρJ 5 mJ , we obtain X X X 15 ΔVJ WIJ 1 ΔVJ WIJ 1 ΔVJ WIJ ; (12.128a) JAΩf

JAΓ

JAΩs

^ it follows from which, when multiplied by the prescribed variable v, X X X ^ J WIJ 1 ^ J WIJ 1 ^ J WIJ : vΔV vΔV vΔV v^ 5 JAΩf

JAΓ

On the other hand, the kernel approximation for the value vI of the variable v gives X X X vJ ΔVJ WIJ 1 vJ ΔVJ WIJ 1 vI 5 v^J ΔVJ WIJ : JAΩf

JAΓ

(12.128b)

JAΩs

JAΩs

(12.128c)

526

Fluid
Solid Interaction Dynamics

Since the particles on Ωs have to satisfy the prescribed boundary condition and the particles on Γ should equal the calculated value vI by kernel approximation, we have ^ J 5 v^J ; vΔV vI 5 vJ ;

JAΩs ; JAΓ;

from which, when Eq. (12.128b) is subtracted from Eq. (12.128c), we obtain P JAΩf ðvJ 2 v^ÞΔVJ WIJ P ; vI 5 v^ 1 1 2 JAΓ ΔVJ WIJ

(12.128d)

(12.128e)

which gives the kernel approximation formulation of the boundary particles with the prescribed boundary condition satisfied. In this equation, the outside particles are not involved.

12.2.5.6 Time integration The discrete SPH equations can be integrated with standard numerical integration approaches. The explicit time integration schemes are subject to the stability condition, which requires that the computational domain of dependence in numerical simulation should include the physical domain of dependence. Therefore the maximum speed of numerical propagation must exceed the maximum speed of the physical propagation, so that the time step is taken as   hi Δt 5 min : (12.129) c Introducing viscous and external force effects into Eq. (12.129), the corresponding time step is given by 2 3 h I 5; Δtd 5 min4 cI 1 0:6ðαcI 1 β max ϕIJ Þ Δtf 5 min

h i1=2 hI f1

;

(12.130)

Δt 5 minðλ1 Δtd ; λ2 Δtf Þ; where f is the magnitude of the force per unit mass (acceleration), and λ1 and λ2 are suggested as 0.4 and 0.25, respectively. Another expression for time step is proposed as  2 h μ (12.131) Δt 5 0:125 ; ν5 ; ρ ν for considering viscous diffusion. Here the kinetic viscosity v is involved.

12.2.5.7 Particle interactions SPH method uses the smoothing function defined in a compact support domain with a finite number of particles to represent a function by particle approximations. These particles are generally referred to as nearest neighboring particles or the integration particles for the involved particle. The all-pair search approach calculates the distance rIJ from a particle I to each and every particle J in the problem domain to find the nearest particles by rIJ , κh. If the symmetric smoothing length is employed, particle I is also within the support domain of particle J. Therefore the computational time in this all-pair search approach is very long, especially for cases with a large number of particles. Monaghan and Gingold (1983) suggested using cells as a bookkeeping device to substantially save the neighboring particle search time, which is called a linked-list algorithm implementing in the following procedure (Monaghan, 1985; Dominguez et al., 2010). For a given problem, a temporary mesh with uniform spacing κh is overlain on the problem domain. For a given particle I, its nearest neighboring particles can be only in the same grid cell or in the immediately adjoining cells. Therefore the search is confined to only 3, 9, and 27 cells for 1D, 2D, and 3D space, respectively, and the search time is greatly reduced (Wro´blewski et al., 2007).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

527

The SPH simulations consume quite a bit of time to identify the interacting particles, after which the simulations can be carried out in the summation process over all the interacting particles. The two pairwise storage arrays can be used to store a particle ID number and its interaction particle ID number, and the two smoothing function arrays can be used to store the calculated weighted function values and its derivative values for each pair particles, based on which following summations can be effectively implemented.

12.2.5.8 Open source smoothed particle hydrodynamics solvers The SPHysics home page (http://www.wiki.machenster.ac.uk/sphsics/indes.php) provides open source SPH solvers (http://dx.doi.org/10.1016/j.cageo.2012.02.028; http://dx.doi.org/10.1016/j.cageo.2012.02.029) that can be used worldwide. SPHysics is a platform of SPH codes inspired by the formulation of Monaghan (1992) and developed jointly by researchers at the Johns Hopkins University (United States), the University of Vigo (Spain), the University of Manchester (United Kingdom), and the University of Rome La Sapienza (Italy). Developed over a number of years primarily to study free surface flow phenomena where Eulerian methods can be difficult to apply, such as waves or the impact of dam breaks on offshore structures, they have announced three available codes. Interested readers may access their webpage for more detailed information and the available codes.

12.2.6

Improved methods based on smoothed particle hydrodynamics

12.2.6.1 Corrective smoothed particle method The SPH method has many advantages in computations as just discussed, but it also has met with difficulties. One problem is the tensile instability (Swegle et al., 1995), in which the particle motion becomes unstable when particles are in a tensile stress state, which could result in particle clumping or a complete blowup in the computation. For a onedimensional von Neumann stability analysis, Swegle et al. (1995) identified a criterion for tensile instabilities in SPH. A sufficient condition for the unstable growth is σ11 W;11 . 0;

(12.132)

where σ11 5 σxx and W;11 5 d2 W=dx2 represent the stress and the second derivative of the smoothing function, respectively. As pointed out by Swegle et al. (1995), tensile instability does not happen for problems with equations of state where no tensile stress is generated. Another difficulty involves the boundaries of the problem. On the boundaries, the truncated integral in the kernel approximation and the insufficient particles in the particle approximations cause one of difficulties called the boundary deficiency problem in SPH (Monaghan, 1992; Randles and Libersky, 1996), which results in lower density near or on the boundary and finally yields spurious pressure gradients on the surface (Morris et al., 1997). To remedy these difficulties, different methods have been proposed (Liu and Liu, 2003); for example, Randles and Libersky (1996) developed a normalization formulation for the density approximation as shown in Eq. (12.78). In these modified approaches, the corrective SPM (CSPM) developed by Chen et al. (1999a,b,c, 2000) is discussed as follows. Using the Taylor series expansion at point I of function f(xi), we obtain    I I x 2 x 2 x x i j     i j f;ijI 1 ?; (12.133a) f ðxi Þ 5 f xIi 1 xi 2 xIi f;iI 1 2 where, as used in this book, the subindex i denotes the tensor index. Multiplying Eq. (12.133a) by the smoothing function Wðxi 2 xIi Þ and integrating over the support domain Ω, as well as neglecting the derivative term, we obtain Ð f ðxÞWðx 2 xI ÞdΩ I I ; (12.133b) f ðxi Þ 5 f ðx Þ 5 ΩÐ I Ω Wðx 2 x ÞdΩ which is the corrective kernel approximation for the function f(xi). The corresponding particle approximation can be directly obtained by considering dΩJ 5 mJ =ρJ and replacing the volume integration by its summation from preceding equation, that is, PN ðmJ =ρJ ÞfJ WIJ I fI 5 f ðx Þ 5 PJ51 : (12.133c) N J51 ðmJ =ρJ ÞWIJ

528

Fluid
Solid Interaction Dynamics

In the similar way to multiplying Eq. (12.133a) by the first derivative W;j ðxi 2 xIi Þ of the smoothing function, the three first derivatives f;iI at point I can be obtained by solving the following three coupling equations: ð ð

      I I I f;i xi 2 xi W;j x 2 x dΩ 5 f ðxÞ 2 f ðxI Þ W;j x 2 xI dΩ: (12.133d) Ω

Ω

This equation can be rewritten in the matrix form 2 a11 a12 4 a21 a22 a31 a32 where

32 3 2 3 f;1 a13 b1 a23 54 f;2 5 5 4 b2 5; f;3 a33 b3

ð aji 5 ðxi 2 xIi ÞW;j ðx 2 xI ÞdΩ; ðΩ bj 5 ½f ðxÞ 2 f ðxI ÞW;j ðx 2 xI ÞdΩ:

(12.133e)

(12.133f)

Ω

Also, we can obtain the corresponding particle approximation in the form 2 32 3 2 3 f;1 a~11 a~12 a~13 b~1 4 a~21 a~22 a~23 54 f;2 5 5 4 b~2 5; f;3 a~31 a~32 a~33 b~3

(12.133g)

where a~ji 5 b~j 5

N X mJ 

J51 N X

ρJ

 xJi 2 xIi W;jIJ ;

mJ J f ðx Þ 2 f ðxI Þ W;jIJ : ρ J51 J

(12.133h)

Equations in Eqs. (12.133a)
(12.133h) give the fundamental formulae in the CSPM used to generate the discrete equations for PDEs of the conservation laws in continuum mechanics. It has been shown that CSPM successfully improves the boundary deficiency and reduces the tensile instability in the SPH method (Chen et al., 1999b).

12.2.6.2 Moving particle semiimplicit method The moving particle semiimplicit (MPS) method was developed by Koshizuka and Oka (1996) and has shown improvement in stability in numerical simulations (Kondo and Koshizuka, 2011). The MPS method is similar to the SPH method (Gingold and Monaghan, 1977; Lucy, 1977) in that both methods provide approximations to the strong form of the PDEs on the basis of integral interpolants. However, the MPS method applies simplified differential operator models based solely on a local weighted averaging process without taking the gradient of a kernel function. In addition, the solution process of the MPS method differs from that of the original SPH method as the solutions to PDEs are obtained through a semiimplicit prediction
correction process rather than the fully explicit process in the original SPH method. The basic formulations of MPS method are summarized as follows. Function approximation The interpolation of a field function f(x) is similar to the normalization formulation in Eq. (12.133b) and is given by f ðxI Þ 5

N X f ðxJ ÞWJI J6¼I

nI

;

where WJI is the kernel function as used in SPH, but here it takes the form  r =r; ð0 # r , re Þ; WJI 5 e 0; ðre , rÞ;

(12.134a)

(12.134b)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

529

5.0 W (r)

0.0

r/re

1.0

FIGURE 12.17 Weight or kennel function W(r).

where re is the radius of the interaction area, r 5 jxJ 2 xI j is the distance between the particles I and J, and N is the particle number in the support domain. This function is called the weight function, which does not satisfy the utility property in the MPS method. The notation nI is defined by nI 5

N X

WJI

(12.134c)

J6¼I

and is called a particle number density in MPS. Actually, the meaning of this parameter denotes the total weighted particle number in the support domain. For point I, the different particle J has a different contribution to the variable value at point I, which, measured by the weight function W as shown in Fig. 12.17, shows that the particle J near particle I has the larger contribution to the variable value at particle I. Furthermore, considering Eq. (12.49) and assuming that each particle mass mJ 5 m is a constant, we have ρI 5 mnI , implying that nI is the equivalent number of particles with the same mass m in the supporting domain around point I. The particle number in a unit volume, or by normalization, can be obtained by the following equation: nI NI 5 Ð ; (12.134d) ΩI WðrÞdΩ where the denominator is the weighed volume of the supporting domain. Therefore Eq. (12.134d) implies that the total weighted particle number is divided by the weighted volume, which gives the particle number in the unit volume around particle I. Assuming that the particles have the same mass m, the mass density of the fluid can be derived as mnI ρI 5 mNI 5 Ð : (12.134e) ΩI WðrÞdΩ Compared with Eq. (12.49), this is the normalization formulation for mass density. Modeling of incompressibility The particle number density is the equivalent number of particles, of which each carries the same mass m, in the supporting domain of the same volume; therefore a constant n0 implies an incompressible fluid, and the continuity equation is satisfied. If the particle number density n is not n0, it is implicitly corrected to n0 by adding a correction value n^ satisfying the mass conservation equation n 2 n0 2 n~ 5 5 rU~v; Δtn0 Δtn0

n 1 n~ 5 n0 ;

(12.134f)

where v~ denotes the correction value of velocity, and its divergence is the volume deformation rate, of which the positive value gives the volume increment rate. If n~ 5 n0 2 n , 0, the particle number increases corresponding to a compression case with a negative volume deformation rate at state (*) to correct it; the divergence of the correction velocity v~ should contribute a positive divergence as shown in Eq. (12.134f). Using the relationship between the pressure gradient and the velocity given in Eq. (4.9a), the velocity correction value is derived from the implicit pressure gradient term as v~ 5 2

Δt n11 rp ; ρf

(12.134g)

530

Fluid
Solid Interaction Dynamics

from which, with Eq. (12.134f), a Poisson equation of pressure is obtained: 

ρ f nI 2 n0 : Δt2 n0

ðr2 pn11 ÞI 5 2

(12.134h)

The right side is represented by the derivation of the particle number density from its constant value. Derivative approximation The derivative of a function f(x) from points I to J is calculated by   df fJ 2 fI ðfJ 2 fI ÞðrJ 2 rI Þ 5 5 ; dr I rJ 2 rI jrJ 2rI j2

(12.135a)

from which the averaged gradient of function f(x) between particle I and its neighboring particles J is weighted as follows: ðrf ÞI 5

N dX ðfJ 2 fI ÞðrJ 2 rI Þ WJI ; n0 J6¼I jrJ 2rI j2

(12.135b)

where d is the number of space dimensions. In the numerical code by Koshizuka and Oka (1996), the final gradient formulation is   N fJ 2 fImin ðrJ 2 rI Þ dX WJI ; (12.135c) ðrf ÞI 5 0 n J6¼I jrJ 2rI j2 where fImin 5 min fJ ;

(12.135d)

for all neighboring points J in the supporting domain. Approximation of Laplacian A time-dependent diffusion equation with the diffusion coefficient v is represented in Laplacian terms as @f 5 νr2 f : @t

(12.136a)

The total variance of the distribution of function by the total particles n0 is proportional to 2dνΔt during time step Δt, where d is the number of space dimensions, of which the averaged contribution by each particle is proportional to 2dνΔt=n0 . The part of quantity fI of particle I is distributed to the neighboring particle J according to the kernel function. Therefore the quantity transferred from particle I to J is assumed to be ΔfIJ 5 fJ 2 fI 5

2dνΔt fI WJI ; n0 λ

(12.136b)

where λ is a coefficient with the dimension of length square due to the second derivatives with respect to the coordinate involved on the right-hand side of Eq. (12.136c), which is taken as the averaged second moment of supporting domain of the kernel function, that is, Ð WðrÞr 2 dΩ : (12.136c) λ 5 ΩÐ I ΩI WðrÞdΩ In this model, the right-hand side in Eq. (12.136b) represents the Laplacian contribution, multiplied by νΔt, of the quantity distributed (output from I) from I to J. Therefore the quantity at point I is reduced, so that we can write Eq. (12.136b) in the following form: ðr2 f ÞI 5 2

2d fI WJI ; n0 λ

(12.136d)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

531

where a negative notation is added on the right-hand side to denote that the quantity at point I is reduced. Similarly, we can represent the quantity distributed from particle J to I (input into I) in the forms ðr2 f ÞI 5

2d 2d fJ WIJ 5 0 fJ WJI ; 0 n λ n λ

(12.136e)

where the positive notation implies the quantity input into particle I. We have used the relation WIJ 5 WJI. Considering all neighboring particles, the approximation of Laplacian contribution by all particles J to I is represented by a summation of Eqs. (12.136d) and (12.136e) overall particles, that is ðr2 f ÞI 5

N X 2d ðf 2 fI ÞWJI : 0λ J n J6¼I

(12.136f)

The approximation formulations in Eqs. (12.134a)
(12.134h), (12.135a)
(12.135d), and (12.136a)
(12.136f) in MPS for any field functions and their derivatives with Laplacian operation provide the fundamental techniques to obtain the conservation laws in the MPS model, which is similar to the SPH process, neglected here. Interested readers may refer to the original paper by Koshizuka and Oka (1996), as well as further investigations and modifications, such as Javed et al. (2013a,b, 2014a,b); Javed (2013); Sun et al. (2015a,b); and Sun (2016).

12.3

Meshfree Galerkin methods

In this subsection, we briefly present the fundamental formulation and techniques of MGMs. As an example, we summarize only the essential concept of the moving least square (MLS) method, which seems to be the basis of many other various approaches. This method uses the MLS method to generate the kernel interpolant and then, based on the Galerkin weak or variational formulations of PDEs, to seek the approximation solutions of the problems. Detailed information on these formulations with mathematical proofs may be read in Liu and Liu (2003) and Li and Liu (2004).

12.3.1

Moving least square representing kernel interpolant

The MLS method has two main technical characteristics: one is the local standard least square procedure, and the other is the moving process from local region to the whole region for a global approximation. The moving process behavior was not clearly stated until it was later elaborated by Liu et al. (1997) and again by Li and Liu (1999a,b). The detailed formulations of the MLS methods can be read in the publications by Li and Liu (1996a, 1997, 2004), which is briefly summarized as follows.

12.3.1.1 Local standard least square approximation Let uðx; x~ Þ be a sufficiently smooth function defined on a simply connected open domain in the space, that is, ΩARn . Considering a fixed point x~ AΩ, we can define a local function  uðx; x~ Þ; xABlρ ð~xÞ l u ðx; x~ Þ 5 (12.137a) 0; otherwise; where Blρ ð~xÞ denotes an effective spherical ball with its dilation parameter ρ 5 ð1B3Þ 3 particle space that is similar to the smoothing length of SPH. According to the well-known Stone-Weierstrass theorem (Rudin, 1976), we are able to approximate u(x) by a local polynomial series:     l X ~ x 2 x~ l T x2x u ðx; x~ Þ 5 Pα Dα ð~x; ρÞ 5 P Dð~x; ρÞ; (12.137b) ρ ρ α51 where

T Dð~x; ρÞ 5 D1 D2 ? Dl ;

T PðxÞ 5 P1 ðxÞ P2 ðxÞ ? P3 ðxÞ :

(12.137c)

D is an unknown vector, and Pα ðxÞ is an independent polynomial basis with P1(x) 5 1. Since the polynomial series can be sufficiently smoothed by choosing a suitable term number l, the local approximation of the function in

532

Fluid
Solid Interaction Dynamics

Eq. (12.137b) is sufficiently smooth. For vectors: PðxÞ 5 1 PðxÞ 5 1 PðxÞ 5 1

examples, in 1D, 2D, and 3D cases, we have the following polynomial basis

T x x2 ? xl21 ;

T x1 x2 x21 x1 x2 x22 ; x1 x2 x3 x21 x22 x23 x1 x2

x2 x3

x3 x1

T

(12.137d) :

The error of the approximation in Eq. (12.137b) for the function u(x) at point I is given by   ~ T xI 2 x Dð~x; ρÞ 2 uðxI Þ; rðxI ; x~ Þ 5 P ρ

(12.138a)

of which the weighted square residual is defined as RðDÞ 5

N X

r 2 ðxI ; x~ ÞWρ ð~x 2 xI Þ;

(12.138b)

I51

where Wρ ð~x 2 xI Þ denotes the kernel function as a weight at point I. To minimum the residual in Eq. (12.138b), the standard least square method can be used to obtain the optimal coefficient vector D, satisfying     N X ~ 2 xI x~ 2 xI T x P (12.138c) uðxI Þ 2 P Dð~x; ρÞ Wρ ð~x 2 xI Þ 5 0; ρ ρ I51 from which we obtain

  N X x~ 2 xI P uðxI ÞWρ ð~x 2 xI Þ; ρ I51

(12.138d)

    N X x 2 xI T x 2 xI P P Wρ ðx 2 xI Þ; ρ ρ I51

(12.138e)

Dð~xÞ 5 M21 ð~xÞ and the moment matrix defined by MðxÞ 5 or in an integration form



ð MðxÞ 5

P Ω

   x2y T x2y P Wρ ðx 2 yÞdΩy : ρ ρ

(12.138f)

The moment matrix is always positive since Pα ðxÞ ðα 5 1; 2; . . .; lÞ are linearly independent, so that it is invertible. Substituting Eq. (12.138d) into Eq. (12.137b), we obtain the local standard weighted least square approximation of the function in the form     N X x~ 2 x x~ 2 xI ul ðx; x~ Þ 5 PT uðxI ÞWρ ð~x 2 xI Þ: M21 ð~xÞ P (12.138g) ρ ρ I51

12.3.1.2 Moving process for global approximation Eq. (12.138f) is optical in a local region Blρ ð~xÞ. To extend this approximation to the whole region, a moving process is ensured, which requires the fixed point x~ -x, that is,   N X uðxÞ 5 lim ul ðx; x~ Þ 5 PT ð0ÞM21 ðxÞ P x 2 xI uðxI ÞWρ ðx 2 xI Þ: (12.139a) x~ -x ρ I51 We define the correction function



 x 2 xI Cρ ðx; xI Þ 5 P ð0ÞM ðxÞP ρ     x 2 xI x 2 xI 5 BT ðx; ρÞ P 5 PT Bðx; ρÞ; ρ ρ T

21

Bðx; ρÞ 5 M21 ðxÞPð0Þ;

(12.139b)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

533

so that we can obtain the kernel function produced by the MLS process κρ ðx 2 xI ; xÞ 5 Cρ ðx; xI ÞWρ ðx 2 xI Þ;

(12.139c)

and the corresponding function approximation uρ ðxÞ 5

N X

uðxI Þκρ ðx 2 xI ; xÞ;

(12.139d)

uðyÞκρ ðx 2 y; xÞdΩ:

(12.139e)

I51

or its integration form

ð uρ ðxÞ 5

Ω

12.3.1.3 An example It would be helpful to examine Eqs. (12.139a)
(12.139e) through a pedagogic example. Now we consider the case of the linear polynomial basis

T

T Pð0Þ 5 1 0 0 0 ; (12.140a) PðxÞ 5 1 x1 x2 x3 ; of which the moment matrix can be derived using Eq. (12.138f), that is, 2 3 1 ð 6 7 x 2 y1 x1 2 y1 x3 2 y3 6 ðx1 2 y1 Þ=ρ1 7 1 1 MðxÞ 5 6 Wρ ðx 2 yÞdΩy ; 7 ρ1 ρ1 ρ3 Ω 4 ðx1 2 y1 Þ=ρ1 5

(12.140b)

ðx3 2 y3 Þ=ρ3 which may be rewritten in compact form:

where

ð  mαβ 5

Ω

m00 6m 6 10 MðxÞ 5 6 4 m20

m01 m11

m02 m12

m21

m22

3 m03 m13 7 7 7 5 MT ðxÞ; m23 5

m30

m31

m32

m33

xα 2 yα ρα

2



 xβ 2 yβ Wρ ðx 2 yÞdΩy ; ρβ

x0 2 y0 5 1; ρ0

(12.140c)

ð0 # α; β # 3Þ:

(12.140d)

The corresponding correction function in Eq. (12.139b) is expressed as 2 3 1 7

A 6 6 ðx1 2 y1 Þ=ρ1 7 Cρ ðx; xI Þ 5 1 0 0 0 6 7 jMj 4 ðx1 2 y1 Þ=ρ1 5 5

A11 

5 PT

A12 A13 jMj

ðx3 2 y3 Þ=ρ3 2 3 1

6 A14 6 ðx1 2 y1 Þ=ρ1 7 7 6 7 4 ðx1 2 y1 Þ=ρ1 5

 x 2 xI Bðx; ρÞ; ρ

ðx3 2 y3 Þ=ρ3

where the matrix A denotes the adjoint matrix of the moment matrix M(x) of which an element     ; α~ 6¼ α; β~ 6¼ β; Aαβ 5 ð21Þα1β mα~ β~  333

(12.140e)

(12.140f)

534

Fluid
Solid Interaction Dynamics

and

Bðx; ρÞ 5

12.3.2

A11

A12 A13 jMj

A14

:

(12.140g)

Shepard interpolant

Eq. (12.139d) is also considered as an interpolation of a function, which may be rewritten in the form uðxÞ 5

N X

uI NI ðx; ρÞ;

NI ðx; ρÞ 5 κρ ðx 2 xI ; xÞ;

(12.141a)

I51

which, when the polynomial series in Eq. (12.137d) is used, satisfies uðxI Þ 5 uI ;

(12.141b)

NI ðxJ Þ 5 δIJ ;

(12.141c)

but the interpolation shape function does not satisfy

which is exclusively reserved for the term interpolation. A special case is l 5 1 in which, when considering Eqs. (12.138a)
(12.138g) and (12.139a)
(12.139e), we have PðxÞ 5 1; N X Mð0Þ 5 Wðx 2 xI Þ;

1

M21 ð0Þ 5 PN

I51

I51

Wðx 2 xI Þ

; (12.141d)

Wðx 2 xI Þ κðx 2 xI ; xÞ 5 PN : I51 Wðx 2 xI Þ It is clear that N X

0 , κðx 2 xI ; xÞ , 1;

κðx 2 xI ; xÞ 5 1:

(12.141e)

I51

This particular interpolant was first studied by Shepard (1968), from which the interpolation in Eq. (12.141a) becomes PN uI Wðx 2 xI Þ uðxÞ 5 PI51 : (12.141f) N I51 Wðx 2 xI Þ As suggested by Shepard, a remedy for the interpolation property is to choose the singular weight function 1 ; jx2xI jα

Wðx 2 xI Þ 5

(12.142a)

where α is a positive even integer. Therefore the Shepard interpolant, based on Eq. (12.142a), is κS ðx 2 xI ; xÞ 5

1

P jx2xI jα NJ51

α

1=jx2xJ j

5

11

PN

J6¼I

1 jx 2 xI jα =jx 2 xJ jα

;

(12.142b)

which satisfies the following equations: κS ðxJ 2 xI ; xJ Þ 5 δIJ ; N X κS ðx 2 xI ; xÞ 5 1; I51

rκS ðx2xI Þx5J 5 0:

0 # κS ðx 2 xI ; xÞ # 1; lim κS ðx 2 xI ; xÞ 5

x-N

1 ; N

(12.142c)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

12.3.3

535

Orthogonal basis for local approximations

To reduce the time consumed in simulations, orthogonal basis would be a good choice. Based on the linear independent basis vector in Eq. (12.137c), the following orthogonal basis can be constructed: q1 ðx; x~ Þ 5 P1 ðx; x~ Þ 5 1;   p2 ; q1 x~  q1 ðx; x~ Þ; q2 ðx; x~ Þ 5 P2 ðx; x~ Þ 2  q1 ; q1 x~     p3 ; q1 x~ p3 ; q2 x~  q1 ðx; x~ Þ 2   q2 ðx; x~ Þ; q3 ðx; x~ Þ 5 P3 ðx; x~ Þ 2  q1 ; q1 x~ q2 ; q2 x~ ^ qα ðx; x~ Þ 5 Pα ðx; x~ Þ 2

  α21 X pα ; qβ x~   qβ ðx; x~ Þ; qβ ; qβ x~ β51

(12.143a)

α 5 1; 2; . . .; l:

Here, the following inner product definition is used ð   f ; g x~ 5 f ð~x 2 xÞgð~x 2 xÞWρ ð~x 2 xÞdΩ:

(12.143b)

Ω

Because of orthogonalization, the local moment matrix is diagonal, and its inversion is 1  : qα ; qα x~

(12.143c)



T Dð~xÞ 5 D1 D2 ? DN ; PN qα ðxI ; x~ ÞuI Wρ ð~x 2 xI Þ  : Dα ð~xÞ 5 I51  qα ; qα x~

(12.143d)

M21 ð~xÞ 5 diagðaα Þ;

aα 5 

Thus the vector in Eq. (12.138d) becomes

The global approximation correction function in Eq. (12.139b) is now given by Cρ ðx; xI Þ 5

l X qα ðx; xÞqα ðxI ; xÞ   ; qα ; qα x~ α51

(12.143e)

which can be used in Eqs. (12.139c) and (12.139d) to give the global approximation formulations based on the orthogonal basis in Eq. (12.143a). The Lagrangian polynomial (Waring, 1779) N x 2 xI Lðx; xI Þ 5 Π (12.144a) I50 xJ 2 xI I 6¼ J

is a typical example of an orthogonal polynomial that satisfies the following interpolation and orthogonal conditions: LðxJ ; xI Þ 5 δIJ ; N N X X ½Lðxα ; xJ ÞLðxα ; xI Þ 5 ½δαJ δαI  5 NδIJ : α51

12.3.4

(12.144b)

α51

Applications of the moving least square reproducing kernels

The kernel function given in Eq. (12.139c) can be used in numerical simulations in the following aspects: 1. By choosing different polynomial series and weight functions, this general equation produces the different kernel functions with sufficiently smooth characteristics to be used in MPMs, such as SPH and MPS discussed in Section 12.2. 2. This general function provides a shape function defined at point I, that is,

536

Fluid
Solid Interaction Dynamics

NI ðx; ρÞ 5 κρ ðx 2 xI ; xÞ;

(12.145a)

which can be used like the Rayleigh
Ritz functions in Section 1.4.3.3 or interpolation functions in FEA to obtain the approximation solutions of PDEs modeled by the variational principles, weighted residual method discussed in Chapter 1, Introduction. The references by Li and Liu (1996a, 1997, 2004) present various problems in continuum mechanics in the weak form of the functionals or Galerkin formulations and then solve them by the assumed solution uðxÞ 5

N X

NI ðx; ρÞQI ;

(12.145b)

I51

where QI denotes the node variables in analysis, which may be a time function for the dynamic problems. For examples, Belytschko et al. (1994a,b) adopted the MLS interpolant and the weak Galerkin formulation to simulate a crack growth problem in a linear elastic solid.

12.4

Mixed finite element—smoothed particle method for fluid
solid interaction problems

Now we discuss the proposed MFE
SPMs for nonlinear FSI problems. Compared with the FE
CFD method discussed in Chapter 11, Mixed finite element
computational fluid dynamics method for nonlinear fluid
solid interactions, the MFE
SPM method is more convenient in simulating problems with violent fluid motions, such as breaking waves and solid
fluid separations, since SPM models the fluids as many individual material particles that are easily traced even in violent flow cases. Here, SPM is used as a generalized term to include various MPMs, such as the SPH, MPS, and MLS methods.

12.4.1

Generalized solution procedure

As discussed in Section 11.3 for solving nonlinear FSI problems using the MFE
CFD approach, we also have two general options—the simultaneous integration and the partitioned iteration—to be chosen for the MFE
SPM dealing with nonlinear FSI problems. But the partitioned iterative procedure is suggested for combining the available powerful FEM codes for solid structures and well developed SPMs codes for fluids. Obviously, as reported in some references, it is possible to use SPM techniques to both fluids and solids modeling. However, the numerical practices show that for a solid continuum system defined in a 2D or 3D domain, it is not difficult to model it by the SPM, but for many complex structures often used in engineering, SPM modeling might not be effective. The combination of FEM for structures and SPM for fluids may be the most appropriate approach to solve nonlinear fluid
structure interaction problems involving violent fluid motions. In this combined method, the structures can be exactly modeled by the nonlinear FEM described in Section 11.1, and the fluids can be exactly simulated by the SPM discussed in Sections 12.2 and 12.3, as well as other by similar methods reported in the references. In this simulation, since both the fluid and the solid undergo their nonlinear mechanical process, the solutions for the fluid and the solid are also based on their iterative calculations. The flowchart of the numerical solution is given in Fig. 12.18, where the parameter α is a relaxing factor to modify the approximate intermediate solutions for numerical convergence purposes, which can be determined for different methods by experience or numerical tests in calculations. Based on theoretical analysis (Aitken, 1937; see also Gianola and Schaffer, 1987), this relaxation factor can be updated to α in the next iteration calculated by the increment, denoted by ΔðÞ, of the variables in the form α 5

2αðΔU ÞT ΔðU 2 UÞ   : ΔðU 2 UÞ

(12.146a)

The simulations consist of the following main steps: 1. Input the initial conditions: geometrical sizes, material constants, initial values of variables such as solid displacement and velocity, fluid pressure and velocity, boundary conditions, etc. 2. Do fluid SPM iterative calculations based on the prescribed initial solid information to obtain the fluid pressure P ðt 1 ΔtÞ, velocity, etc. 3. Do solid FE iterative calculations based on the fluid information obtained in step 2 to obtain the solid displacement U  ðt 1 ΔtÞ, velocity, and accelerations, etc.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

Start t=0

537

Yes Stop No Time end

Initial P(t),U(t)

Fluid SPM calculation P* (t + Δt)

U = U + α(U*–U)

t = t + Δt Yes

⏐⏐U*–U⏐⏐< ε

Solid FEM calculation U* (t + Δt)

No

FIGURE 12.18 Flowchart of the partitioned iteration of mixed FEM
SPM for nonlinear fluid
structure interaction dynamics. FEM, FE method. FE, finite element; SPM, smoothed particle method.

4. Make a convergence check to determine whether the error is less than the allowed error ε, that is, :U  2 U: , ε. Time goes to next time step; otherwise, modify the solid displacement to U 5 U 1 αðU  2 UÞ to repeat the calculations in steps 2 and 3 until convergence reached. 5. Check the time period of interest in the FSI problem. If it has not been reached, go to next time step to repeat steps 2
4 until it is reached. The detailed numerical implementation process depends on the models adopted for solids and fluids. For nonlinear fluid
structure interaction problems in which the structure undergoes motions with no large rigid motions but only large elastic deformations, the updated FE formulation could be effective in simulating its motions. If the structure undergoes large rigid motions with small, relatively elastic deformation, FEM may be not effective. For example, for a structure undergoing motions with no elastic deformation, there is no elastic deformation energy anywhere in the rigid body, so that there is no need to use elements describing the local deformation at every point. Therefore for nonlinear rigid body fluid interactions, it is more effective to use the methods in rigid body dynamics to model the motions in the body. In Section 11.3.1.2, we gave the mathematical formulations governing rigid body motions coupled with fluid flows. In maritime structures
water interactions, such as nonlinear ship
water interactions or missile into or out of the water problems, the structure normally undergoes very large rigid motions with relatively small elastic deformation. To model this type of nonlinear FSI problem, it is convenient to separate the rigid motion and small elastic motion to establish the governing equations describing the solid motion, discussed in the next subsection.

12.4.2 Modeling of fluid
solid interaction involving large rigid motions with small elastic deformation The generalized mathematical modeling of the MFE
SPM for violent water
structure interactions, with freak waves and with the structure undergoing a large rigid motion plus a small elastic deformation, is presented by Xing et al. (2016). In this model, the total displacement of the structure is expressed by a rigid motion of 6 degrees of freedom (DOF) plus a small displacement relative to the rigid motion. Since the elastic deformation is small, its deformation is governed by a linear elastic theory, and a mode summation approach is valid in reducing the final size in the numerical equations. The elastic mode shapes of elastic structure can be obtained by FEA or by the available theoretical mode functions for some structural members, such as rods, beams, shafts, plates, etc. The water is assumed to be inviscid and incompressible, and its motion is governed by a nonlinear NS equation. On the coupling interface where no FS

538

Fluid
Solid Interaction Dynamics

separation occurs, the equilibrium and consistency conditions are required. The governing equations of the nonlinear FSI systems are as follows.

12.4.2.1 Governing equations As shown in Fig. 12.19, a solid occupies a material domain Ωs of boundary S 5 ST , SU , Σ with its outside normal vector denoted by ηi . A traction force T^ i and a displacement U^ i are, respectively, given on part of the boundary ST and SU. The solid body is floating on the surface of the water, and no displacement boundary SU is given, of which the dry structure is considered a free
free body in the space. The solid body interacts on its wet interface Σ with the water in the domain Ωf of boundary Γ 5 Γf , Γb , Σ , Γv1 , Γv2 with its outside normal vector ν i . The dynamic pressure on the free surface Γf is assumed to be the atmosphere denoted by p 5 0, and the sea bottom Γb is considered a fixed boundary with the water velocity vi 5 0. The possible velocities v^i may be given on the boundaries Γv1 or Γv2 ; for example, an inflow velocity is given on Γv1 , while Γv2 is treated as an infinite boundary with no disturbance of v^i 5 0. Three coordinate systems are defined. A Cartesian system o 2 x1 x2 x3 with positive x3 in the vertical direction is fixed at a point o on the free surface in the static equilibrium state of the system. A moving system O 2 y1 y2 y3 , with its axis O 2 yi parallel to the axis o 2 xi , (i 5 1, 2, 3), respectively, is fixed at the mass center O of the solid, which moves with the mass center. Beside these two coordinate systems, a principal inertial material system O 2 X1 X2 X3 is defined, of which the origin and three axes are located at the same positions of the system o 2 x1 x2 x3 at time t 5 0. The position coordinates of a material point at time t 5 0 are denoted by Xj to identify this material particle, which moves to a new position xj at time t. The motion of this material point Xj is denoted by a summation of a translation Uio of the mass center, a rigid rotation θi about three axes, and small elastic displacement Uj of the structure. Using the notations of Cartesian tensors and summation convention, such as the Kronecker delta δij and permutation symbol eijk , and denoting the rigid rotation by a tensor Rij ðθk Þ with its time derivative and partial derivative with respect to θi given in Eqs. (10.156) and (10.157), that is, @Rij 5 eimk Rkj ; @θm

R_ij 5 eimk θ_ m Rkj ;

(12.146b)

from which, along with the preceding assumptions, we can derive the dynamic equations of the system as follows. Since the material coordinate system is fixed at the mass center of the body and it is chosen as the principal inertial system, its first-order inertial moment and the second-order cross-inertial products vanish, and the inertial matrix J is a diagonal matrix, that is,  ð ð 0; j 6¼ i ρs Xi dΩ 5 0; ρs Xi Xj dΩ 5 Jij 5 (12.146c) JII ; j 5 I 5 i: Ωs Ωs The total mass M of the body is calculated by the following integration over the body volume: ð ρs dΩ 5 M:

(12.146d)

Ωs

x3y3 X3

y3

θ3

X2

X3

P

p

x

x2y2 X2

O

FIGURE 12.19 Scheme of nonlinear FSI system involving large structural rigid motions with small elastic deformation (Sun, 2016). FSI, Fluid
solid interaction.

t=0

ρgg

x1y1 X1

Ωs

θ1 y1

Ωf Γv1

y2

O Σ

Γb

θ2

t=t1 X1

Γv2

ρfg Fluid field modeled by particle method

Γf

u

x

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

539

Solid domain Displacement and velocity fields The position and velocity of particle Xj at time t, respectively, are xi 5 Xi 1 Uio 1 Rij Xj 1 Rij Uj ; o x_i 5 U_ i 1 R_ij Xj 1 R_ij Uj 1 Rij U_ j o 5 U_ i

(12.147a)

1 eimk θ_ m Rkj ðXj 1 Uj Þ 1 Rij U_ j :

Here the elastic displacement Ui of the material point Xj is measured in the material coordinate system, so that a rigid rotation is required to transfer it into the fixed coordinate system. Introducing a mode transformation matrix Φi with its related generalized coordinate q(t), the small elastic displacement can be expressed in the following summation form: 2 3 q1 6

6 q2 7 7 (12.147b) Uj 5 Φj q 5 φj1 φj2 ? φjN 6 7; 4 ^ 5 qN where φjα ðXÞ and qα(t) represent, respectively, the αth elastic mode shape of the dry structure and the corresponding generalized coordinate, and N denotes the retained mode number. These elastic modes satisfy the orthogonal relationships ( ð 0; α 6¼ β T ρs φiα φiβ 5 mαβ 5 ; m 5 diagðmα Þ; mα ; α 5 β Ωs ( ð 0; α 6¼ β (12.147c) ϕTiα;j Cijkl φkβ;l dΩ 5 kαβ 5 ; k 5 diagðkα Þ; k α5β α; Ωs kα m21 k 5 diagðω~ 2α Þ; ω~ 2α 5 : mα Here the diagonal matrices m and k, respectively, are the generalized mass and stiffness matrices of the elastic structure, and ω~ 2α is the square of the αth natural frequency of the structure. The elastic mode shape and the corresponding generalized coordinate are the functions of material coordinate X and time t, respectively. These mode shape functions can be obtained by traditional FEA or theoretical mode functions of the dry structure as previously mentioned. Based on Eq. (12.147b), the elastic structure with an infinite number of DOF is now approximated by a discrete system of DOF N, and the integrated body is now modeled by a system of DOF N 1 6 consisting of 6 rigid DOF and N elastic DOF for numerical simulation. For a free
free elastic structure with no any external forces applied on it, when it vibrates in an elastic natural mode φjα , its inertial force must be in dynamic equilibrium at each time, which implies ð ð 2 ρs φiα q€α dΩ 5 2 q€α ρs φiα dΩ 5 0; Ωs ð Ωs (12.147d) ρs φiα dΩ 5 0: Ωs

Substituting Eq. (12.147b) into Eq. (12.147a), we obtain the displacement and velocity field of the body in the following mode summation form ui 5 xi 2 Xi 5 Uio 1 Rij Xj 1 Rij Φj q; o Vi 5 U_ i 1 R_ij Xj 1 R_ij Φj q 1 Rij Φj q_ o 5 U_ i

(12.147e)

_ 1 eimk θ_ m Rkj ðXj 1 Φj qÞ 1 Rij Φj q:

Eq. (12.147e) can be expressed in the matrix form 2 03 2 03 _

U

U V 5 I A RΦ 4 θ_ 5; u 5 I R RΦ 4 X 5; q q_

U 5 U1

U2

U3

T

5 Φq;

(12.147f)

540

Fluid
Solid Interaction Dynamics

where u 5 u1

u2

X 5 X1

X2

Φ 5 ΦT1 2 6 6 A56 4

T

u3

X3

ΦT2

U0 5 U10

;

T

ΦT3

U20

θ 5 θ1

;

T

θ2

I 5 IT1

;

0

I3 RðX 1 ΦqÞ

0

0

I2 RðX 1 ΦqÞ 0

I1 5 1 0 0 ; I2 5 0 1

U30 θ3

T

T

;

;

IT3 ;

IT2

0

3

(12.147g)

7 7 I1 RðX 1 ΦqÞ 7; 5 0



0 ; I3 5 0 0 1 :

Strain field and strain energy density The linear strain tensor of the small deformation of the elastic body can be calculated by Eq. (4.50), that is, Eij 5

1 1 ðUi;j 1 Uj;i Þ 5 ðΦi;j 1 Φj;i Þq; 2 2

(12.148a)

from which the elastic strain energy density of the body is given by e5

1 1 1 Eij Cijkl Ekl 5 Ui;j Cijkl Uk;l 5 qT ΦTi;j Cijkl Φk;l q; 2 2 2

(12.148b)

where Cijkl is the elastic tensor as given in Section 3.4, and its symmetry has been used to derive Eq. (12.148b). The total strain energy of the elastic body is given by ð 1 1 ~ ΠðqÞ 5 edΩ 5 qT kq 5 q~ T Kq; 2 2 Ωs 2 3 0 0 0

T 6 7 (12.148c) q~ 5 U0T θT qT ; K 5 4 0 0 0 5; 0 0 k ð k 5 ΦTi;j Cijkl Φk;l dΩ 5 kT : Ωs

Here k is the generalized elastic stiffness matrix of the body as defined by Eq. (12.147c). Kinetic energy density Based on Eq. (12.147f), the kinetic energy density of the body per unit volume is calculated by k5

5

1 ρ VT V 2 s 1h 2

T θ_

_ 0T U 2

I

1 T6 5 q_~ 4 AT 2 ΦT RT

2 3 I i 6 T 7 q_ T 4 A 5 I ΦT R T A AT A ΦT R T A

RΦ

A 3

7 AT RΦ 5q_~ ΦT Φ

2 03 _ U

6 7 RΦ 4 θ_ 5 q_

;

(12.149a)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

541

where the rotation matrix satisfying RT R 5 I has been used. The total kinetic energy of the body is given by ð 1 T _ 0 _ _ _ 5 kdΩ 5 q_~ Mq; ~ TðU ; θ; θ; q; qÞ 2 Ωs 2 3 (12.149b) ð I A RΦ T T T 4 5 A Mðθ; qÞ 5 ρs A A A RΦ dΩ; Ωs ΦT RT ΦT RT A ΦT Φ where the matrix Mðθ;qÞ is the mass matrix. Based on Eqs. (12.146c), (12.147c), and (12.147d), we obtain ð ð ρs AdΩ 5 0; ρs AT dΩ 5 0; Ωs Ωs ð ð ρs RΦdΩ 5 0; ρs ΦT RT dΩ 5 0; Ωs Ωs ð ð ~ I Þ; ~ ρs AT AdΩ 5 diagðA ρs AT RΦdΩ 5 B; Ωs Ωs ð ~ 1 5 ρ ðqT ΦT 1 XT ÞRT IT I2 RðX 1 ΦqÞdΩ; A s 2 Ω ð s ~ 2 5 ρ ðqT ΦT 1 XT ÞRT IT I3 RðX 1 ΦqÞdΩ; A s 3 Ω ð s ~ 3 5 ρ ðqT ΦT 1 XT ÞRT IT I1 RðX 1 ΦqÞdΩ; A s 1 Ωs ð ρs ΦT ΦdΩ 5 diagðmα Þ 5 m;

(12.149c)

Ωs

so that the mass matrix Mðθ;qÞ is given by

2

MI

6 Mðθ; qÞ 5 4 0 0

0

3 0 ~ 7 B 5 5 MT ðθ; qÞ:

~ IÞ diagðA T ~ B m

(12.149d)

Generalized force We assume that the solid body is subjected a traction force T^ i on its surface ST and that a traction force caused by the water pressure p, that is, fip 5 2 pηi , on the wet surface Σ. The components of these two traction forces are defined in the material coordinate system, so that their components relative to the Eulerian system are obtained by the corresponding rotation, that is,

^ ^ 5 T^ 1 T^ 2 T^ 3 T ; T i 5 Rij T^ j ; T 5 RT; T (12.150a) p p f i 5 2 pRij ηj ; f 5 2 pRη: The body force of the solid is the gravitational force defined in the Eulerian system, which is represented by

(12.150b) f g 5 0 0 2ρs g : To obtain the generalized force vector Q, we use the virtual work δW produced by all forces when the system is subjected to a virtual displacement 2 03 δU

6 7 δu 5 vδt 5 I A RΦ 4 δθ 5; (12.150c) δq

542

Fluid
Solid Interaction Dynamics

to obtain the virtual work

ð δW 5

ð δu f dΩ 1 T g

Ωs

5 δU0T 5 δU0T

δθT

ð

T p

δu f dS 1 δuT TdS ST ð  ð ð

T g T p T T H f dΩ 1 H f dS 1 H TdS δq Σ

Ωs

Σ

δqT Q 5 δq~ T Q;

δθT

HT 5 I

and the generalized force vector

ST

A

RΦ

(12.150d)

T

ð

Qðp; θ; qÞ 5 F 1 F 1 F ; F 5 HT TdS; ST ð ð ð g p T g T p F 5 H f dΩ; F 5 H f dS 5 2 pHT RηdS: g

p

T

Ωs

T

Σ

(12.150e)

Σ

This generalized force vector is a function of the fluid pressure, rotation, and elastic deformation of the structure. Dynamic equation for numerical simulation When the potential energy in Eq. (12.148c), the kinetic energy in Eq. (12.149c), and the generalized force in Eq. (12.150e) are substituted into the Lagrangian
Hamilton equation (Meirovitch, 1997),       d @T @T @Π 2 1 5 Q; (12.151a) dt @q_~ @q~ @q~ ~ as given in Section 2.2, to obtain the following we can use the derivative rule of a scalar with respect to vector q, dynamic equation in the numerical simulation form: _~ q_~ 1 Kq~ 5 Qðp; θ; qÞ; Mðθ; qÞq€~ 1 Cðθ; q; qÞ  dMðθ; qÞ  T @Mðθ; qÞ C θ; q; q_~ 5 2 q_~ dt 2@q~ 2 3 @Mðθ; qÞ T @Mðθ; qÞ 5 5 q_~ T @Mðθ; qÞ 2 5 q~_ 4 @q~ 2@q~ 2@q~ 5

(12.151b)

N 16 X @M q_~I : @q~I I51

Obviously, Eq. (12.151a) is a nonlinear matrix equation with the DOF N 1 6 consisting of three rigid translations of the mass center, three rigid rotations about the mass center, as well as N elastic degrees describing the elastic deformation of the structure. The mass matrix involving the rigid rotation and the elastic deformation and the damping matrix are affected not only by the geometrical position of the structure but also by its velocity. The generalized force is a function of the geometrical position and the fluid dynamic pressure, which is caused by the FSI. Fluid domain Navier
Stokes equation The water is considered an incompressible perfect fluid, so that its divergence of velocity vanishes and its motion is governed by Eq. (3.135) with the body gravity force fi 5 2 gδi3 in the coordinate system defined in Fig. 12.19: p;i in Ωf ; v_i 5 2 gδi3 2 ; vi;i 5 0; (12.153a) ρf or in the matrix form

Dv rp 5g2 ; Dt ρf

rUv 5 0;

where v and g denote the velocity vector and the body force vector as follows:

T

T g 5 0 0 2g : v 5 v1 v2 v3 ;

(12.153b)

(12.153c)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

543

Based on the MPS method given in Section 12.2.6.2, the fluid velocity, the gradient, and the Laplacian of the pressure with their values at position rn of the time step n can be determined by Eqs. (12.134a), (12.135c), and (12.136f), respectively, that is, vn ðxI Þ 5

N X vn ðxJ ÞWJI

nI

J6¼I

;

 n  min N pJ 2 pn; ðrJ 2 rI Þ dX I ðrp ÞI 5 0 WJI ; n J6¼I jrJ 2rI j2 n

ðr2 pn ÞI 5

N X  2d  n pJ 2 pnI WJI : 0 n λ J6¼I

(12.153d)

(12.153e)

(12.153f)

Now, using the projection method discussed in Section 12.2.5.4, the intermediate velocity marked by *, without considering the pressure gradient contribution, and the corresponding intermediate position of particle can be calculated as v 5 vn 1 Δtg;

r 5 rn 1 Δtv :

(12.154a)

To include the pressure contribution, the velocity at time step n 1 1 is calculated by vn11 5 v 2 Δt

rpn11 ; ρf

(12.154b)

from which, when a diligence operation is taken and considering the incompressibility condition in Eq. (12.153b) valid at time n 1 1, we obtain a pressure Poisson equation: r2 pn11 5

ρf rUv ; Δt

(12.154c)

which has its particle form in Eq. (12.134h), that is, 

ðr p

2 n11

ρ f n 2 n0 ÞI 5 2 2 I 0 : Δt n

(12.154d)

For pressure calculation accuracy, this equation is further modified to include the numerically accumulated error in the mass density at time step n (Sun et al., 2015a,b; Sun, 2016), so that the final Poisson equation is given by 0 1  0 n 0 ρ n 2 n n 2 n f ðr2 pn11 ÞI 5 2 2 @ I 0 1 α I 0 A Δt n n 2 3 (12.154e)  ρf ð1 1 αÞn0 2 nI 2 αnnI 5; 5 24 Δt n0 where

 8   0 > n  > n 2 n >  1 ΔtjrUvn j; ðn0 2 nn ÞrUvn $ 0;  > > > <  n0   α 5   > 0 n >  n 2 n > > ;  > ðn0 2 nn ÞrUvn , 0: > :  n0 

(12.154f)

Physically, ðn0 2 nn ÞrUvn . 0 implies that both (n0
nn) and rUvn are negative or positive. The negative case implies the fluid is compressed at time n and is further compressed due to the motion with rUvn , 0, while the positive case means the reverse: the fluid expands at time n, and further expands by the flow of rUvn . 0; therefore an additional term ΔtjrUvn j is added in the Poisson equation. However, when ðn0 2 nn ÞrUvn , 0, the particle density (n0
nn) and the velocity divergence rUvn have different developing directions that amend each other, and no term ΔtjrUvn j is added in the modified parameter α.

544

Fluid
Solid Interaction Dynamics

Combining Eqs. (12.153f) and (12.154e), we obtain   N X  ρf nI 1 αnnI 2 ð1 1 αÞn0 2d  n11 n11 W p 2 p 5 ; JI J n0 λ I Δt2 n0 J6¼I which is the equation to solve the pressure expressed in the following matrix form: 2 3n 2 3n11 2 3n p1 N1 H11 ? H1N & ^ 5 4 ^ 5 5 4 ^ 5 5 N; HP 5 4 ^ HN1 ? HNN pN NN

(X

where

N

L6¼I

HIJ 5

2d WLI ; ðn0 λÞI

(12.155a)

(12.155b)

J 5 I;

2 2d WJI ; J 6¼ I; ðn0 λÞI "  # ρf nI 1 αnnI 2 ð1 1 αÞn0 NI 5 2 : Δt n0

(12.155c)

This equation is a simultaneous equation expressed by a linear symmetric coefficient matrix H, of which the solution gives the pressure at time n 1 1. When the pressure is obtained, the velocity and the location of the fluid can be updated by vn11 5 r

n11

vn 2 Δtrpn11 ; ρf

5 r 1 Δtv n

n11

(12.156)

:

Fluid boundary conditions For the fluid boundaries without FSI, such as free surface and fixed solid boundaries, the treatment approaches described in Section 12.2.5.5 can be used. For the FSI boundary, the motion of the solid body is also not known, so a special consideration is required. Fluid
solid interaction boundary condition Pressure Neumann condition On the FSI boundary, both a force equilibrium and a geometrical consistency must be satisfied. The force equilibrium condition has been considered by the solid dynamic equation in Eq. (12.151b), where the generalized force Qðp; θ; qÞ is a function of the fluid pressure p. Therefore we only need to consider the geometrical consistency condition given in Eq. (4.16a), which is rewritten in the reference system defined in Fig. 12.19, where the outside normal vector of fluid ν 5 2 η (for the solid): _ νn11 Urpn11 5 ρf ðνn11 Ug 2 νn11 UV

n11

Þ:

(12.157a)

_ n11 and the unit normal vector νn11 of the FSI interface are also In this equation, both the solid acceleration V unknown at time step n 1 1, so that an iteration process is needed. As an approximation, their values at time n may be _ n11 can be derived by using Eq. (12.147f) while the used in the corresponding explicit form. The solid acceleration V gradient of the pressure can be expressed in particle form as given by Eq. (12.135.6) in numerical calculations. Laplacian operator compensation near fluid
solid interaction interface For the fluid particles near the FSI interface, the Laplacian operator is modified to be consistent with the Newman condition in Eq. (12.157a). As shown in Fig. 12.20, a virtual particle in the solid is positioned along the normal direction ν and away from the FSI interface by a distance dr0. This virtual particle is also included in the discretization equation of the Laplacian operator. Based on Eq. (12.157a), the pressure of the fluid pF is derived as _ 0; pF 5 pS 1 ρf ðνUg 2 νUVÞdr where pF denotes the pressure of the solid particle.

(12.157b)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

re

545

Virtual particle

ν dr0

r

Fluid particle under concern dr0

Boundary particle

FIGURE 12.20 Pressure modification on FSI interface (Sun, 2016). FSI, Fluid
solid interaction.

Intermediate velocity of fluid
solid interaction boundary The intermediate velocity v obtained from Eq. (12.154b) should guarantee the velocity consistency condition on the FSI interface, as suggested by Brown et al. (2001). Eq. (12.157a) is just the geometrical condition in the pressure gradient form, from which, when it combines with Eq. (12.154b), it follows _ νUv 5 νU½vn11 1 Δtðg 2 V

n11

Þ;

(12.158a)

from which, when the fluid velocity equals the solid velocity, that is, vn11 5 Vn11 , it follows _ νUv 5 νU½Vn11 1 Δtðg 2 V

n11

Þ;

(12.158b)

_ which may be used to modify the intermediate velocity v . Actually, if we approximate ΔtV Eq. (12.158b) further reduces to

n11

νUv 5 νUðΔtg 1 Vn Þ:

5 Vn11 2 Vn , (12.158c)

This equation is just the projection of Eq. (12.154a) on the unit normal direction. Therefore if the fluid velocity exactly satisfies vn11 5 Vn11 on the FSI boundary, the modification by Eq. (12.158b) seems unnecessary. However, in numerical simulations, the condition vn11 5 Vn11 can never be exactly satisfied due to the numerical error at each time step, so that the modification in Eq. (12.158b) is needed. For a perfect fluid, along the tangent direction of an FSI interface, the fluid velocity and the solid velocity can be different, and there is no geometrical condition in the tangent direction. For the viscous fluid, the nonslip condition requires the fluid and solid to have the same velocity on the FSI interface; as a result, Eq. (12.158b) should be modified to _ n11 Þ: v 5 Vn11 1 Δtðg 2 V

12.4.2.2 Integrated coupling equations Here we list the integrated coupling equations used in numerical simulations.

(12.158d)

546

Fluid
Solid Interaction Dynamics

Solid domain Eqs. (12.147a), (12.147f), and (12.151b) should be satisfied at time step n 1 1; they are rewritten as _~ n11 q_~ n11 1 Kq~ n11 5 Qðp; θ; qÞn11 ; 1 Cðθ; q; qÞ 2 0 3n11 U

6 7 n11 n11 n11 u 5 I R R Φ 4X5 ;

Mðθ; qÞn11 q€~

n11

Vn11 5 I An11

q 2 0 3n11 _ U

6 7 n11 R Φ 4 θ_ 5 ; q_

(12.159a)

(12.159b)

Un11 5 Φqn11 ; xn11 5 X 1 U0n11 1 Rn11 X 1 Rn11 Un11 ; in which, the elastic mode matrix Φ, the elastic stiffness matrix K, and the material coordinate vector X depend on only the geometrical and physical properties of the body and do not change with time. Fluid domain Eqs. (12.154a), (12.155b), (12.55c), (12.153e), and (12.156) are satisfied at time n 1 1 as follows: v 5 vn 1 Δtg; r 5 rn 1 Δtv ; 2 3n 2 3n11 2 3n p1 H11 ? H1N N1 & ^ 5 4 ^ 5 5 4 ^ 5 5 Nn ; Hn Pn11 5 4 ^ HN1 ? HNN pN NN

(X N

HIJn

L6¼I

5

2d WLI ; ðn0 λÞI

ðrp

 n11  min N pJ 2 pn11; ðrJ 2 rI Þ dX I ÞI 5 0 WJI ; 2 n J6¼I jrJ 2rI j vn11 5 r

n11

vn 2 Δtrpn11 ; ρf

5 r 1 Δtv n

n11

(12.160b)

J 5 I;

2 2d WJI ; J 6¼ I; ðn0 λÞI 2 3  n 0 ρ n 1 αnI 2 ð1 1 αÞn 5 f ; NIn 5 2 4 I Δt n0

n11

(12.160a)

(12.160c)

(12.161a)

(12.161b)

:

Fluid
solid interaction interface The coupling conditions in Eqs. (12.157a) and (12.158b) are rewritten as

 _ n11 ; νn11 Urpn11 5 ρf νn11 Ug 2 νn11 UV h

i _ n11 : νn11 Uv 5 νn11 U Vn11 1 Δt g 2 V Initial conditions At the initial time t 5 0, the initial position and the velocity of the solid and fluid are given by

(12.162a) (12.162b)

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

547

Solid body Initial position: Mass center: Material point:

x0o 5 0; x0 5 X:

(12.163a)

Initial displacement and velocity: Mass center: Rigidrotation: Elastic deformation:

_ ; U0 5 0 5 U _ θ 5 0 5 θ; _ U 5 0 5 U: 0

(12.163b)

Fluid domain Initial position, velocity and pressure: Position vector: r 5 r0 ; Velocity: v 5 v0 ; Pressure: p 5 2 xi δ3i g:

(12.164)

12.4.2.3 Numerical simulation process As discussed in Chapter 11, Mixed finite element
computational fluid dynamics method for nonlinear fluid
solid interactions, for nonlinear FSI problems, the partitioned iteration solution is suggested. The numerical process follows the flowchart given in Fig. 12.18 and entails three main steps. 

1. Starting from the initial position of the system (n 5 0), in which the position u~ n11 5 un with the corresponding velocity and acceleration on the FSI interface prescribed, use the fluid equations in Eqs. (12.160a)
(12.160c), (12.161a), (12.161b), (12.162a), and (12.162b) to solve a fluid problem with the prescribed initial fluid pressure and  velocity fields, as well as its boundary conditions in order to obtain an intermediate pressure P ðt 1 ΔtÞ 5 p n11 ,   velocity v n11 , and position vector r n11 .  2. The resultant pressure p n11 is substituted into the solid Eqs. (12.159a) and (12.159b), and the solid equations are    solved to obtain the solid displacement u n11 , velocity u_ n11 , and acceleration u~ n11 using the time integration approaches.   3. Compare the displacement u n11 and the testing displacement u~ n11 by checking whether the convergence  obtained   condition u n11 2 u~ n11  , ε is reached. If it is reached, go to next time step; otherwise modify     u~ n11 5 u~ n11 1 αðu n11 2 u~ n11 Þ to repeat the preceding process. The relaxation factor α can be updated by Eq. (12.146a).

12.4.3

Application examples

Here, based on the numerical methods formulated in this chapter, some application examples of nonlinear FSI problems completed by our PhD students, colleagues, and the author (see Sun et al., 2011, 2012, 2013, 2014, 2015a,b; 2016, 2017; Sun, 2013, 2016; Javid et al., 2013a,b, 2014a,b, 2016; Javed, 2015) are presented to demonstrate the proposed numerical methods.

12.4.3.1 Rigid wedge dropping on the water Water impact is an important problem in marine and offshore engineering (Khabakhpasheva and Korobkin, 2003, 2013). Wedge dropping tests have been used to study the loads of a ship slamming into water to determine safety considerations in the design of marine structures. As the velocity of the dropping wedge depends on FSI interactions, simulations could be difficult if grid-based methods are employed to treat breaking waves on the free surface and moving FSI boundaries. Oger et al. (2006) applied the WCSPH method to wedge dropping simulation using denser particles in the impact area, and the smoothing length was changed depending on the requirement of accuracy to ensure an acceptable level of density fluctuation in the fluid. Gong et al. (2009) proposed an alternative method by using a sponge layer on the bottom of the tank to adjust the density calculation of the fluid particles. When Shao’s ISPH method was applied to water entry of a free-falling

548

Fluid
Solid Interaction Dynamics

wedge, mirror particles were used on the moving solid (Shao, 2009). With Lee’s ISPH method, the proposed two boundary treatments can be applied, which simplifies the model generation and reduces computation time. In this example, simulated by Sun et al. (2011) and Sun (2013), as shown in Fig. 12.21, a symmetric wedge with a dead-rise angle θ 5 π=6 (an angle measured upward from a horizontal plane at keel level) dropping into water is simulated using ISPH with the denser wall particle boundary treatment. In the simulation, the particle spacing is 0.01 m, but near the boundary it is 0.005 m, and the time step is 0.001 second. The weight of the dropping wedge is 241 kg, the width of the wedge is 0.5 m, and the length is 1 m. The tank size is 2 m 3 1 m. In the simulation, the wedge is placed just above the free surface of the calm water with a dropping velocity of 6.15 m/s given from the 2D experiment done by Zhao et al. (1997). The wedge is allowed to move only in the vertical direction, and its motion follows the equation of motion for rigid body. The resultant water pattern generated in the simulation is shown in Fig. 12.22, which is compared with the photo from the experiment by Tveitnes et al. (2008). When the falling wedge hits the water surface, the surface breaks because of the strong impact, and the water is pushed up around the wedge. The breaking wave pattern obtained from the SPH simulation is marked with a black line for easy comparison with the experiment. The velocity of the wedge

1m 0.5 m

θ

1m

2m

FIGURE 12.21 Sketch of a rigid wedge dropping into the water (Sun, 2013).

FIGURE 12.22 Breaking wave pattern compared with the experiment photo taken by Tveitnes et al. (2008) and Sun (2013).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

549

6.50 Experimental (Zhao et.al., 1997)

6.00

Vy(m/s)

5.50 5.00 4.50 4.00 3.50 3.00 0.00

0.01

0.02

0.03

Time (s)

FIGURE 12.23 Wedge dropping velocity compared with the experiment result given by Zhao et al. (1997) and Sun (2013).

Experimental (Zhao et.al., 1997) Analytical

7000.0 6000.0

Fy/N

5000.0 4000.0 3000.0 2000.0 1000.0 0.0 0.00

0.01

0.02

0.03

Time (s)

FIGURE 12.24 Impact force on the wedge compared with the experiment/analytical results given by Zhao et al. (1997) and Sun (2013).

after impact and the impact forces on the wedge from the water are compared with the experimental data and analytical results given by Zhao et al. (1997) in Figs. 12.23 and 12.24, respectively. Fig. 12.23 shows good agreement between the velocity obtained by the current SPH method and the experiment result, and it is slightly lower than the experimental values in the later stage with a maximum error of 2%. Also, it is seen that the wedge dropping velocity decreases more rapidly after 0.007 second since the water buoyancy produces an upward acceleration after the wedge enters the water. As observed in Fig. 12.24, the fluid force initially increases steadily and then slows down before reaching the peak at around 0.015 second, after which it starts to decrease. The computed force values are slightly overpredicated at first, and then it is underpredicated for a short period of time, but it is higher than experimental values at the last stage. Overall, both the dropping velocity and the vertical fluid force obtained from the proposed SPH method agree well with experimental values. To investigate the effect of different parameters on the water entry process and to provide some guidelines for engineering applications, the following three cases are studied by Sun (2013): 1. Case 1: Different wedge masses with a dead-rise angle of 30 degrees and initial dropping velocity of 6.15 m/s 2. Case 2: Different initial dropping velocities with a wedge mass of 241 kg and dead-rise angle of 30 degrees 3. Case 3: Different dead-rise angles with wedge mass of 241 kg and initial dropping velocity of 6.15 m/s Fig. 12.25 shows the results affected by different wedge masses, from which it is clear that, as the wedge mass increases, both the wedge velocity and the water force on the wedge are increased. This may be explained by the momentum theorem. Assume that during the impact time period Δt, the wedge of mass M is subject to the fluid impact force F, and its velocity reduces from the initial value V0 to the ending value Vt. The momentum theorem for this impact process is given by MðVt 2 V0 Þ 5 ðMg 2 FÞΔt; from which it follows

 F Vt 5 V0 1 g 2 Δt; M

(12.165a)



(12.165b)

Vertical velocity (m/s)

550

Fluid
Solid Interaction Dynamics

FIGURE 12.25 Time history of wedge velocities (top) and fluid forces (bottom) affected by different wedge masses (Sun, 2013).

M=141 kg M=241 kg M=341 kg

6.5 6 5.5 5 4.5 4 3.5 3 0

0.005

0.01

0.015

0.02

0.025

0.03

Time (s)

Fluid force (N)

7000

M=141 kg M=241 kg M=341 kg M=441 kg

6000 5000 4000 3000 2000 1000 0 0

0.005

0.01

0.015

0.02

0.025

0.03

Time (s)

Vertical velocity (m/s)

6.5

FIGURE 12.26 Time history of dropping velocities (top) and fluid forces (bottom) affected by different initial wedge velocities (Sun, 2013).

V0=6.15 m/s V0=5.15 m/s V0=4.15 m/s

6 5.5 5 4.5 4 3.5 3 0

0.005

0.01

0.015 Time (s)

0.02

Fluid force(N)

6000

0.025

0.03

V0=6.15 m/s V0=5.15 m/s V0=4.15 m/s

5000 4000 3000 2000 1000 0 0

0.005

0.01

0.015

0.02

0.025

0.03

Time (s)

  V 0 2 Vt F5M g1 5 Mðg 1 at Þ: Δt

(12.165c)

Eq. (12.165b) indicates that for the same impulse FΔt, the wedge velocity Vt increases as the mass increases, while Eq. (12.165c) shows that for the same velocity change rate at, the impact force also increases with a larger mass. Fig. 12.26 shows that the decreasing rate of the wedge velocity increases with increasing initial dropping velocity before time 0.0125 but that, after this time, the three curves tend to be parallel with almost the same rate. For the fluid force, larger initial velocity generates higher fluid force at the early stage with a larger peak value. By contrast, the force values become almost the same at the later stage, which is tending to the static buoyance with the wedge velocity further reduced.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

FIGURE 12.27 Time history of dropping velocities (top) and fluid forces (bottom) affected by different dead-rise angles (Sun, 2013).

6.5 Vertical velocity (m/s)

551

6 5.5 5 4.5

Dead-rise angle=30º

4

Dead-rise angle=45º

3.5

Dead-rise angle=60º

3 0

0.005

0.01

0.015 Time (s)

0.02

Fluid force (N)

6000

0.025

0.03

Dead-rise angle=30º Dead-rise angle=45º

5000

Dead-rise angle=60º 4000 3000 2000 1000 0 0

0.005

0.01

0.015 Time (s)

0.02

0.025

0.03

Fig. 12.27 shows that fluid force acting on a wedge is smaller for the wedge with a larger dead-rise angle and that consequently the velocity decreases much more slowly. This is because the vertical force on the wedge is a projection of fluid pressure on the vertical line, which is proportional to the cosine of the dead-rise angle, so that a larger dead-rise angle leads to a smaller force.

12.4.3.2 Flexible wedge dropping on the water Further simulation of a wedge dropping problem including wedge elastic deformation was completed by Sun et al. (2015a) and Sun (2016). Fig. 12.28 shows the initial configuration of the problem, of which the tree different flexible cases with the corresponding parameters listed in Table 12.1 are considered. In the simulations, the symmetrical wedge is assumed to be constructed by a rigid frame marked by shading lines and the two inclined elastic beams on the bottom surface. The integrated system is considered a 2D symmetrical one with its symmetrical motions. The mass center of the wedge undergoes a large rigid motion, and its bottom elastic beams move with the mass center and undergo their small elastic deformation modeled by the mode summation method. For each beam, three theoretical modes are taken into account, and the corresponding first three circular frequencies are: 96.210, 602.943, and 1688.238. The MPS method described in Section 12.2.6.2 is adopted in modeling the fluid motion. The fluid field is discretized by particles with the initial space of 0.005 m, which results in 38,400 fluid particles and a total 40,122 particles including solid particles on the FSI interface. The time interval is determined by the Courant
Friedrichs
Lewy condition with a maximum limit of 0.0002 second. It is worth mentioning that the majority of the computational time is used for the fluid solver, that is, MPS part. The time used for the structure solver is negligible since the number of DOF of the solid is only 4, consisting of a rigid and three elastic ones. The convergence error is set to 1025 using the averaged coupling root mean square (RMS) (CRMS) as the error criterion. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u NFSI u 1 X   (12.166) ðu n11 2 u~ n11 ÞI Uðu n11 2 u~ n11 ÞI ; CRMS 5 t NFSI I51 



where NFSI denotes the particle number on the FSI interface, and ðu n11 2 u~ n11 ÞI denotes the difference of the coupling var  iable values of particle I at time step n 1 1 between the starting value u~ n11 and the ending value u n11 in the iterations. The accelerations of the flexible wedge, calculated by the modified MPS, are compared with the results from experimental, rigid body simulation and Wagner’s theory (Panciroli et al., 2012, 2013; Panciroli, 2013, respectively) in

552

Fluid
Solid Interaction Dynamics

FIGURE 12.28 Initial configuration of flexible wedge dropping on the water (Sun, 2016).

Y 0.8 m

Equivalent dropping height H*

0.8 m 0.6 m

X

0 1.6 m

TABLE 12.1 The parameters of three cases of flexible wedge. Case 1 Material

Case 2

E-glass (woven)/epoxy

Case 3 Aluminum

Young modulus, E (GPa)

30.3

68

Density, ρ (kg/m )

2015

2700

3

Mass of the rig (kg/m)

22

Length of each bottom, L (m)

0.3

Thickness (mm)

2.2

Dead-rise angle, β (degree)

2

30

Entry speed (m/s)

4.29

5.57

4

Equivalent height, H* (m)

0.938

1.5813

0.8155

Fig. 12.29. As shown in Fig. 12.29A and B for elastic wedges, the numerical results coincide with the experimental data for the main trend; the first impact pressure peak as well as the oscillating tend to a constant value with time until the end of the simulation. It can be seen that for the two elastic cases there is a trough at about 0.025 second on both the experimental and the numerical curves, which is not shown on the curves for the rigid cases given in Fig. 12.29C and D. The numerical curves show an overshot of the second acceleration peak at the time earlier than the experimental curve. This is probably caused by a 2D beam approximation in the simulations for the 3D experiment; the real 3D bottom is a plate with different natural frequencies from those of the beam. The pressure and velocity contours are shown in Figs. 12.30 and 12.31, respectively. Due to the flexibility of the wedge bottom, the cavitation starts to develop from about t 5 0.02 second and vanishes until about t 5 0.04 second. Because the current model involves only the water phase, the dynamics caused by the entrapped air between the wedge bottom and water could not be captured correctly. Fig. 12.32 shows the deformation of the flexible bottoms at some typical time instants. During the initial stage of the impact, the bottoms are bended down by the coupling effect of the inertia and the concentrated impact force near the wedge tip. After about t 5 0.02 second, the deformation of the beam starts to bounce back toward the symmetry line of the wedge.

(A) 250

(C) 250 Rigid - Modified MPS+rigid-body dynamics Rigid - Wagner theory

Flexible - Modified MPS+CRMS Flexible - Experiment

200

150

acc (m/s2)

acc (m/s2)

200

100

150 100

50 50 0 0

0.01

0.02

0.03

0.04

0.05

Time (s)

(B)

0

0.06

0.02

0.03

0.04

0.05

0.06

Time (s) 400

Flexible - Modified MPS+CRMS Flexible - Experiment

350

Rigid - Modified MPS+rigid-body dynamics Rigid - Wagner theory

350

300

300

250

acc (m/s2)

acc (m/s2)

0.01

(D)

400

200 150 100

250 200 150

50

100

0

50

–50

0

0

0.01

0.02

0.03

Time (s)

0.04

0.05

0.06

0

0

0.01

0.02

0.03

0.04

0.05

0.06

Time (s)

FIGURE 12.29 Acceleration of the flexible/rigid wedge obtained by MPS and compared with the experiment or Wagner’s theory (Panciroli et al., 2012, 2013; Panciroli, 2013): (A) elastic case 1, (B) elastic case 2, (C) rigid with case 1 entry speed 4.29 m/s, (D) rigid with case 2 entry speed 5.57 m/s (Sun, 2016). MPS, Moving particle semiimplicit.

0.6

0.6

0.6

0.5

0.5

0.5

0.4

y (m)

0.7

0.4 0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0 0.3

0.4

0.5

0.6

0.7

0.8 x (m)

0.9

1

1.1

1.2

1.3

0 0.3

0.4

0.5

0.6

0.7

0.8

t=0.02 s

0.8 x (m)

0.9

1

1.1

1.2

1.3

0.6

0.6

0.5

0.5

0.5

y (m)

0.6

y (m)

0.7

0.4 0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

0 0.3

0.4

0.5

0.6

0.7

0.9

1

1.1

1.2

1.3

0.8

t=0.035 s

0.6

0.5

0.5

y (m)

0.6

0.5

y (m)

0.6

0.4

0.4 0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0.5

0.6

0.7

0.8 x (m)

0.9

1

1.1

1.2

1.3

0

0.3

0.4

0.5

0.6

0.7

0.8 x (m)

FIGURE 12.30 Water pressure contours at different times obtained by case 2 simulation (Sun, 2016).

0.9

1

1.1

1.2

1.3

t=0.03 s

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

0.9

1

1.1

1.2

1.3

1

1.1

1.2

1.3

t=0.06 s

0.4

0.3

0.4

0.8

0.8

t=0.04 s

0.7

0.3

0.7

x (m)

0.7

0

0.6

x (m)

0.8 0.7

0.8

0.5

0.4

0.3

0.4

0.4

0.8

t=0.025 s

0.7

0.3

0.3

x (m)

0.7

0.4

t=0.015 s

0.4

0.3

x (m)

y (m)

0.8

t=0.01 s

0.7

0.8

y (m)

0.8

t=0.005 s

0.7

y (m)

y (m)

0.8

0

0.3

0.4

0.5

0.6

0.7

0.8 x (m)

0.9

0.8

0.8

y (m)

y (m)

0.6

0.5

0.4

0.3

0.7

0.6

0.6

0.5

0.4

0.5

0.6

0.7

0.8 x (m)

0.9

1

1.1

1.2

0.4

0.3

0.3

1.3

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

0.2

1.3

0.8

0.6

0.5

y (m)

0.6

y (m)

0.6

0.5

0.4

0.4

0.3

0.3

0.3

0.6

0.7

0.8 x (m)

0.7

0.9

1

1.1

1.2

0.2

1.3

0.8

0.2 0.3

0.4

0.5

0.6

0.7

0.8

t = 0.035 s 0.7

0.7

0.6

0.6

0.8 x (m)

0.9

1

1.1

1.2

1.3

1

1.1

1.2

1.3

1

1.1

1.2

1.3

0.5

0.4

0.5

0.6

t = 0.03 s 0.7

0.4

0.5

t = 0.025 s 0.7

0.3

0.4

0.8

t = 0.02 s 0.7

0.2

0.3

x (m)

0.8

y (m)

0.5

0.4

0.2 0.3

t = 0.015 s

0.7

y (m)

0.7

0.2

0.8

t = 0.01 s

t = 0.005 s

0.8 x (m)

0.9

1

1.1

1.2

1.3

t = 0.04 s

0.3

0.4

0.5

0.6

0.7

0.8

0.8 x (m)

0.9

t = 0.06 s

0.5

y (m)

y (m)

y (m)

0.7

0.5

0.6

0.5 0.4

0.4

0.3

0.3

0.4

0.3 0.2

0.2

0.3

0.4

0.5

0.6

0.7

0.8 x (m)

0.9

1

1.1

1.2

1.3

0.3

0.4

0.5

0.6

0.7

FIGURE 12.31 Velocity contours at different times obtained by case 2 simulation (Sun, 2016).

0.8

x (m)

0.9

1

1.1

1.2

1.3

0.3

0.4

0.5

0.6

0.7

0.8 x (m)

0.9

0.75 0.8

0.75 0.7

0.75

t = 0.01 s

0.7

0.65

0.6

0.6

0.55

0.5

0.7

t = 0.015 s

0.65

y(m)

y(m)

y(m)

t = 0.005 s

0.65

0.6

0.6

0.7

0.8

0.9

1

0.55

0.5

1.1

0.6

0.7

x(m)

0.8

0.9

1

1.1

0.5

0.6

0.7

x(m)

0.8

0.9

1

1.1

0.9

1

1.1

1

1.1

x(m) 0.75

0.75 0.75 0.7 0.7

t = 0.025 s

0.7

0.65

t = 0.03 s

0.65 y(m)

y(m)

y(m)

t = 0.02 s 0.65

0.6

0.6 0.6 0.55 0.55 0.55 0.5 0.5 0.5 0.5

0.6

0.7

0.8 x(m)

0.9

1

1.1

0.5

0.6

0.7

0.8 x(m)

0.9

1

1.1

0.5

0.7

0.65

y(m)

y(m)

y(m)

t = 0.06 s

0.6

0.5

0.5

0.45

0.5

0.5

0.6

0.7

0.8 x(m)

0.9

1

1.1

0.55

0.55 0.55

0.45 0.5

0.8 x(m)

0.6

t = 0.04 s

0.65

0.6

0.7

0.65

0.7

t = 0.035 s

0.6

0.6

0.7

0.8 x(m)

FIGURE 12.32 Beam elastic deformation at different times obtained by case 2 simulation (Sun, 2016).

0.9

1

1.1

0.5

0.6

0.7

0.8 x(m)

0.9

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

(A)

557

200 SPH Modified MPS+CRMS Experiment

0 –200

Strain (x10–6)

–400 –600 –800 –1000 –1200 –1400 –1600 –1800

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

t (s) (B)

800 SPH Modified MPS+CRMS Experiment

600

Strain (x10–6)

400 200 0 –200 –400 –600 –800 –1000

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

t (s) FIGURE 12.33 The beam strains (case 3) at points (A) 30 mm and (B) 120 mm from beam tip by MPS: (A) at 30 mm from beam tip and compared with the ones by SPH and experiment (Panciroli, 2013; Sun, 2016). MPS, Moving particle semiimplicit; SPH, smoothed particle hydrodynamics.

This process is also reflected in the time history of strain that is monitored at two different locations on the upper surface of each bottom beam, that is, 30 and 120 mm away from the wedge tip respectively, as shown in Fig. 12.33. The positive part of the strain in Fig. 12.33B, monitored at 120 mm (about the middle of the bottom), represents the initial downbending stage of the beam shown in Fig. 12.32. After that, the strain remains negative because of the pressure of the surrounding water. The numerical strain record matches well with the experimental data.

12.4.3.3 Two-dimensional water dam breaking impact on rigid/flexible beams Spring-supported rigid beam Fig. 12.34 shows the sketch of a water dam break impact to a rigid uniform beam of length L 5 2 m, total mass M 5 1 kg, and supported by a rotation spring of stiffness K 5 1500 N m at the base. The rotation angle θ of this rigid beam about it base is governed by its dynamic equation derived by the Newton’s second law, that is, 1 I θ€ 1 Kθ 2 MLg sin θ 5 T; 2 T5

NFSI X J51

(12.167) pJ YJ ;

558

Fluid
Solid Interaction Dynamics

FIGURE 12.34 Sketch of a 2D water dam break impact to a spring-supported rigid beam (Sun, 2016). 2D, Twodimensional.

2m

θb

0.6 m

1.2 m

3.22 m

where I 5 ML2/3 is the rotation moment of inertia of the beam about the beam base point. The torque T of fluid pressure is obtained by a summation of the product of each particle pressure pJ and distance YJ of its impact point along the beam to the base point, since the fluid pressure is always perpendicular to the beam line. The geometrical parameters of the initial water dam and the flow domain are indicated in the figure. The initial particle distance is chosen as 0.01 m, which results in 7200 fluid particles and a total 8652 particles, including the solid ones on the FSI interface. The MPS method was used to simulate this problem by Sun (2016). The pressure contour and free surface profiles at selected time instants are shown in Fig. 12.35, from which it can be seen that the distribution of the pressure is quite smooth in the space domain. Two major impacts are found in two durations, t 5 0.75
2 seconds and t 5 5.1
5.9 seconds, during which the large fluid pressure generated from the falling of the water column pushes the beam to relatively large angles. Fig. 12.36A shows the time history of the rotational angle of the beam, on which the two significant angle impulses are consistent with the time of the two major impacts observed in Fig. 12.35. Except for these violent interactions, the water pressure applied on the beam is relatively much smaller. In the video recording of the simulation results, it was seen that the beam oscillates in a period 0.2 second, which is larger than the period 0.1879 second of its natural frequency due to the FSI effect with additional water mass on the beam. In order to investigate the effect of the beam rotation on the pressure field, the time history of the pressure monitored at a beam point 0.16 m above the right corner of the beam is compared in Fig. 12.36B: case 1 for the springsupported rotational beam and case 2 for the nonrotational fixed one. From this figure it can be seen that the time history during the first major impact in case 1 is basically same as the one in case 2 for the rigid fixed beam case. However, the peak pressure value in case 1 during the second major impact is much larger than in case 2. This is because at the beginning of the second impact, the beam rotates back, so that the water front and the beam move toward each other and have a more vigorous impact. After t 5 6 seconds for the rotational beam case, both the rotational angle and the water pressure show an oscillation of the frequency near the natural frequency of the rotational beam, which is not observed for case 2, fixed rigid beam shown in Fig. 12.36B. Fixed elastic beam Fig. 12.37 shows a 2D elastic beam
water dam FSI system with the related geometrical parameters given. The elastic beam is fixed at its base, and its parameters are chosen as: Young’s modulus E 5 0.2 GPa, thickness D 5 0.006 m, mass density m 5 47.16 kg/m, and the moment of inertia of the cross section I 5 1:8 3 1028 m4 . Based on the FE
MPS method proposed for a nonlinear FSI, for the water domain, the initial particle distance is 0.004 m, which results in 1250 fluid particles and a total 1736 particles for the integrated system; the beam’s motion is modeled by the linear FEM. The pressure contours and the beam deformation at some selected time instants are shown in Fig. 12.38. The deflection X of the beam top is shown in Fig. 12.39A, from which it can be seen that the period of its oscillation is about 0.6 second, very close to its first natural period 0.5971 second, with the modal shape shown in Fig. 12.37 by the dashed line The time histories of water pressure monitored at point Y 5 0.02 m of both the rigid and the elastic wall cases are given in Fig. 12.39B, which shows the very small fluctuation of pressure for both cases with the negligible difference between them. This implies that the small elastic deformation of the beam will not greatly change the fluid motion and its pressure field.

2.5

2.5

2.5

2.5

2

2

2

2

t=0.75 s

t=1.25 s

t=1.56 s 1.5

y (m )

y (m )

y (m )

t=2 s

1.5

1.5

y (m )

1.5

1

1

1

1

0.5

0.5

0.5

0.5

0

0.5

1

1.5

2

2.5

3

3.5

0

0

0.5

1

1.5

2 x (m)

x (m) 2.5

2.5

3

3.5

0

4

2.5

2

2

t=2.2 s

0

0.5

1

2.5

3

3.5

0

4

0

2.5

2

2

t=5.1 s

1.5

1

1

1

1

0.5

0.5

0.5

0.5

0

0

0.5

1

1.5

2 x (m)

2.5

3

3.5

0

0 0

0.5

1

1.5

2

2.5

3

x (m)

FIGURE 12.35 Pressure contour and free surface profiles at selected time instants (Sun, 2016).

3.5

0.5

1

1.5

2 x (m)

2.5

3

3.5

4

3

3.5

t=5.9 s

1.5

y (m )

y (m )

y (m )

1.5

2 x (m)

2.5

t=3 s

1.5

1.5

y (m )

0

0

0.5

1

1.5

2 x (m)

2.5

3

3.5

4

0

0

0.5

1

1.5

2 x (m)

2.5

560

Fluid
Solid Interaction Dynamics

(A) 0.4

Arc degree

0.3 0.2 0.1 0 –0.1

0

1

2

3

4

5

6

7

8

9

10

Time (s) (B) 8000 Elastic beam Rigid beam

p (Pa)

6000 4000 2000 0

0

1

2

3

4

5

6

7

8

9

10

Time (s) FIGURE 12.36 Time history of the rotation angle (A) and the pressure at monitor point (B): case 1, elastic beam for spring-supported rotational one; case 2, rigid beam for fixed rigid one (Sun, 2016).

y

η (s)

0.304 m

0.2 m

x 0.1 m

0 0.352 m

FIGURE 12.37 Sketch of water dam break impact to an elastic beam wall (Sun, 2016).

12.4.3.4 Flow-induced vibration of two-dimensional cylinder As reported by Javed (2015) and Javed et al. (2016), a coupled meshfree
mesh-based fluid solver was employed for flow-induced vibration problems. As shown in Fig. 12.40, fluid domain is modeled by a hybrid grid consisting of a body conformal meshfree nodal cloud around the solid object and a static Cartesian grid surrounding the meshfree cloud in the far-field. The meshfree nodal cloud (Fig. 12.41) provides flexibility in dealing with solid motion by moving and

t=0.24 s 0.3

0.3

0.3

0.25

0.25

0.25

0.25

0.2

0.2

0.2

0.2

0.15

0.15

0.15

0.1

0.1

0.1

0.1

0.05

0.05

0.05

0.05

0.1

0.15

0.2 x (m)

0.25

0.3

0.35

0 0.34

0.4

0.3

0.25

0.25

0.2

0.2

0.15

0.1

0.05

0.05

0.05

0.1

0.15

0.2 x (m)

0

0.37

0.25

0.3

0.35

0.4

t=0.9705 s

0.15

0.1

0

0.36

0

0.05

0.1

0.15

x (m)

0.3

0

0.35

0 0.34

0.35

0.36

0.37

0.2 x (m)

0.25

0.3

0.35

0.4

0 0.34

0.3

0.3

0.25

0.25

0.2

0.2

y (m)

0.05

y (m)

0

y (m)

y (m)

y (m)

y (m)

y (m) 0.15

0

t=0.65 s

y (m)

0.3

0.15

0.1

0.05

0.05

0

0

0.05

0.1

FIGURE 12.38 Pressure contour, free surface profiles, and beam deformations at selected time instants (Sun, 2016).

0.15

0.2 x (m)

0.25

0.3

0.35

0.4

0.36 x (m)

0.37

t=1.32 s

0.15

0.1

x (m)

0.35

0 0.34

0.35

0.36 x (m)

0.37

562

Fluid
Solid Interaction Dynamics

X coordinate

(A) 0.38 0.37 0.36 0.35 0.34

0

1

2

3

4

5 Time (s)

6

7

8

9

10

(B) 5000 Elastic beam Rigid beam

p (Pa)

4000 3000 2000 1000 0

0

1

2

3

4

5 Time (s)

6

7

8

9

FIGURE 12.39 Time histories of: (A) the beam deflection at its top and (B) the water pressure at beam point Y 5 0.02 m (Sun, 2016).

FIGURE 12.40 Hybrid grid configuration in fluid domain (Javed, 2015).

FIGURE 12.41 Nodes in a support domain of point x1 in the meshfree zone (Javed, 2015).

10

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

563

FIGURE 12.42 A 2D cylinder vibration induced by an incoming flow (Javed, 2015).2D, Two-dimensional.

morphing along with the solid boundary without necessitating remeshing. The Cartesian grid, on the other hand, provides improved performance by allowing the use of the computationally efficient mesh-based method. Flow equations, in ALE formulation, are solved by local radial basis function (RBF) in FD mode on moving meshfree nodes. Conventional FD is used over static Cartesian grid for flow equations in the Eulerian formulation. The equations for solid motion are solved using the classical Runge
Kutta method. Closed coupling is introduced between fluid and structural solvers by using a subiterative prediction
correction algorithm. In order to reduce computational overhead due to subiterations, only near-field flow in the meshfree zone is solved during inner iterations. The full fluid domain is solved during outer iterations only when the convergence at the solid
fluid interface has already been reached. In the meshfree zone, adaptive sizing of the influence domain is introduced to maintain a suitable number of neighboring particles. A hybrid grid has been found to be useful in improving computational performance by faster computing over a Cartesian grid, as well as by reducing the number of computations in the fluid domain with FSI. The solution scheme was tested for problems relating to flow-induced cylindrical and aerofoil vibrations. The results are found to be in good agreement with previous work and experimental results. Here, as an example, we discuss the problem on flow-induced vibration of a 2D cylinder as follows. Motion equation of two-dimensional cylinder Fig. 12.42 shows the sketch of a 2D cylinder in the flow, of which the motion equations can be derived using Newton’s second law in the following forms: mx€ 1 cx x_ 1 kx x 5 DðtÞ; my€ 1 cyðy_ 1 ky y 5 LðtÞ; DðtÞ 5

Σ

σ1j ηj dS;

ð LðtÞ 5

Σ

σ2j ηj dS;

(12.168)

where m, c, and k denote the mass, viscous damping coefficient, and the stiffness with subindex x and y to distinguish two directions of the cylinder. D(t) and L(t) represent the resultant fluid drag and lift forces on the cylinder, obtained by the integration with respect to the Cauchy stress σij on the cylinder surface Σ of outside unit normal vector ηj . Fluid domain In the active meshfree zone, the fluid velocity is represented by an interpolation based on nodes shown in Fig. 12.41 and using a set of interpolation functions called RBFs (Kansa, 1990). The spatial derivatives of the velocity can be directly obtained by conducting the derivative operation of the RBFs, from which the original PDE of flows can be transformed into an ODE of time. Actually, this method belongs in the category discussed in Section 12.3. In the nonactive meshfree zone, the meshfree notes are overlapped by a Cartesian grid in order to match the stationary Cartesian grid in the far-field, so the ALE method is adopted to account for nodal movement.

564

Fluid
Solid Interaction Dynamics

For numerical simulations, the following nondimensional parameters are defined: Reynolds number: Reduced velocity: Lift coefficient: Drag coefficient:

ρUD ; μ U vr 5 ; fN D L CL 5 ; ρU 2 D D CD 5 ; ρU 2 D

Re 5

Mass ratio:

m 5

Frequency ratio:

fNx ; fNy

(12.169)

m ; displace fluid mass

where ρ, U, and μ denote the mass density, free stream velocity, and dynamic viscosity of the fluid, respectively; D, fNx, and fNy represent the diameter and the horizontal and vertical natural frequencies of the cylinder. In the numerical simulations, the Reynolds number 150, the frequency ratio 2, the mass ratio 2, and reduced velocity range 1
12 are chosen, but the damping coefficients for cylinder are set to zero. As shown in Fig. 12.43, the chosen parameters with the full fluid domain size are the same as the reported investigation and experiment by Dahl et al. (2010) and Zhou and Tu (2012), from which the dual resonance was found. The size of the active meshfree zone around the cylinder is set at 3D 3 3D, outside of which an overlapped meshfree zone extends by a length of 1.5D in all four directions. Solutions are obtained for both ordered and randomized meshfree nodal arrangements. Resultant amplitudes of crossflow and inflow vibrations, RMS values of the lift coefficient, and mean values of the drag coefficient are shown in Fig. 12.44 along with numerical solutions obtained by Zhou and Tu (2012) and experimental results by Dahl et al. (2010). It can be observed that the results do not change significantly with randomization of meshfree nodes. Vibration amplitudes and lift and drag coefficients tend to increase dramatically as the resonance conditions are approached near to vr 5 6, at which the crossflow amplitude achieves its maximum value 0.908. However, inflow and crossflow vibration amplitudes are very low outside the resonance range from vr 5 6. These observations are in agreement with the results by Dahl et al. (2010) and Zhou and Tu (2012). Fig. 12.45 shows the vortex structures produced by the vibrations of a 2-DOF cylinder with a frequency ratio fNx =fNy 5 2:0 at different reduced velocities.

FIGURE 12.43 The geometrical configuration of the FSI system (top), and the active and inactive meshfree grid around 2D cylinder and Cartesian grid far from the cylinder (bottom) in the numerical simulation (Javed, 2015). FSI, Fluid
solid interaction; 2D, two-dimensional.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

0.35

1.4

0.3

1

0.25

0.8

0.2

0.6

0.15

0.4

0.1

0.2

0.05 2

8 4 6 →Reduced velocity (vr)

10

0

12

4.5

1.2

4

1

3.5

0.8 0.6

0.2

1.5 1 4 6 8 →Reduced velocity (vr)

8

10

12

10

12

2.5 2

2

6

3

0.4

0

4

→Reduced velocity (vr)

1.4

0

2

0

mean CD

rms(CL)

Xmax/D

Ymax/D

1.2

0 0

565

10

12

0

2

8 4 6 →Reduced velocity (vr)

FIGURE 12.44 Nondimensional vibration amplitudes of cylinder in crossflow Ymax/D and inflow Xmax/D, RMS lift coefficient CL, and mean drag coefficient CD obtained in the simulation with parameters (m 5 2 5 fNx =fNy , Re 5 150) and meshfree nodes: ordered 2 O 2 and random 2Δ 2 , which is compared with the results
*
by Zhou and Tu (2012) based on the same parameters, and the experimental result (----) based on parameters (m 5 5:7; fNx =fNy 5 1:9, Re 5 15,000
60,000) by Dahl et al. (2010) and (Javed, 2015). RMS, Root mean square.

FIGURE 12.45 Vortex structure produced by 2-DOF cylinder vibrations in the frequency ratio fNx =fNy 5 2:0 and at different reduced velocities (Javed, 2015). DOF, Degrees of freedom.

566

Fluid
Solid Interaction Dynamics

12.4.3.5 Wave/wind energy harvesting systems Fundamental principle for wave/wind energy harvesting device design With increasing requirements for green energy to reduce environment pollution, scientists and engineers have put more effort into extracting energy from sea waves and wind. Investigations on wave energy harvesting devices have attracted wide interest around the world (for example, see Thorpe, 1999; Falnes, 2002; Office of Naval Research, 2003; Bedard et al., 2005; Rhinefrank, 2005; Wave Dragon, 2005; U.S. Department of the Interior, 2006; JAMSTEC, 2006; Ocean Power Technologies, 2006; Ocean Power Delivery Ltd., 2006; Wave Plane Production, 2006). As discussed in a review paper by Xing (2016), among these designs, the fundamental principle is to use waves/winds to excite mechanical motions of energy harvesting devices and then to convert mechanical energy to storable energies. Therefore the motions of energy harvesting devices excited by fluids are required to be as large as possible. To reach this aim, the following physical mechanisms may be adopted. Interdisciplinary research on fluid
solid interaction and electric
mechanical interaction As it is well-known that sea wave/wind energy harvesting devices are operated on or in fluids, the total system is typically an FSI. The dynamic behaviors of devices designed in dry cases are affected by fluids. A minor change of the characteristic frequency of a device can cause a large difference of dynamic response, as demonstrated by Xing et al. (2009b). Therefore FSIs have to be considered to design effective energy harvesting devices. Further as studied in the papers by Xing et al. (2009b, 2011), electric units in energy harvesting devices cause mechanical behaviors affected by electromagnetic dynamics involving electric
mechanical interactions (EMIs). As a result, for efficient energy harvesting designs, interdisciplinary research concerning linear/nonlinear dynamics, fluids, solids, and electric systems, as well as their interactions, is necessary. Resonance Resonance is a vibration to be avoided in engineering designs of structures, since a very small excitation force from the environment causes a very large motion of the system, so that the designed structure can fail in its dynamic strength or cannot work normally in the unwanted dynamic environment. However, in a reverse consideration, the resonance phenomenon can be adopted in order to design a linear device with its natural frequency close to the wave frequency, whose resonance is expected to create large mechanical motions by excited waves for energy harvesting devices. In considering FSI and using the developed numerical method (Xing et al., 1991, 1996) with computer code fluid
structure interaction analysis program (Xing, 1992a,b, 1995a,b), a wave energy harvesting device
water interaction system was subjected to the wave maker excitation in a towing tank (Xing et al., 2009b). The results demonstrated that FSI changes the natural frequency of the device designed in dry cases and therefore obviously affects the system efficiency. More importantly, this investigation revealed that the collected energy plays a role like an active damping added to the system, so that it is useful to keep a stable oscillation at its resonance state of the system. In Section 7.4.9.7, we have discussed this linear oscillator
water interaction system based on the resonance of FSI system to harvest large amounts of wave energy. Periodical orbit of nonlinear system Another idea is to design a nonlinear energy harvesting system and to use its inherent large stable orbit motion to extract energy. A number of researchers (Xu et al., 2005, 2007; Litaka et al., 2008, 2010; Lencia et al., 2008; Horton et al., 2011; Nandakumar et al., 2012; Pavlovskaia et al., 2012) proposed and investigated the possible rotational motions of a nonlinear pendulum subject to different base motion excitations aiming to harvest wave energy, but unfortunately they did not consider FSI effect. To amend this, Xing et al. (2011) proposed a mathematical model to study the dynamics of nonlinear oscillators, pendulums, and smooth and discontinuous (SD) oscillators (Cao et al., 2006, 2008a,b) coupling with water, as well as an electric-magnetic energy converter aiming to harvest wave energy. In the model, the pressure in the finite water domain is chosen to describe the water motions, and the paper qualitatively discusses the integrated FSI/EMI characteristics and proposed numerical solution approaches. Moreover, using the potential of velocity as a variable and based on the Green function satisfying the free surface wave condition and the Green identity, in Section 9.3.4.2 we discussed a nonlinear pendulum
water interaction system (Xing et al., 2011) to use its periodical solution for harvesting energy. By this investigation, we have further demonstrated that the energy converter in integrated systems acts as a damping mechanism to keep the excited large periodical motions in a stable state of motion. The nonlinear energy flow theory (Xing, 2015a,b) can be used to analyze more complex system based on energy conservation laws. Aerofoil flutter system Flutter is a very harmful vibration for airplanes that has to be avoided (Bisplinghoff et al., 1955; Bisplinghoff, 1958; Bisplinghoff and Ashley, 1962; Fung, 1955, 1969). However, harmful flutter vibration could be used to harvest wave/wind energy. The idea is to design an aerofoil with the energy converter excited by flows. Since the energy converter has a damping effect (Xing et al., 2009b; 2011; Xing 2016), the flutter generated is a stable oscillation to be practically used to harvest wave/wind energy.

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

567

FIGURE 12.46 Schematic model of oscillating aerofoil power extraction device (Yang, 2013).

Here, we present the two types of flutter aerofoil energy harvesting systems: the active and the semipassive designs. As discussed in the review papers (Young et al., 2014a,b; Xiao and Zhu, 2014), the full- and semipassive designs require an excitation by an input energy from another source to keep the device in the expected motion, while the full active one does not need this extraexcitation energy so that is more conducive to harvesting natural energies. Active aerofoil The energy harvesting device studied by Yang et al. (2011) is shown in Fig. 12.46. A 2-DOF aerofoil is placed in a uniform airflow of velocity U, with h denoting heaving displacement and α being positive nose up, representing pitching angle. The moving coil is rigidly connected to aerofoil at its elastic axis by a massless bar. The static equilibrium position where h 5 0 5 α is taken as a reference position. The aerofoil is considered a unit span in the direction perpendicular to the paper plane. Based on quasisteady aerodynamic theory, its governing equations of motion are written as follows (Fung, 1969; Dovel and Ligamov, 1988; Lee and Jiang, 1999). Equations of aerofoil motion mh€ 1 Sα€ 1 ch h_ 1 kh h 5 F 2 FM ;

Sh€ 1 Iα α€ 1 cα α_ 1 kα α 5 M;

kh 5 kh1 1 kh2 h2 ;

kα 5 kα1 1 kα2 α2 ; (12.170)

where m, S, and I are the total mass of the aerofoil and electric coil, the static and second moments of the aerofoil mass, respectively; c and k, respectively, represent the damping and stiffness coefficients identified by subindex h for vertical motion and α for pitching rotation; F and M denote the resultant lift and its moment acting on the aerofoil by the flow, and FM is the magnetic force. The stiffness coefficients of springs are assumed nonlinear, denoted by Eq. (12.170). Eq. (12.170) are written in a nondimensional form  2  ω ω  1 Fm b Cl ðτ Þ 1 ξv 1 xα αv 1 2ζ ξ  ξ 0 1 ξ 1 β ξ ξ3 5 2 ; (12.171a)  U U πμ mU 2  xa ζ 1  2 ξv 1 αv 1 2 α α0 1  2 α 1 β α α3 5 Cm ðτ Þ; U πμra2 ra2 U where h U Ut ωξ ; ω5 ξ 5 ; U 5 ; τ5 ; b bωα b ωα sffiffiffiffiffiffiffi kα1 ch cα ωα 5 ; ζξ 5 ; ζα 5 ; Iα 2mωξ 2Iα ωα

sffiffiffiffiffiffi kh1 ; ωξ 5 m μ5

m ; πρb2

xα 5

S ; mb

rα 5

sffiffiffiffiffiffiffiffi Iα ; mb2

(12.171b)

βξ 5

kh2 ; kh1

βα 5

kα2 : kα1 (12.171c)

568

Fluid
Solid Interaction Dynamics

Here, a prime ðÞ0 denotes derivative of ðÞ with respect to nondimensional time τ, and μ, xα, and rα are the mass ratio, the distance of mass center after the elastic axis, and the radius of gyration about the elastic axis, respectively; the values of lift force and moment were obtained by adopting a quasisteady aerodynamic theory (Fung, 1969), that is, F 5 2 ρbU 2 Cl ðτ Þ;

(12.172a)

M 5 2ρb U Cm ðτ Þ;

(12.172b)

2

2

where Cl ðτ Þ and Cm ðτ Þ are the nondimensional lift and pitching moment coefficients, respectively. Given incompressible flow, they are given by the expressions as in Fung (1969) and Lee and Jiang (1999), that is,       ðτ 1 1 ðσÞ 0 0 ðσÞ 0 2 ah α ð0Þ φðτ Þ 1 2π φðτ 2 σÞ α ðσÞ 1 ξvðσÞ 1 2 ah α αv dσπ Cl ðτ Þ 5 2π αð0Þ 1 ξ ð0Þ 1 2 2 (12.172c) 0 1 ðξv 2 ah αv 1 α0 Þ;

      π 1 π 0 π 1 1 ah ðξv 2 ah αvÞ 2 2 ah α 2 αv 1 π 1 ah αð0Þ 1 ξ0 ð0Þ 1 2 ah α0ð0Þ φðτ Þ 2 2 2 16 2 2    ð τ  1 1 1 ah 2 ah αvðσÞ dσ; 1π φðτ 2 σÞ α0 ðσÞ 1 ξvðσÞ 1 2 2 0

C m ðτ Þ 5

(12.172d)

where the Wagner function φðτ Þ is expressed as (Fung, 1969) φðτ Þ 5 1 2 Ψ1 e2E1 τ 2 Ψ2 e2E2 τ :

(12.172e)

Here, the constants’ values Ψ1 5 0:165, Ψ2 5 0:0335; E1 5 0:0455, and E2 5 0:3 are given by Jones (1940). Electromagnetic power generation Electric power generation obeys the Laplace theorem describing electromagnetic phenomena (Kittel, 1967; Xing et al., 2009b): a voltage will be induced by a moving coil of effective length l in the magnetic field of intensity B, that is, _ eðtÞ 5 Blh:

(12.173a)

When the electric coil is in an electric circuit with resistance R and inductance L, a dynamic current i is introduced. Assuming the electrical conductance is zero, the equation for the electric circuit is obtained by applying a Kirchhoff’s voltage law: diðtÞ 1 RiðtÞ 2 e 5 0: dt Combining Eqs. (12.173a) and (12.173b), we obtain a differential equation with respect to τ: L

i0 1

Rb Blb 0 i2 ξ 5 0: LU L

(12.173b)

(12.173c)

Since the coil moves in the magnetic field, an electromagnetic force will be introduced FM 5 iBl:

(12.173d)

pg ðτ Þ 5 Ri2 :

(12.173e)

The power harvested by resistor R is

Solution approach The integral terms in Eqs. (12.172c) and (12.172d) make it difficult to study the system’s dynamical behaviors analytically. To assist in investigation, as used by Lee and Jiang (1999), four new variables are introduced: ðτ ðτ 2E1 ðτ2σÞ αðσÞdσ; w2 5 e2E2 ðτ2σÞ αðσÞdσ; w1 5 e 0 ðτ ð 0τ (12.174) 2E1 ðτ2σÞ w3 5 e ξðσÞdσ; w4 5 e2E2 ðτ2σÞ ξ ðσÞdσ: 0

0

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

569

Based on this definition, Eqs. (12.171a), (12.171b), and (12.173b) can be rewritten as c0 ξv 1 c1 αv 1 c2 ξ 0 1 c3 α0 1 c4 ξ 1 c5 ξ3 1c6 α 1 c7 w1 1 c8 w2 1 c9 w3 1 c10 w4 1 c11 i 5 f ðτ Þ; 0

0

d0 ξv 1 d1 αv 1 d2 α 1 d3 α 1 d4 α 1 d5 ξ 1d6 ξ 1 d7 w1 1 d8 w2 1 d9 w3 1 d10 w4 5 gðτ Þ; 3

0

0

n1 ξ 1 n2 i 1 n3 i 5 0; where 2 f ðτ Þ 5 μ

(12.175a) (12.175b) (12.175c)

     1 2 ah αð0Þ 1 ξ ð0Þ ψ1 E1 e2E1 τ 1 ψ2 E2 e2E2 τ ; 2

(12.175d)

ð1 1 2ah Þf ðτ Þ ; 2ra2

(12.175e)

gð τ Þ 5 2 with the following coefficients:

1 ah ω 2 ; c1 5 xα 2 ; c2 5 2ζ ξ  1 ð1 2 Ψ1 2 Ψ2 Þ; μ U μ μ

2 2 1 c3 5 ½1 1 ð1 2 2ah Þð1 2 Ψ1 2 Ψ2 Þ; c4 5 Uω 1 ðE1 Ψ1 1 E2 Ψ2 Þ; μ μ 8 9   < =

2 2 1 2 a h ð E 1 Ψ 1 1 E 2 Ψ2 Þ ; c5 5 Uω β ξ ; c6 5 ð 1 2 Ψ1 2 Ψ2 Þ 1 ; μ: 2 2 3 2 3     2 1 2 1 2 ah 5; c8 5 E2 Ψ2 41 2 E2 2 ah 5 ; c7 5 E1 Ψ1 41 2 E1 μ 2 μ 2 c0 5 1 1

c9 5 2

2E21 Ψ1 2E2 Ψ2 ; c10 5 2 2 ; μ μ

xa ah 8a2 1 1 ζ 1 2 2ah 2 2 ; d1 5 1 1 h 2 ; d2 5 2 α 2 2 8μra ra μra U 2μra2 8 9   = 1 1 1 2ah < 1 2 a d3 5  2 ð Ψ 1 E Ψ Þ 2 2 ð 1 2 Ψ 2 Ψ Þ ; E h 1 1 2 2 1 2 ; U 2μra2 : 2 d0 5

d4 5

βα ð1 2 Ψ1 2 Ψ2 Þð1 1 2ah Þ ;  2 ; d5 5 2 μra2 U

d6 5 2

d8 5 2 d10 5

(12.175f)

   1 2 ah E1 Ψ1 ð1 1 2ah Þ 1 2 E1 2

E 1 Ψ 1 1 E 2 Ψ2 ð1 1 2ah Þ; d7 5 2 μra2    1 2 ah E2 Ψ2 ð1 1 2ah Þ 1 2 E2 2 μra2

μra2

; d9 5

E21 Ψ1 ð1 1 2ah Þ μra2

E22 Ψ2 ð1 1 2ah Þ ; μra2

n1 5 2

Blb Rb ; n2 5 1; n3 5 : L LU

Introducing nine new parameters, x1 5 α; x2 5 α0 ; x3 5 ξ; x4 5 ξ 0 ; x5 5 w1 ; x6 5 w2 ; x7 5 w3 ; x8 5 w4 ; x9 5 i;

(12.176a)

570

Fluid
Solid Interaction Dynamics

we transform Eqs. (12.175a)
(12.175c) into a set of first-order differential equations: 0

x1 5 x2 ; c 0 H 2 d0 P 0 x2 5 ; d0 c1 2 c0 d1 0

x3 5 x4 ; 0

x4 5 2

c 1 H 1 d1 P ; d0 c 1 2 c 0 d1

0

x 5 5 x 1 2 E1 x 5 ;

(12.176b)

0

x 6 5 x 1 2 E2 x 6 ; 0

x 7 5 x 3 2 E1 x 7 ; 0

x 8 5 x 3 2 E2 x 8 ; n1 x 4 1 n3 x 9 0 x9 5 2 ; n2 where P 5 c2 x4 1 c3 x2 1 c4 x3 1 c5 x33 1 c6 x1 1c7 x5 1 c8 x6 1 c9 x7 1 c10 x8 1 c11 x9 2 f ðτ Þ;

(12.176c)

H 5 d2 x2 1 d3 x1 1 d4 x31 1 d5 x4 1 d6 x3 1d7 x5 1 d8 x6 1 d9 x7 1 d10 x8 2 gðτ Þ:

(12.176d)

By using a fourth-order Runge
Kutta method, those first-order differential equations can be numerically integrated. Solution results Fig. 12.47 shows the phase portrait obtained, which clearly indicates the limit cycle oscillation of the aerofoil. It is interesting to notice that incorporated electric generator results in decreases of amplitudes of both pitching and heaving motions, implying a damping effect on the integrated system. While the pitching amplitude reduces slightly, that of the heave decreases sharply. The reason is that the electromagnetic force created by the coil acts as resistance for the heaving motion. The differences in displacement and velocity between the two cases, with and without considering the electrical
magnetic system, can be used to estimate the energy harvested. Fig. 12.48 shows the time history of the dynamic current of the system, which oscillates at a constant amplitude of about 8.63 A. This finding demonstrates that, by adding a magnetic energy converter, the unstable flutter aerofoil system behaves as a stable oscillation system due to the energy converter’s damping, which is very important to guarantee its practical applications. The maximum instantaneous power generation, which is approximately 75.58 W, can thus be obtained by using Eq. (12.173e). Assuming that the current variation is sinusoidal, we obtain the average power of 37.79 W. Noticing that in all formulations, only a unit span value is considered, if the span of the aerofoil is 5 m, the average power would be about 188.95 W. Semiactive aerofoil For the preceding active aerofoil calculations, the quasisteady aerodynamic theory in terms of the prescribed fluid lift and moment functions in Eqs. (12.172a)
(12.172d) has not accurately given the dynamic forces in the process.

FIGURE 12.47 Phase portrait showing the limit cycle oscillations of the system: (A) pitch motion and (B) heave motion. Solid line: without electro magnetic generator; dashed line: with electromagnetic generator ðU 5 6:5; ω 5 0:2; ζ ξ 5 0:01; ζ α 5 0:01; μ 5 100, xα 5 0:25, rα 5 0:5, ah 5 2 0:5, β ξ 5 1, β α 5 1:5; B 5 0:5 T, l 5 20 m, R 5 1 Ω, L 5 0:05 H, b 5 1 m; ρ 5 1:293 kg=m3 Þ (Yang, 2013).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

571

FIGURE 12.48 Dynamic current at steady state considering the electromagnetic generator. Parameter values given in the caption for Fig. 12.46 (Yang, 2013).

FIGURE 12.49 Aerofoil flutter system studied by Javed (2015) and Javed et al. (2016).

Javed (2015) and Javed et al. (2016) adopted a coupled meshfree
mesh-based fluid solver described in Section 12.4.3.4 to simulate a semiactive aerofoil system, as shown in Fig. 12.49. In this case, the aerofoil is subjected to a prescribed pitching motion about its elastic axis and is allowed to move freely along the heave axis due to fluid forces. The aerofoil is mounted on a translational spring-damper system. When the aerofoil is subjected to periodic pitch oscillation, it causes a corresponding variation in the fluid forces over time. These time-varying fluid forces induce heaving motion. Such mechanisms have recently gained attention for their potential application in tidal and wind energy extraction systems (Wu et al., 2014; Deng et al., 2015), in which the pitching motion is prescribed by other excitation and both the vertical and the pitching motions are used to harvest the flow energy. Equations of aerofoil motions Aerofoil NACA0015 is chosen of chord length c with its electric axis located at a distance c/3 from the leading edge to complete the numerical simulation. The vertical h(t) and pitching θðtÞ motions of the aerofoil in the fluid are governed by the nonlinear PDEs similar to the ones given by Dubcova et al. (2009), that is, 2 ^ mh€ 1 Sθ€ cos θ 2 Sθ_ sin θ 1 ch h_ 1 kh h 5 LðtÞ 1 L; ^ Sh€ cos θ 1 Iθ θ€ 1 cθ θ_ 1 kθ θ 5 MðtÞ 1 M;

(12.177a)

which, when the pitching angle is small, reduces to ^ mh€ 1 Sθ€ 1 ch h_ 1 kh h 5 LðtÞ 1 L; € _ € ^ Sh 1 Iθ θ 1 cθ θ 1 kθ θ 5 MðtÞ 1 M;

(12.177b)

572

Fluid
Solid Interaction Dynamics

which is similar to Eq. (12.170), except that no magnetic force is included. Here, L^ and M^ denote external excitations to produce the prescribed pitch motion, and other notations are the same as the ones in Eq. (12.170), except α is replaced by the prescribed rotation angle θ. The produced lift force L(t) is calculated by Eq. (12.168), and the moment M(t) of the lift force is given by ð MðtÞ 5 e3ki r k σij ηj dS; (12.177c) Σ

where rk denotes the position vector from the elastic axis to a lift point on the surface of aerofoil. The prescribed pitching motion is defined by θðtÞ 5 θ0 cos ωt;

(12.178)

where θ0 and ω denote the amplitude and frequency of this harmonic motion, respectively. Fluid motion The fluid of mass density ρ and the free stream velocity U are assumed to be viscous with their chosen Reynolds number Re 5 ρUc=μ 5 1100; 1000, of which the flow is predominantly laminar (Kinsey and Dumas, 2008). The grid arrangement of the fluid domain is given by Fig. 12.40, for which Fig. 12.50 gives more details. As sketched in Fig. 12.50B, the aerofoil is placed at a distance of 4c from the inlet and 12c from outlet, and the width of fluid domain is set to 10c. The dimensions of active meshfree zone around the aerofoil are set to 1.35c 3 1.6c. The velocity of the flow particle on the aerofoil surface is required to equal the velocity of aerofoil motion, and the pressure value at the boundary is obtained by solving the pressure Poisson equation based on the Neumann boundary condition, as discussed in Eqs. (12.162a) and (12.162b). A total of 300 nodes are placed on the aerofoil surface, and the computational domain comprises a total of 25,880 meshfree and 65,354 Cartesian nodes. The movement of mesh is accomplished by displacing the meshfree zone according to the prescribed pitching motion. The time step is chosen as 0.001 second. Solution The solution process is the same as used for the flow-induced vibrations given in Section 12.4.3.4, of which the details can be found in the thesis by Javed (2015) and the paper by Javed et al. (2016). The aerofoil equation is solved in a nondimensional form by using the following parameters: m 5

m ; 0:5ρc2

c 5

c ; 0:5ρUc

k 5

k ; 0:5ρU 2

f 5

ω : 2π

(12.179)

FIGURE 12.50 Hybrid grid around aerofoil: (A) arrangement of Cartesian grid and meshfree particles, (B) full domain, (C) leading edge, (D) trailing edge (Javed, 2015).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

573

TABLE 12.2 Comparison of the simulation result by Javed et al. (2016) (J) with that by Wu et al. (2014) (W) and that by Deng et al. (2015) (D). Mechanical parameters θ0 15

(degree)

m*

c*

k*

f*

1

10

10

0.2

π

30

75

75

0.1022

0

Re

Source

CLmax

CD

1100

W

0.704

0.179

J

0.690

0.170

0.1

0.12

0.22

1000

max CM

W

0.905

0.345

J

0.885

0.334

D

2.0

0.33

J

2.017

0.31

D

2.8

0.6

J

2.550

0.56

The simulations are run at two chosen Reynolds Re 5 1100, 1000 with four different mechanical parameters. The simulation results are compared with the ones reported by Wu et al. (2014) and Deng et al. (2015) in Table 12.2, from which a good agreement is found. The time histories of heave/pitch displacements and lift/drag coefficients of the aerofoil subjected to the prescribed pitch harmonic motion in Eq. (12.178) with its amplitude θ0 5 76:33 degrees and frequency ω 5 0:28π are shown in Fig. 12.51. The obtained vorticity profiles at selected time instants are given in Fig. 12.52. From these two figures, it is observed that with the pitch angle increasing in the initial phase t , T/8 of oscillation, the lift achieves its maximum value, and the flow remains largely attached with the aerofoil’s top surface as shown in Fig. 12.52A. The first peak in lift profile appears at around t 5 T/8. Increasing lift also causes an increase in drag, and therefore the drag coefficient also increases. This initial rise in lift is followed by flow separation close to the leading edge, as shown in Fig. 12.52B, causing a reduction in lift. Subsequently, the detached leading edge vortex reattaches with the aerofoil close to its trailing edge (Fig. 12.52C) at about t 5 3T/8, causing a second peak in the lift profile. However, as the leading edge vortex leaves the aerofoil from the trailing edge and moves farther downstream, a sharp decline in lift is observed between 3T/ 8 , t , 5T/8. The lift coefficient reduces to zero and then shows a similar profile in the negative direction. These results are in good agreement with the reference values by Kinsey and Dumas (2008). Fig. 12.53 shows the two time histories of lift coefficient affected by the excitation frequency. It can be observed that the peak value of the lift coefficient increases at a higher frequency f  . Similar behavior was observed by Wu et al. (2014). Harvested energy As mentioned, this semiactive energy harvesting system requires the external force M^ to keep the aerofoil in the given pitch motion, so that external input power is needed, which reduces the net harvesting power of the system. In this simulation, the main aim is to investigate the dynamic behavior of the aerofoil energy harvesting system based on the NS PDEs solved by th meshfree particle approach, and the energy converter is not included. We have learned that the energy converter in the system acts as a damping mechanism; therefore we can consider the damping coefficients ch and cθ in Eqs. (12.177a)
(12.177c) to reflect the energy collection damping for the heave and pitch motions, respectively, of the aerofoil. From this assumption, we can estimate the harvested power of the system as follows. Instantaneous energy flow equilibrium equations The dynamic equations in Eqs. (12.177a)
(12.177c) are now rewritten in matrix form:         h LðtÞ m S cos θ h€ kh 0 L^ ch 2Sθ_ sin θ h_ 5 1 ^ : (12.180) 1 1 _ € MðtÞ 0 kθ θ S cos θ Iθ θ M 0 cθ θ

Premultiplying Eq. (12.180) by the row matrix h_ θ_ , we obtain its energy flow equilibrium equation (Xing and Price, 1999; Xing, 2015) in the form

1

100

h(t) θ(t)

50

0

0

–0.5

–1

θ(t)

h(t)

0.5

–50

0

T/4

T/2 Time (t)

3T/4

T

2

–100

4

CL

1

3

0

2

–1

1

–2

0

T/4

T/2 Time (t)

3T/4

T

CD

CL

CD

0

FIGURE 12.51 Time histories of the heave displacement h(t), pitch angle θðtÞ, lift coefficient CL, and drag coefficient CD of the NACA0015 aerofoil subjected to the prescribed pitch motion in Eq. (12.178) at Re 5 1100, θ0 5 76:33 degrees, and ω 5 0:28π (Javed, 2015).

FIGURE 12.52 Instantaneous vorticity profiles around NACA0015 at Re 5 1100, Re 5 1100, θ0 5 76:33 degrees, and ω 5 0:28π (Javed, 2015).

FIGURE 12.53 The time history of lift coefficients of NACA0015 aerofoil subject to the pitch oscillation with its amplitude θ0 5 75 degrees and two frequencies, f  5 0:12; 0:22 (Re 5 1000, m 5 0:1022, c 5 π, k 5 0) (Javed, 2015).

Mixed finite element—smoothed particle methods for nonlinear fluid
solid interactions Chapter | 12

  



^ dT dΠ 1 m S cos θ h_ ^I 5 h_ θ_ L ; 1 PH 1 5 P^I ðtÞ 1 PF ðtÞ; T 5 ; P h_ θ_ S cos θ Iθ θ_ M^ dt dt 2     

kh 0

ch 0

LðtÞ 1 h h_ h θ ; PF 5 h_ θ_ ; PH 5 h_ θ_ ; Π5 0 cθ θ_ 0 kθ θ MðtÞ 2 where we have used the following equation in deriving this result: " #" #  " #



m S cos θ h€ dT 1 0 2Sθ_ sin θ h_ _ _ _ _ 5 h θ 1 h θ dt 2 S cos θ Iθ θ_ 2Sθ_ sin θ 0 θ€ ( " # " # " # )

m S cos θ h€ 0 2Sθ_ sin θ h_ _ _ 1 : 5 h θ S cos θ Iθ θ_ θ€ 0 0

575

(12.181a)

(12.181b)

Here, T and Π, respectively, denote the kinetic energy and potential energy of the system, while PH is the harvested power, and P^I and PF represent the input power of the external excitation force and the power done by the aerodynamic forces, respectively. Time-averaged energy flow equilibrium equation Taking a time average of Eqs. (12.181a) and (12.181b) by the operation for a time function P, ð hPi 5 1T^ Pdt 0T^

over a vibration period T^ 5 2π=ω, we obtain the time-averaged energy flow equilibrium equation as         dT dΠ 1 hPH i 1 5 P^I ðtÞ 1 PF ðtÞ : dt dt Net harvested averaged power

  hPN i 5 hPH i 2 P^I

(12.182a)

(12.182b)

(12.182c)

Here it is necessary to remember the following important points for designing nonlinear energy harvest systems: 1. The dynamic equation shown in Eq. (12.180) is a nonlinear equation, in which the related matrices, especially the fluid dynamic forces, depend on the motion of the system, so that the excited motion and the fluid forces of the nonlinear system produced by the harmonic pitch oscillation are not necessarily harmonic variables with the same frequency as the pitch oscillation. The energy flow variables in Eqs. (12.181a), (12.181b), and (12.182a)
(12.182c) generally must be calculated by numerical integrations. 2. For a nonlinear system, the time-averaged kinetic/potential energies do not necessarily vanish (Xing, 2015), imply^ and ΠðtÞ 6¼ Πðt 1 TÞ; ^ therefore the time change rates of the kinetic and ing that, in general cases, TðtÞ 6¼ Tðt 1 TÞ potential energies in Eq. (12.182b) remain. Physically, this represents that some energy exchange occurs between the solid system and the flow and/or excitation forces during a given time period. 3. For the semiactive aerofoil system, the arrangement of the vertical/rotational stiffens, providing the possibility of adjusting the natural frequencies of the solid system in order to obtain the required large motions and the harvested energy. To confirm this, analysis of the net harvested averaged power affected by the design parameters is needed.