Hydrodynamic Loads Computation Using the Smoothed Particle Methods

The study of wave fluid flows is now under special consideration in view of serious effects, caused by dams breaking and consequent formation of moving waves, their interaction with solids and structures, uprush on shore, etc. Thereby solving the problem of hydrodynamic loads estimation is important for designing the shape and stiffness of the structures, interacting with oncoming waves. Such problems, due to large deformations of free surfaces, are very complex, and meshless methods proved to be the most suitable for numerical simulation of them. Particle methods form the special class of meshless methods, which mainly based on the strong form of governing equations of gas dynamics and fluid dynamics. The peculiar representatives of particle methods are Smoothed Particle Hydrodynamics (SPH) (Lucy, 1977; Gingold & Monaghan, 1977) and Incompressible SPH (ISPH) (Cummins & Rudman, 1999; Shao & Lo, 2003; Lee et al., 2008). Large amount of papers, devoted to numerical simulations of free surface flows using SPH or ISPH, demonstrated a high degree of efficiency of both methods in obtaining the kinematic characteristics of flows, though it has been revealed, that ISPH shows a larger particle scattering at the stages, following the water impact, in comparison with the classic SPH, where particles are more ordered. However, dynamic characteristics of flows are still hard to compute, especially it concerns the classic SPH. The objective of the chapter is to analyze the capacity of the methods to compute pressure fields and hydrodynamic loads subsequently.


Introduction
The study of wave fluid flows is now under special consideration in view of serious effects, caused by dams breaking and consequent formation of moving waves, their interaction with solids and structures, uprush on shore, etc. Thereby solving the problem of hydrodynamic loads estimation is important for designing the shape and stiffness of the structures, interacting with oncoming waves.Such problems, due to large deformations of free surfaces, are very complex, and meshless methods proved to be the most suitable for numerical simulation of them.Particle methods form the special class of meshless methods, which mainly based on the strong form of governing equations of gas dynamics and fluid dynamics.The peculiar representatives of particle methods are Smoothed Particle Hydrodynamics (SPH) (Lucy, 1977;Gingold & Monaghan, 1977) and Incompressible SPH (ISPH) (Cummins & Rudman, 1999;Shao & Lo, 2003;Lee et al., 2008).Large amount of papers, devoted to numerical simulations of free surface flows using SPH or ISPH, demonstrated a high degree of efficiency of both methods in obtaining the kinematic characteristics of flows, though it has been revealed, that ISPH shows a larger particle scattering at the stages, following the water impact, in comparison with the classic SPH, where particles are more ordered.However, dynamic characteristics of flows are still hard to compute, especially it concerns the classic SPH.The objective of the chapter is to analyze the capacity of the methods to compute pressure fields and hydrodynamic loads subsequently.

Governing equations
The governing equations of fluid dynamics, including the Navier-Stokes equations and the continuity equation, in the case of the Newtonian viscous compressible fluids, are of the following form: (2) www.intechopen.com Hydrodynamics -Optimizing Methods and Tools

52
where ab 123  -numerical indices of coordinates, a v -components of the velocity vector, a F -components of the vector of volumetric forces density, ab  -Kronecker symbols, p and  -pressure and density of the fluid, correspondingly.Here the Einstein summation convention is assumed.The viscous stress tensor components are calculated by the formula (  -dynamic viscosity): (3) For enclosing the system (1)-( 3) one should make some assumptions about fluid properties.The original SPH method assumes the fluid to be weakly compressible, and therefore is applied to the system (1)-( 3) with certain equation of state for enclosure.The most often used equation of state is the Theta form equation for barotropic processes (Monaghan et al., 1994): Selecting the coefficient of volume expansion B one can obtain the effect of incompressible fluid.
The ISPH method in contrast to the original SPH uses the model of incompressible fluid, what means dd t /0   .In that case the equation of state shouldn't be considered and the enclosed system of governing equations takes the following form:

The basis of the methods
The key idea of smoothed particle methods lies in discretization of the problem domain into a set of Lagrangian particles, which play the role of nodes in function approximation.For construction of approximation formulas in smoothed particle methods the exact integral representation with the Dirac  -function is used: The Dirac  -function is changed here by a compactly supported function W , called the kernel function, what allows to obtain the integral formula about the bounded domain: The value h determines a size of support domain D of the function W and is called a smoothing length.Having a set of particles scattered about the problem domain  we can estimate the value of the above integral with the quadrature (Lucy, 1977;Gingold & Monaghan, 1977): where n is a number of particles, determined as "nearest neighbours" of the i -th particle within the support domain D .Two particles i and j are called neighbouring or interacting particles, if the distance between their centers does not exceed kh , where k depends on the type of kernel function and -radius-vector, mass and density of the j -th particle, correspondingly.A simple formula for the gradient of a function has the form:

Kernel function
As kernel function is a keystone of smoothed particle methods a great attention is paid to construction of new types of kernels.Till now a large amount of different types of kernel functions have been developed.All of them should satisfy the following basic conditions: Here for the problems, simulated with SPH, the original Monaghan's cubic spline is utilized (Monaghan et al., 1994): where q h    rr .
As it was pointed out (G.R. Liu & M.B. Liu, 2008) the approximations of functions based on the kernels that haven't smooth second derivative are too sensitive to particle scattering.It plays a crucial role for the ISPH method as elliptic Poisson equation is solved for obtaining a pressure field.That is why in numerical simulations using ISPH the fourth-order spline has been used (Morris, 1996;Lee et al., 2008): www.intechopen.com Hydrodynamics -Optimizing Methods and Tools 54

Approximation of governing equations
For approximation of gradient terms in equations (1) or (5) the original formula (10) may be applied.However, it is usually implemented for derivation of new forms of gradient approximations.In numerical simulations the following form is commonly used: This formula has an advantage of being symmetric in relation to interacting particles and thus conserves total momentum of a system of particles, representing the problem domain.
Besides it gives more stable results of numerical simulations in comparison to (10).
For a divergence of a velocity field in the continuity equation ( 2) the following expression is usually applied: The above form gives a zero-valued first derivatives for a constant field.Using ( 13) for approximation of gradient of a function one can obtain the following discrete representation for viscous term in equation ( 1): Normal and tangent components of viscous stress tensor T i are defined by following expressions similar to (14) (G.R. Liu & M.B. Liu, 2008): As it will be pointed out in section 3.4 the pressure Poisson equation need to be solved in the ISPH method.There are some ways to obtain the approximations of second derivatives in smoothed particle methods.One way consists in directly deriving the formula in a similar manner as for first derivative (10).The idea of the other is in subsequent implementation of a gradient formula ( 13) and a divergence of vector field ( 14).However these ways proved to be too sensitive to inhomogeneous particle distribution and result in non-physical oscillations of pressure field.So the approximation of the first derivative in terms of the SPH method and its finite difference analogue are usually applied together according to Brookshaw's idea (Brookshaw, 1985).Based on it some different forms of Laplacian operator were derived (Cummins & Rudman, 1999;Shao & Lo, 2003;Lee et al., 2008).Here for numerical simulations the form of Lee (Lee et al., 2008) is used: www.intechopen.com The approximation formulas for viscous forces in ISPH are obtained in a similar way and may take different forms (Cleary & Monaghan, 1999;Shao & Lo, 2003).Here for numerical simulations the following viscous term by Morris (Morris et al., 1997) is utilized:

Time integration
In the original SPH method for time integration the "predictor-corrector" scheme is commonly used: "predictor": "corrector": The new radius-vectors of particles on (1 ) n  -th time step are calculated using the Euler integration scheme: For time integration of motion equations in the ISPH method the split step scheme is applied (Yanenko, 1960;Chorin, 1968).According to its idea time integration process is splitted into convection-diffusion and pressure contribution.So the first step of the scheme for preliminary velocity values takes the from: Projecting the preliminary velocity values onto a null-divergence field one can obtain: provided the pressure field on n (1 )  -th time step is calculated through the pressure Poisson equation (Lee et al., 2008): www.intechopen.com where the velocity divergence at right hand side of above equation is calculated using formula ( 14).The radius-vectors of particles on n (1 )  -th time step can be get out of the following formula according to Euler explicit integration scheme: The equation ( 24) is reduced to the system of linear algebraic equations with symmetric matrix.For solving this system the preconditioned generalized minimum residual method (PGMRES) is applied (Saad, 2003).

Free surface
For identification of particles on the free surface, one can apply some different ways.One of such ways is using the geometrical Dilts algorithm (Dilts, 2000), based on the fact, that each particle has its size, commonly determined by the smoothing length.
The other way is detection of particles, satisfying the inequality (Lee et al., 2008): as free surface particles have less nearest neighbors in comparison with the inner ones.
Here the Dilts algorithm is used for the original SPH method and the formula (26) for ISPH.
For free surface particles the Dirichlet condition is imposed: p 0  .For the original SPH it means that free surface particles has the zero pressure, not the pressure obtained out of the equation of state as for inner fluid particles.As the pressure Poisson equation (PPE) is solved in the ISPH method for obtaining pressure field, the Dirichlet condition is embedded into the matrix of the system of linear algebraic equations (SLAE), which is the discrete representation of PPE.This procedure conserves the symmetry of matrix of SLAE.

Solid boundary
In smoothed particle methods the most commonly way of imposing conditions at solid boundaries is the virtual particle method, divided into two basic types.The first type -Monaghan virtual particles method (Monaghan et al., 1994).The virtual particles are located along the solid boundary in a single line, don't change their characteristics in time, and effect on the fluid particles by means of a repulsive force, based on certain interaction potential.The most popular among researchers is the Lennard-Jones potential.The second type -Morris virtual particles (Morris et al., 1997).These particles are located along the solid boundary in several lines.The number of the lines depends on the smoothing length of particles of the fluid.This allows solving one of the main problems of the SPH method -asymmetry of the kernel function near the boundaries.The effect of the Morris particles on the fluid ones differs from the effect of Monaghan particles by the fact, that there is no need in using any interaction potential.Instead of this, values of the characteristics in the Morris particles are calculated on the basis of their values in particles of the fluid.Here for imposing solid boundary conditions on velocity the Morris virtual particles are used for both methods.In ISPH the Morris virtual particles are also implemented for imposition of Neumann boundary conditions on solid walls, that is /0 pn   (Koshizuka et al., 1998;Lee et al., 2008).The procedure of embedding these conditions into the matrix of SLAE breaks its symmetry.Therefore, as it was mentioned in section 3.4, the PGMRES solver is utilized.

www.intechopen.com
Hydrodynamic Loads Computation Using the Smoothed Particle Methods 57

Pressure field in the original SPH method
In the SPH method barotropic condition for pressure pp ()   is supposed .For the first time Monaghan (Monaghan et al., 1994) applied equation for pressure in the Theta form: where Bg H 200 /    -gravitational constant,  -density, 0  -initial density, H -initial height of fluid, 7   .
Monaghan applied this equation for free surface flow simulations, such as breaking dam problems.But research of the calculation of pressure by ( 27) shows that pressure field in fluid has a significant oscillations.
To reduce pressure oscillations we smooth density field.For free surface problems in the case of the system being at rest under the action of gravity force at the initial time the hydrostatic pressure distribution is true: p gH y 00 () Then we can define the corrected value for the initial density from equation of state ( 27): Besides in time integration scheme for density computation the equation for density smoothing is added based on the formula (9) following Chen's idea (Chen et al., 2001): Thus the pressure at solid boundary particles is calculated using the values of the pressure at neighbouring fluid particles by formula (9).

Hydrodynamic loads
Hydrodynamic loads onto the solid boundary  is the integral characteristic of the wave pressure.Here the following formula is used (Afanasiev & Berezin, 2004): where p(0) is initial pressure and p T () is the pressure at any other moment T . www.intechopen.com Hydrodynamics -Optimizing Methods and Tools 58 In the numerical computations the value of the integral ( 31) is estimated by the formula: where b P is a set of solid boundary particles.

Nearest neighbour search
In numerical simulations using the smoothed particle methods it is necessary to determine for every particle j its interacting particles, as all physical characteristics of the fluid are estimated over the values at neighbouring particles according to the formula (9).For each fluid particle j its smoothing length j h is set, determining the radius of interaction with neighbours.As it is clear from section 3.1 in smoothed particle methods if particle i interacts with particle j then particle j interacts with particle i too, so forming the interacting pair.Thus it is necessary to solve a geometrical problem of determination of points which are in the circle of radius kh with the center at the point j (fig. 1 a The idea of the method consists in construction of a grid on each time step which fully covers the problem domain.The linear size of grid cells is constant and equals to: where 0   and 1   .At next step for each particle its belonging to one of the cells of a grid is defined.Then nearest neighbours for particle j are determined using direct search algorithm but only within the adjacent cells (fig. 1 b).
In fig. 2 the results of testing the speed of both algorithms are presented (X-axis corresponds to total number of particles and Y-axis corresponds to full time search procedure).
Test calculations were carried out on uniprocessor system: AMD Athlon 2000+, 512 Mb RAM.Time of nearest neighbour particle search depending on number of particles for 1000 time steps was measured.It can be noted that grid algorithm is very efficient and, for example, calculations with 8000 fluid particles gives acceleration of about 100.For the grid algorithm it is shown that its analytic time complexity is about ON () operations (Afanasiev et al., 2008) that agrees well with obtained numerical data (see fig. 2 b).

Poiseuille flow
This problem is one of the classical tests for viscous fluid flows, because of well-known analytical solution for velocity profile.Here two-dimensional non-stationary viscous fluid flow between two parallel solid walls is considered.Initially the fluid in the infinite channel, bounded with solid walls Г 2 and Г 4 , is at rest.Motion of fluid particles occurs in rectangular domain  , representing the infinite channel, due to difference of pressure at opposite open boundaries Г 1 and Г 3 (fig.3).On horizontal solid walls Г 2 and Г 4 the slip condition is set (the zero-valued velocity vector).
where in out P P , -the pressure at Г 1 and Г 3 accordingly; , L ,   are the density, dynamic viscosity and the channel length, H is the height of the channel.The infinity of the channel is simulated by cyclic returning of particles, passed through the right open boundary Г 3 , on left boundary Г 1 with the obtained physical characteristics.Pressure difference is simulated by the horizontal volumetric force F, directed from Г 1 to Г 3 : The analytical solution of above ordinary differential equation with slip boundary conditions takes the form (Leonardo et al., 2003): where dH /2  is the half-height of the channel, /    is the kinematic viscosity and the first term in the right hand side is the stationary velocity in the channel when t .

Laminar fluid flow along the infinite inclined plane
The problem is of special interest because it is one of few problems for viscous free surface flows, that have an analytic solution.The problem domain is shown in fig. 5 a.The fluid flow takes place in a rectangular infinite region  , bounded with solid wall Г 1 inclined at an angle  to the horizontal surface.Г 2 is free surface and initially fluid flow is at rest.
Fluid flow occurs under gravity force, directed vertically to the horizontal surface.On solid boundary Г 1 the slip condition is set.
The formulation can be simplified by performing rotation of the coordinate axes by angle  so that the X -axis coincides with the horizontal surface.Considering that the velocity of the fluid depends only on the vertical coordinate y : x vvy ()  , the action of gravity can be replaced by volumetric horizontal force F , which is the projection of gravity onto X -axis.Thus, the problem domain is changed to shown in fig. 5 As was mentioned above the problem has stationary analytic solution, that has the following form (Slezkin, 1995):

Droplet problem
At initial moment of time the problem domain  is a circle of incompressible fluid with radius r 1  and with center located at the origin of the coordinates (Ovsyannikov, 1967).Deformation of a circle into ellipse with semi-axes at () (along y 0  ) and bt () (along x 0  ) is initiated by the non-zero velocity field: Incompressibility is provided by constancy of ellipse's square for any moment of time, that is at bt () () 1  .So at any moment of time the form of ellipse is described with the following equation: x at y at where at () is taking from the system of ordinary differential equations: with appropriate initial conditions: In simulations using the ISPH method parameter  (in free surface detecting algorithm, see section 3.5) and parameter hd x  were varied.The best results were obtained for hd x 11   and 0 75   .The results are provided for these values of mentioned parameters.
Fourth-order Runge-Kutta method was used for obtaining sample results.The results of numerical simulations are provided only for the ISPH method as the results of the original SPH are very similar (Afanasiev et al., 2006).

Conclusion
As it follows from the results of simulation of dam breaking problem using SPH and ISPH, the both methods demonstrate good results in pressure field calculation.The time charts of hydrodynamic loads on the solid walls of the basin show good agreement, what proves their correctness.It can be concluded that the smoothed particle methods allow obtaining correct hydrodynamic loads in the case of large deformations of free surface and can be used for problems of that type.As some differences between hydrodynamic loads obtained on one side by SPH and on the other side by the ISPH method is observed, explanation of this fact becomes a subject for future work.It is also planned to analyse the influence of turbulent forces onto the hydrodynamic loads.
at solid boundary particles can be determined out of the following expression:

Fig. 1 .
Fig. 1.Nearest neighbour search: a) search area, b) cells for search Direct search algorithm has time complexity about ON 2 ( ) operations for procedure of determination of all interacting pairs, where N is the total number of particles in problem domain.Here the efficient algorithm, based on rectangular grid construction is implemented.The idea of the method consists in construction of a grid on each time step which fully covers the problem domain.The linear size of grid cells is constant and equals to:

Fig. 3 .
Fig. 3. Problem domain for Poiseuille flow Within the problem domain  the fluid motion is described with the simplified momentum equation: Fig. 4. Velocity profile for Poiseuille flow: a) t s 0.02  , b) t s 0.6 

Fig. 5 .
Fig. 5. Problem domain: a) initial, b) simplified The infinity of the channel is modeled by the algorithm described for Poiseuille flow in section 6.1.Provided x Fg sin  , equation of motion with slip boundary conditions is written as: v g y 2

Fig. 6 .
Fig. 6.Velocity profiles for laminar flow along the incline plane: a) t s 0.075  , b) t s 4 

Table 1 .
Numerical errors by SPH and ISPH for different sets of particles

Table 2 .
Numerical errors by SPH and ISPH for different sets of particles