- Email: [email protected]
A meshless Hermite-Cloud method for nonlinear fluid-structure analysis of near-bed submarine pipelines under current
Engineering Structures 26 (2004) 531–542 www.elsevier.com/locate/engstruct
A meshless Hermite-Cloud method for nonlinear fluid-structure analysis of near-bed submarine pipelines under current Hua Li a,, J.Q. Cheng a, T.Y. Ng a,b, Jun Chen a, K.Y. Lam a a
b
Institute of High Performance Computing, National University of Singapore, 1 Science Park Road, #01-01, The Capricorn, Singapore Science Park II, Singapore 117528, Singapore School of Mechanical & Production Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore Received 28 June 2003; received in revised form 1 December 2003; accepted 1 December 2003
Abstract Through higher-order Hermite construction of the interpolation functions, a novel true meshless numerical technique—the Hermite-Cloud method is developed for the analysis of nonlinear fluid-structure interaction of near-bed submarine pipelines under a current. Following the development of the discrete point collocation technique throughout the computational domains, the present method constructs the approximations of both the unknown functions and their first-order derivatives based on the classical reproducing kernel particle method, except that a fixed reproducing kernel approximation is employed instead. For partial differential equations with boundary conditions, certain auxiliary conditions are required due to the use of the Hermite theorem. After validation of the presently developed Hermite-Cloud method by several two-dimensional elasticity examples, the method is used for the failure analysis of a near-bed submarine pipeline under a horizontal current with consideration of a nonlinear fluid–structure interaction effect. The numerical result indicates that the distance from the pipelines to seabed and the current velocity have remarkable influences on the deformation behavior of the pipelines. It also shows that various critical current velocities corresponding to different failure patterns exist, and the presently simulated relation of the critical velocities agrees with the results in the previous publications. # 2003 Elsevier Ltd. All rights reserved. Keywords: Meshless method; Hermite-Cloud; RKPM; Point collocation; Nonlinear fluid–structure interaction; Submarine pipelines
1. Introduction The development of the finite element method (FEM) has formed the basis for the development of commercial computational mechanics based software for the modeling and simulation of various mechanical computer aided engineering (MCAE) problems. In this endeavor however, the requirement of mesh generation in the FEM results in certain disadvantageous features, such as iterative remeshing for the tracking of dynamic processes in large deformation problems and the requirement for large storage and memory due to the large number of element nodes involved in FEM discretization. To overcome these FEM associated deficiencies, various meshless numerical techniques [1] have
Corresponding author. Tel.: +65-64191561; fax: +65-64191380. E-mail address: [email protected] (H. Li).
0141-0296/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2003.12.005
been recently developed, for example, the backgroundmesh type of approaches including the Galerkin-based techniques [2–12] and the direct domain discretization such as point collocation-based techniques [13–18]. Some interesting work should also be noted [19–23]. To explore more efficient collocation-based techniques, this paper develops a novel true meshless technique — the Hermite-Cloud method, which is based on the Hermite interpolation theorem and point collocation as well as reproducing kernels, and represents a significant development in reproducing kernel techniques. It employs the Hermite theorem for the construction of the interpolation functions, where the shape functions are constructed to correspond respectively to the unknown functions and the first-order derivatives of these unknown functions. It is formulated on the basis of the classical reproducing kernel particle method (RKPM) but fixed kernels are used here
532
H. Li et al. / Engineering Structures 26 (2004) 531–542
instead [17]. Furthermore, for a given set of the partial differential equations (PDEs) with the Dirichlet and/or Neumann boundary conditions, certain differential-type auxiliary conditions are derived based on the Hermite theorem to generate a complete set of equations for the governing partial differential boundary value (PDBV) problem. Following the scattering of a set of points in the computational domain including the edges, the point collocation technique is applied for the PDBV discretization. The approximate solutions of both the unknown functions and their first-order derivatives are thus expressed in terms of both the shape functions and unknown point values leading to a complete set of discrete algebraic equations being described. Finally, they are solved with respect to the unknown point values and the numerical solutions of the PDBV problem can be computed in a straightforward manner. To validate the developed Hermite-Cloud method, several two-dimensional (2-D) elasticity examples, including a higher-order patch subjected to a uniform unity-intensity load, cantilever beams subjected to a linearly varying axial and a shear end loads respectively, and a coupled thermo-elasticity problem, are analyzed numerically and the results are in good agreements with either exact results or existing numerical meshless solutions. For the application of the present Hermite-Cloud method in the submarine engineering, a nonlinear fluid–structure interaction of near-bed submarine pipelines under a current is studied here. The behavior of near-seabed submarine pipelines is greatly affected by the environment factors, such as ocean currents, earthquakes and underwater explosions. This paper investigates the influence of current on the submarine pipelines. When the distance between a supported submarine pipeline and the seabed is very large, the current around the pipeline is hardly affected by the seabed and the flow is symmetrical to the pipeline. In the vertical direction, there is no net force applied on the pipeline. However, when the distance comes to be small, the velocity of current between the pipeline and seabed becomes higher, compared with that above the pipeline. Then based on Bernoulli’s equation, the current pressure above the pipeline is higher than that between the pipeline and seabed. A downward net force caused by the pressure difference is thus produced and the pipeline is expected to deform under the external load. The deflection of pipelines enlarges with the increase of current velocity. When the velocity reaches a critical value, the pipeline will fail to work. In this paper, the two types of failure patterns are considered. One is caused by instability, in which the pipeline will rest on the seabed completely, and the other results from the limit of material properties, in which the stress/deflection caused by the external load reaches the allowable stress/deformation of the constitutive
materials. The variation of critical velocities with the distance from pipelines to seabed in various failure patterns is obtained. The distribution of the pipeline stress and the relation between the pipeline deflection and current velocity are simulated, respectively. The discussion on these computed results is also made in detail, which indicates a good agreement with the published work [19].
2. Development of meshless Hermite-Cloud method 2.1. Theoretical formulation As mentioned earlier, as one in the class of mathematical operators that can reproduce a function by integration transform over a computational domain, the reproducing kernel method for a 2-D real function f ðx;yÞ in a domain X is expressed as follows: ð f ðx;yÞ ¼ Uðx p;y qÞf ðp;qÞdpdq; ð1Þ X
where Uðx;yÞ is a real window function. According to the definition of the reproducing kernel method, the ideal window functions should be orthogonal, and its integration over X should be unity, so as to exactly reproduce f ðx;yÞ. The selection of window functions differentiates the various reproducing kernel methods. However, it is not often easy to select suitable functions as the window function must satisfy both the conditions to exactly reproduce the unknown function f ðxÞ. In the classical RKPM [3,4], an approximate solution is constructed by introducing a window function Uðx;yÞ as the product of a correction function Cðx;y;p;qÞ and a kernel function Kðx;yÞ, that is, Uðx p;y qÞ ¼ Cðx;y;p;qÞKðx p;y qÞ:
ð2Þ
Substituting Eq. (2) into Eq. (1), the approximation ~f ðx;yÞ of the unknown function f ðx;yÞ can be represented as follows ð ~f ðx;yÞ ¼ Cðx;y;p;qÞKðx p;y qÞf ðp;qÞdpdq: ð3Þ X
Further to that mentioned above, the selection of kernel functions also differentiates the various reproducing kernel methods. If the fixed kernel is used as the present kernel function Kðx;yÞ, Eq. (3) can be rewritten as ð ~f ðx;yÞ ¼ Cðx;y;p;qÞKðxk p;yk qÞf ðp;qÞdpdq ð4Þ X
where the point ðxk ;yk Þ is the center point of the fixed kernel corresponding to the kernel function Kðxk p;yk qÞ. Further, depending on different PDBV problems, the kernel function may be constructed by
H. Li et al. / Engineering Structures 26 (2004) 531–542
different forms of weighted window functions, such as the Gaussian functions, spline functions or radial basis functions. In the present Hermite-Cloud method, we employ a cubic spline function to construct the kernel function Kðxk p;yk qÞ ¼ W ððxk pÞ=DxÞW ððyk qÞ=DyÞ=ðDxDyÞ;
conditions, Eq. (9), is rewritten in discrete form as NT X Cðx;y;pn ;qn ÞKðxk pn ;yk qn Þbi ðpn ;qn ÞDSn
bi ðx;yÞ ¼
n¼1 NT X ¼ Bðpn ;qn ÞC ðx;yÞKðxk pn ;yk qn Þ n¼1
ð5Þ
where W ðzÞis a cubic-spline window function and is given as 8 0 j zj 2 > > > > < ð2 jzjÞ3 1 j zj 2 ð6Þ W ðzÞ ¼ 6 > > 2 j z j > > j zj 1 : z2 1 3 2 in which z ¼ ðxk pÞ=Dx or z ¼ ðyk qÞ=Dy. Dx and Dy denote the cloud size with the center point ðxk ;yk Þ. These are adjusted according to the point distribution and accuracy requirement, due to the consistency conditions, of the reproducing kernel approach. It is well-known mathematically that a continuous function can be represented as a sum of the linearly independent functions. As a result, the correction function Cðx;y;p;qÞ in Eq. (4) may be expressed as
bi ðpn ;qn ÞDSn
BT ðx;yÞ ¼ Aðxk ;yk ÞC ðx; yÞ
C ðx; yÞ ¼ A1 ðxk ; yk ÞBT ðx; yÞ;
ð12Þ
where Aðxk ;yk Þ is a symmetric constant matrix at the fixed-cloud center point ðxk ; yk Þ, Aij ðxk ;yk Þ ¼
NT X
bi ðpn ;qn ÞKðxk pn ;yk qn Þ
bj ðpn ;qn ÞDSn
ði; j ¼ 1; 2; . . . ; bÞ:
ð13Þ
By substituting Eqs. (5), (7) and (12) into Eq. (4), the approximate expression f~ðx;yÞ of the unknown real function f ðx;yÞ can be obtained as ð ~ f ðx;yÞ ¼ Bðp;qÞC ðx;yÞKðxk p;yk qÞf ðp;qÞdpdq ðX ¼ Bðp;qÞA1 ðxk ;yk ÞBT ðx;yÞ X
Kðxk p;yk qÞf ðp;qÞdpdq;
ð14Þ
which can further be discretized as ð8Þ
X
ði ¼ 1; 2; . . . ; bÞ:
ð11Þ
we obtain
ð7Þ
In the coefficient vector CT ðx;yÞ=fc1 ;c2 ; . . . ;cb g of Eq. (7), the coefficients ci ði ¼ 1; 2; . . . ; bÞ are unknown and can be determined through the following consistency conditions of the reproducing kernel techniques, ð bi ðx;yÞ ¼ Cðx;y;p;qÞKðxk p;yk qÞ bi ðp;qÞdpdq
ð10Þ
n¼1
where C ðx;yÞ is a bth-order column coefficient vector and Bðp;qÞ a bth-order row basis-function vector, in which bi ðp;qÞ ði ¼ 1; 2; . . . ; bÞ are linearly independent basis functions (b denotes the degree of the polynomials of the basis function). Usually their definitions depend on the PDBV problem to be solved. For example, for a 1-D PDE system, we may take BðpÞ ¼ f1; p; p2 ; . . . ; pb1 g. For a 2-D second-order PDE system, a quadratic basis function Bðp;qÞ is required Bðp;qÞ ¼ fb1 ðp;qÞ;b2 ðp;qÞ; . . . ;bb ðp;qÞg ¼ f1;p;q;p2 ;pq;q2 g ðb ¼ 6Þ:
ði ¼ 1;2; . . . ;bÞ;
where NT is the total number of scattered points covering the interior domain X and the surrounding edges. The subscript n represents the nth scattered point. DSn is defined as the cloud area of the nth point. It is obvious that Eq. (10) is a set of linear algebraic equations with respect to the coefficients ci ði ¼ 1; 2; . . . ; bÞ. Thus, if it is rewritten in the following matrix form with respect to the coefficient vector CT ðx;yÞ ¼ fc1 ; c2 ; . . . ; cb g,
Cðx;y;p;qÞ ¼ Bðp;qÞC ðx;yÞ ¼ fb1 ðp;qÞ;b2 ðp;qÞ; . . . ;bb ðp;qÞgfc1 ;c2 ; . . . ;cb gT ;
533
ð9Þ
If a set of points is scattered arbitrarily in the integral domain X and along its edges, the point collocation technique can be applied for PDBV discretization. By substituting Eq. (7) into Eq. (9), the consistency
f~ðx;yÞ ¼
NT X ðBðpn ;qn ÞA1 ðxk ;yk ÞBT ðx;yÞ n¼1
Kðxk pn ;yk qn ÞDSn Þfn ¼
NT X Nn ðx;yÞfn ;
ð15Þ
n¼1
where fn is defined as the point value of the nth point. Nn ðx;yÞ are the shape functions that consist of the basis functions and are simply polynomials in x and y. As such, any derivatives of the shape functions can be obtained easily by differentiation of the basis function vector Bðx;yÞ. Furthermore, it is evident that the present shape functions Nn ðx;yÞ satisfy the consistency
534
H. Li et al. / Engineering Structures 26 (2004) 531–542
conditions, Eqs. (9) or (10), for all independent basis functions bi ðx;yÞ ði¼ 1; 2; . . . ; bÞ. In particular, when b1 ðx;yÞ ¼ 1:0 (i¼ 1), b2 ðx;yÞ ¼ x (i¼ 2) and b3 ðx;yÞ ¼ y (i¼ 3) are taken for the discretized consistency conditions Eq. (10), we have 1:0 ¼
NT X
Nn ðx;yÞ; x ¼
n¼1
NT X Nn ðx;yÞxn ; n¼1
NT X Nn ðx;yÞyn : y¼
ð16Þ
n¼1
According to the Hermite theorem, the first-order derivatives of the unknown real function f ðx;yÞ with respect to the variables x and y are gx ðx;yÞ ¼
@f ðx;yÞ @f ðx;yÞ ; gy ðx;yÞ ¼ : @x @y
ð17Þ
n¼1
These result in additional unknown functions. If they are also discretized by imposing Eq. (15) in a similar manner, their approximate expressions can be written as ~ gx ðx;yÞ ¼
NS X
Mm ðx;yÞgxm
differential governing equations, consists of all three unknown point-value vectors ffn g1NT , fgxm g1NS and fgym g1NS , thus making it possible to directly compute the approximate solutions of both the unknown functions and their first-order derivatives. In addition, the Hermite-based interpolation approximation, Eq. (20), is also used for the construction of the following necessary auxiliary conditions. Due to the inclusion of additional unknown functions gx ðx;yÞ and gy ðx;yÞ, certain auxiliary conditions are required to generate a complete set of PDBV equations. Based on Eq. (17) for the defined gx ðx;yÞ and gy ðx;yÞ, the following auxiliary conditions are developed naturally by considering Eqs. (16), (18) and (20) ! NS NT NT X X X Nn;x ðx;yÞfn ðNn;x ðx;yÞxn Þ Mm ðx;yÞgxm
ð18Þ
m¼1
ð19Þ
where NS ð NT Þ is the total number of scattered points, and gxm and gym are the unknown point values for the mth point. Mm ðx;yÞ are the shape functions corresponding to the unknown first-order differential functions gx ðx;yÞ and gy ðx;yÞ, but their basis functions are now Bg ðx;yÞ 2 Rb1 . Based on the Hermite interpolation theorem, a true meshless approximation f~ðx;yÞ of the unknown real function f ðx;yÞ can now be finally constructed in the following form f~ðx;yÞ ¼
NT X Nn ðx;yÞfn n¼1
þ
NS X m¼1
þ
NS X m¼1
x
NT X
ðNn;x ðx;yÞyn Þ Mm ðx;yÞgym ¼ 0
ð21Þ
m¼1 NS NT X X m¼1
n¼1
!
ðNn;y ðx;yÞxn Þ Mm ðx;yÞgxm ¼ 0;
ð22Þ
n¼1
where the variable in the subscript after a comma indicates partial differentiation with respect to that variable. These auxiliary conditions are a necessary requirement in the implementation of the proposed Hermite-Cloud method. The above formulation thus defines the HermiteCloud method. In summary, it is based on the Hermite interpolation theorem, and consists of the approximation f~ðx;yÞ of the unknown function f ðx;yÞ, the approximations ~gx ðx;yÞ and ~gy ðx;yÞ of the first-order derivatives gx ðx;yÞ and gy ðx;yÞ, and further combines the necessary auxiliary conditions. 2.2. Numerical implementation
!
In general, engineering PDBV problems can be written as
Nn ðx;yÞxn Mm ðx;yÞgxm
n¼1
! NT X y Nn ðx;yÞyn Mm ðx;yÞgym :
!
n¼1
n¼1
m¼1
n¼1
! NS NT NT X X X Nn;y ðx;yÞfn ðNn;y ðx;yÞyn Þ Mm ðx;yÞgym
m¼1 NS X ~ gy ðx;yÞ ¼ Mm ðx;yÞgym ;
m¼1 NS NT X X
ð20Þ
n¼1
The presently constructed Hermite-based interpolation approximation has many computational advantages. Most notably, the computational accuracy at scattered discrete points in the domain is much refined not only for the approximate functions, but also for their firstorder derivatives. This can be easily understood from the subsequent formulation, Eqs. (29) and (32), where the final unknown vector of the set of the algebraic governing equations, derived by the discretized partial
Lf ðx;yÞ ¼ Pðx;yÞ PDEs in interior domain X
ð23Þ
f ðx;yÞ ¼ Qðx;yÞ Dirichlet boundary condition on CD ð24Þ @f ðx;yÞ ¼ Rðx;yÞ Neumann boundary condition on CN @n ð25Þ where L is a differential operator and f ðx;yÞ an unknown real function. By the point collocation technique and taking f~ðx;yÞ as the approximation of f ðx;yÞ, the problem is discretized and expressed
H. Li et al. / Engineering Structures 26 (2004) 531–542
535
approximately by Lf~ðxi ;yi Þ ¼ Pðxi ;yi Þ i ¼ 1;2; . . . ;NX
ð26Þ
f~ðxi ;yi Þ ¼ Qðxi ;yi Þ i ¼ 1;2; . . . ;ND
ð27Þ
@ f~ðxi ;yi Þ ¼ Rðxi ;yi Þ i ¼ 1;2; . . . ;NN ; @n
ð28Þ
where NX, ND and NN are the numbers of scattered points in the interior domain, and along the Dirichlet and Neumann edges, respectively, and the total number of scattered points is thus NT ¼ ðNX þND þNN Þ. Substituting the approximations, Eqs. (18) and (20), into Eqs. (26) and (28), and further combining the auxiliary conditions of Eqs. (21) and (22), followed by rearrangement of the resulting equations, a set of discrete algebraic governing equations with respect to the unknown point values fi , gxi and gyi , is obtained and can be written in matrix form as ½Hij ðNT þ2NS ÞðNT þ2NS Þ fFi gðNT þ2NS Þ1 ¼ fdi gðNT þ2NS Þ1 ; ð29Þ where fdi g and fFi g are ðNT þ 2NS Þ-order column vectors, with fFi gðNT þ2NS Þ1 ¼ fffi g1NT ; fgxi g1NS ; fgyi g1NS gT
ð30Þ
fdi gðNT þ2NS Þ1 ¼ ffPðxi ;yi Þg1NX ; fQðxi ;yi Þg1ND ; fRðxi ;yi Þg1NN ; f0g12NS gT
ð31Þ
and [Hij] is a ðNT þ2NS Þ ðNT þ2NS Þ coefficient square matrix 2
6 LNj ðxi ;yi Þ NX NT 6 6 6 Nj ðxi ;yi Þ ND NT 6 6 ½0 NN NT ½Hij ¼ 6 6 6 6 Nj;x ðxi ;yi Þ NS NT 6 6 4 Nj;y ðxi ;yi Þ NS NT
Fig. 1. Geometry and point distribution for the higher-order patch subjected to a uniform unidirectional stress of unity magnitude.
3. Validation study of Hermite-Cloud method To validate the numerical stability and accuracy of the present Hermite-Cloud method, numerical comparisons are made in this section for several 2-D planestress elasticity problems. These include a higher-order patch subjected to a uniform unity-intensity unidirectional stress, cantilever beams subjected to linearly varying axial loads or shear end loads, and a 2-D thermo-elasticity analysis. In order to measure the numerical accuracy of the present Hermite-Cloud method, a refined global error n is defined [5] for the numerical comparison: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u NT X 1 u t 1 ð33Þ ðf~i fi Þ2 : n¼ jfmax j NT i¼1 3.1. Higher-order patch For an isotropic higher-order patch subjected to uniform unidirectional stress of unity magnitude, as
NT P L xi Nn ðxi ;yi Þxn Mj ðxi ;yi Þ n¼1
n¼1
n¼1
NX NS
½0ND NS Mj ðxi ;yi Þ NN NS N PT Nn;x ðxi ;yi Þxn Mj ðxi ;yi Þ NS NS n¼1 NT P Nn;x ðxi ;yi Þxn Mj ðxi ;yi Þ NS NS
3
NT P L yi Nn ðxi ;yi Þyn Mj ðxi ;yi Þ NX NS
½0ND NS Mj ðxi ;yi Þ NN NS N PT Nn;x ðxi ;yi Þyn Mj ðxi ;yi Þ NS NS n¼1 NT P Nn;x ðxi ;yi Þyn Mj ðxi ;yi Þ NS NS
7 7 7 7 7 7 7: 7 7 7 7 7 5
n¼1
ð32Þ
Solving numerically the complete set of linear algebraic equations, Eq. (29), (NT þ NS ) point values {Fi}, consisting of the NT point values {Fi} and 2NS point values {gxi} and {gyi}, are obtained and the approximate solution f~ðx;yÞ and its first-order derivatives ~ gx ðx;yÞ and ~ gy ðx;yÞ of the PDBV problem can be computed through Eqs. (20), (18) and (19).
shown in Fig. 1, the plane-stress equilibrium governing equations is given in terms of the displacements u and v as follows: 2 E @ u 1 l @2u 1 þ l @2v þ þ ð34Þ þ Xx ¼ 0 1 l2 @x2 2 @y2 2 @x@y 2 E @ v 1 l @2v 1 þ l @2u þ þ ð35Þ þ Xy ¼ 0; 1 l2 @y2 2 @x2 2 @x@y
536
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 2. Numerical comparison for the higher-order patch subjected to a uniform unidirectional stress of unity magnitude: (a) displacement u, (b) displacement v.
where E is the elastic modulus and l the Poisson’s ratio. Xx and Xy denote the body forces in the x and y directions, respectively. If the material constants E ¼ 1 and l ¼ 0:25 are taken, the exact solution of this problem is uðx;yÞ ¼ x;
vðx;yÞ ¼ 0:25y:
ð36Þ
After implementation of the presently developed Hermite-Cloud method, the displacements u and v are obtained numerically and their profiles are plotted in Fig. 2(a,b). It is observed that, for a point distribution 53, comparisons with the exaction solution reveal global errors of n ¼ 3:19 107 for the displacement u, and n ¼ 3:32 107 for the displacement v. It is clearly evident that the present Hermite-Cloud method can achieve almost exact results at the scattered discrete points even with a coarse 53 point distribution for the present patch test problem. This example clearly demonstrates the efficiency and ease of implementation of the present Hermite-Cloud method are much enhanced over that of existing generic meshless techniques, such as the EFG method, since the present Hermite-Cloud method does not require a background mesh and has better consistency characteristics in the construction of the shape functions.
The exact solution for the case shown in Fig. 3 is given by uðx;yÞ ¼ xy; vðx;yÞ ¼
ð4x2 þ y2 Þ : 8
ð37Þ
Applying the Hermite-Cloud method for this elastic analysis, the numerical displacements u and v are computed. The displacement profiles are shown in Fig. 4(a,b). It is again observed that, for the same coarse 53 point distribution, comparison with the exact solution reveals global errors of n¼ 1:61 106 for the displacement u and n¼ 3:86 106 for the displacement v. In addition, the exact solution for the cantilever beam subjected to a shear end load, as shown in Fig. 5, is given by uðx;yÞ ¼
Pðy 0:5DÞ ð6L 3xÞx þ ð2 þ lÞðy2 2DyÞ 6EI ð38Þ
P 3lðy2 2Dy þ 0:5D2 ÞðL xÞ 6EI þð1 þ 1:25lÞD2 x þ ð3L xÞx2 ;
vðx;yÞ ¼
ð39Þ
3.2. Cantilever beam subjected to various loadings Here, a 2-D plane-stress cantilever beam is considered with the same equilibrium governing equations as Eqs. (34) and (35). It is subjected to a linearly varying axial end load as shown in Fig. 3, and subjected to a shear end load as shown in Fig. 5, respectively.
Fig. 3. Geometry and point distribution for a cantilever beam subjected to a linearly varying axial end load.
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 4.
537
Numerical comparison for the cantilever beam subjected to a linearly varying axial end load: (a) displacement u, (b) displacement v.
where I ¼ D3 =12 and D and L are respectively the height and length of the cantilever beam. The displacement fields are solved numerically using the Hermite-Cloud method, and the computed displacement profiles are shown in Fig. 6(a,b). Further, in order to examine the convergence of the Hermite-Cloud method, Table 1 compares the convergence characteristics between the present Hermite-Cloud method and
the meshless hM-DOR method [18]. It is observed from the table that global errors of both the displacement u and first-order derivative u,x in the present HermiteCloud method decrease monotonically with increasing number of scattered points, and the global-errors are lower than those of the hM-DOR method [18]. These numerical results thus show that the computational accuracy of the Hermite-Cloud method at scattered discrete points in the domain is much refined not only for the approximate solutions, but also for the first-order derivatives of these solutions. 3.3. Coupled thermo-elasticity
Fig. 5. Geometry and point distribution for a cantilever beam subjected to a shear end load.
Fig. 6.
As shown in Fig. 7, a plane-stress analysis of coupled thermo-elasticity is carried out by the developed Hermite-Cloud method. The partial differential governing equations that couple the heat conduction equation with the generic thermo-stress elastic equations, given
Numerical comparison for the cantilever beam subjected to a shear end load: (a) displacement u, (b) displacement v.
538
H. Li et al. / Engineering Structures 26 (2004) 531–542
Table 1 Numerical comparisons between the present Hermite-Cloud method and the hM-DOR [18] method for a cantilever beam subjected to a shear end load (E ¼ 3:0 107 , D ¼ 1, L ¼ 8 and l ¼ 0:25) Points distribution(NxNy)
511 1111 2111 6511 8111
Global error (n) for u
Global error (n) for u,x
hM-DOR [18] (%)
Hermite-Cloud (%)
hM-DOR [18] (%)
Hermite-Cloud (%)
35.9 17.9 6.16 4.28 4.01
35.9 17.7 5.92 1.95 0.384
30.9 15.0 5.38 3.92 3.90
30.9 14.9 5.19 2.13 1.92
in terms of the displacements u and v and temperature T, are @2T @2T þ 2 ¼h @x2 @y
ð40Þ
@2u 1 l @2u 1 þ l @2v @ðT T0 Þ ð1 þ lÞa ¼0 þ þ 2 2 @x 2 @y 2 @x@y @x ð41Þ @2v 1 l @2v 1 þ l @2u @ðT T0 Þ ð1 þ lÞa ¼ 0; þ þ @y2 2 @x2 2 @x@y @y ð42Þ where h is the heat source, T0 the initial reference temperature and a the coefficient of thermal expansion. The considered boundary conditions are uðx;yÞ ¼ 0; rxy ðx;yÞ ¼ 0; Tx ðx;yÞ ¼ 0; at x ¼ 0 and 1;
ð43Þ
vðx;yÞ ¼ 0; rxy ðx;yÞ ¼ 0; Tðx;yÞ ¼ 0; at y ¼ 0
ð44Þ
rxy ðx;yÞ ¼ 0; ryy ðx;yÞ ¼ 0; Tðx;yÞ ¼ 0; at y ¼ 1:
ð45Þ
7
If we take T0 ¼ 0, h ¼ 20, E ¼ 3:0 10 , a ¼ 0:001 and l ¼ 0:25, the exact solutions of this coupled thermo-elasticity problem are Tðx;yÞ ¼ 10y2 ; uðx;yÞ ¼ 0; vðx;yÞ ¼ 10að1 þ lÞy3 =3 ¼ 0:0125y3 =3:
Fig. 7. Geometry of the 2-D thermoelasticity problem.
Again by using the Hermite-Cloud method, the variable fields Tðx;yÞ, uðx;yÞ and vðx;yÞ and their firstorder derivatives are obtained numerically, with a regular 55 point distribution. The numerical results show that the standard error of the displacement u is of the order of 1011 and for the temperature Tðx;yÞ, the global error is less than 1.0107. A probable reason for the excellent results is that the exact solutions of u and T can be found in the basis functions. However, the same cannot be said for the numerical results of the displacement v since its exact cubic solution is not found in any of the basis functions. Compared with exact solution, Fig. 8 presents the accuracy variation of the displacement v with the increase of the regularly scattered points, for NT ðtotal point numberÞ ¼ 25 ð5 5Þ; 121 ð11 11Þ and 441 ð21 21Þ. It is evident that the computed accuracy of displacement v is significantly refined with the increase of the number of scattered points. Fig. 9 compares the convergence characteristics of the present displacement vðx;yÞ and the first-derivatives v;y ðx;yÞ, with that obtained by the meshless Finite Cloud method [17]. Present convergence rates are 2.12 for vðx;yÞ and 1.89 for v;y ðx;yÞwhile that for the Finite Cloud method are 2.2 for vðx;yÞ and 1.95 for v;y ðx;yÞ. It is therefore clear that the present Hermite-Could method attains better
ð46Þ
Fig. 8. Variation of the numerical displacement v with the point distribution density for the thermo-elasticity case.
H. Li et al. / Engineering Structures 26 (2004) 531–542
539
ditions for the above problem are given as
Fig. 9. Convergence comparison between the present HermiteCloud method and Finite- Cloud method [17] for the thermo-elasticity case (n, global error; h, point distance).
computational accuracy for both the field-variable solution and the first-order derivative than the Finite Cloud method.
; on CW ; wðx0 Þ ¼ w
ð49Þ
hðx0 Þ ¼ h; on Ch @h ; on CM Mðx0 Þ ¼ EI ¼M @x x¼x0 @w ; on CV V ðx0 Þ ¼ GAks h þ ¼V @x x¼x0
ð50Þ ð51Þ ð52Þ
in which Cw , Ch , CM and CV are the boundaries where w, h, M and V satisfy, respectively. In Eq. (47), f is the fluid force caused by the current pressure difference. Although it has an analytical solution, it converges slowly. In order to speed up the computation, Lam et al. used the boundary element method (BEM) to construct the approximation of the exact solution. They are in a very good agreement with each other [24]. The expression of fluid force f(x) is given as f ðxÞ ¼
1 qAU02 cðdÞ; 2
ð53Þ
where qis the mass density, U0 the current velocity, c and d are dimensionless coefficient and have the relation as 4. Numerical study of a near-bed submarine pipeline under current When a submarine pipeline is placed near the seabed, the horizontal current is not symmetrical to the pipeline due to the existence of seabed. As a result, a downward external force is produced and applied on the pipeline, resulting in the deformation of pipeline. This problem can be simplified as a beam with different boundary support conditions. In the previously published work [24], Lam et al. studied the steady behavior of the pipeline under a steady nonlinear external force, in which the pipeline was simplified as a Bernoulli– Euler beam by the assumption that there is no transverse shear strain. However it is not applicable in the case that the ratio of length to diameter of the pipeline is not very large, where the effect of shear deformation should be included. In order to break through the limit, the Timoshenko beam theory [25] is introduced here to include the transverse shear deformation. Based on the Timoshenko beam theory, the general governing equations of the pipeline are written as @ @w GAks þh þf ¼0 ð47Þ @x @x @ @h @w EI þ h ¼ 0; ð48Þ GAks @x @x @x where w is the deflection of the pipeline, h the rotation, G the shear modulus, A the cross section area, ks the shear correction coefficient, E the elasticity modulus, I the moment of inertia. The general boundary con-
cðdÞ ¼
2:23d 2 þ 12:54d þ 0:02 0:77d 3 þ 0:44d 2 þ 0:02d
ð54Þ
in which d is defined as d ¼ ðD0 wðxÞ Rs Þ=ð2Rs Þ, and D0 is the distance between the central line of pipeline at initial status and the seabed, Rs the pipe outer radius. It is noted from the Eqs. (53) and (54) that the fluid force f(x) is a nonlinear term related to the deflection w(x). The governing Eqs. (47) and (48) for the pipeline deformation thus are the coupled nonlinear functions. It is also observed that, when the material and dimension of pipelines are given, the fluid force f(x) is a function of both the current velocity U0 and the gap D0. Namely these two parameters have significant influences on the deformation behavior of the near-bed submarine pipelines. In order to demonstrate the influences in detail, a practical example is designed here. A circular steel pipeline is placed near the seabed under a steady current, as shown in Fig. 10. Several physical and material parameters are given, the pipeline length L¼ 10 m, Rs ¼ 0:5 m, the thickness of steel pipeline ts ¼ 0:02 m, D0 ¼ 0:7 m, Es ¼ 2:11 1011 N=m2 , Poisson ratio m ¼ 0:3, qs ¼ 7800 kg=m3 , steel yield stress Ys ¼ 2:5 108 N=m2 , and allowable deflection/span ratio Yw ¼ 0:004. It is noted that the ratio L/Rs is not very large. As such, the pipeline can be simplified as a thick beam fixed at two ends. By the Timoshenko beam theory and HermiteCloud meshless method, numerical analysis is carried out and discussions are also presented in succession.
540
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 10. Schematic diagram of a submarine pipeline’s deformation under a current.
The variation of mid-point deflection wL=2 with the current velocity U0 is illustrated in Fig. 11. It is observed that, with the increase of current velocity, the deflection enlarges until a critical value Ucb. When the velocity is smaller than the critical value, the deflection converges with the iterations and the pipeline reaches the stable equilibrium. However, if the velocity is over the critical one, the deflection does not converge, which is defined as the instability, i.e. the pipeline has fallen into the seabed. It is the first failure pattern of pipelines—the instability failure, as mentioned above. The influence of the distance D0 between the pipeline and seabed under the critical current velocities caused by instability Ucb is revealed in Fig. 12. It is shown that Ucb is a monotonically increasing function with D0. As indicated in Eqs. (53) and (54), if other variables are fixed, the fluid force f(x) decreases with the enlargement of D0 and increases with the increasing U0. Therefore, in order to obtain the same critical fluid force causing the pipeline to lose the stability, when the gap D0 becomes large, it is reasonable to require a higher critical velocity Ucb. It is also expected that the critical velocity is infinite when the gap becomes infi-
nite, since the effect of seabed on the current is then zero and the flow symmetry is not changed at all. As well known, sometimes the failure due to material properties happens before the instability failure. Thus it is necessary to discuss material failure. The distribution of stress along the pipeline is plotted in Fig. 13,
Fig. 11. Variation of the deflection of mid-point of the pipeline with respect to the current velocity U0 (when D0 ¼ 0:7 m).
Fig. 13. Distribution of the stress along the pipeline (when U 0 ¼ 10 m=s and D0 ¼ 0:7 m).
Fig. 12. Effect of the gap D0 on the critical velocities Ucb of the instability failure.
H. Li et al. / Engineering Structures 26 (2004) 531–542
541
Fig. 14. Critical velocities of (a) strength failure, (b) deflection failure (when D0 ¼ 0:7 m).
in which the curve is found to be similar to that of a two-end fixed beam under uniform load, but actually they are different. For the beam under uniform load, the maximal absolute value of stresses (moments) along the beam occurs at the mid-point, yet on the contrary, the present maximal absolute value of stresses (moments) appears at the end points, as shown in Fig. 13. This phenomenon results from the different external load distribution. In this problem, the fluid force f(x) is a decreasing function of pipeline deflection w(x), and it is clear that the deflection of two ends is smaller than that of mid-point. This implies that the fluid force at two ends is larger than that of mid-point. Obviously it is different from the case of uniform load distribution. Fig. 14 presents the critical velocities caused by the two types of material failures—the strength failure and the deflection failure. As mentioned above, the endpoint stress and mid-point deflection firstly reach the
respective critical values along the whole pipeline. Thus their corresponding critical current velocities are defined as the critical values of the whole system. In the Fig. 14(a,b), both the end-point stress and midpoint deflection enlarge with the increase of current velocity U0, and when they reach the yield stress and allowable deflection, the corresponding critical velocity Ucs and Ucw are obtained. Similar to the instability failure, the distance between the pipeline and seabed D0 has an influence on the critical velocities of material failure. In order to compare these two failure patterns, the distributions of respective critical velocities (Ucb due to instability failure, Ucs due to strength failure and Ucw due to deflection failure) with respect to the gap D0 are synthetically plotted in Fig. 15. It is found that there are bifurcations at point A and B. In other words, beyond point A, the pipeline will not break down due to the stability loss, because the strength failure always appears before the instability failure. The results agree with Lam’s work [24], in which the pipeline is simplified as a Bernoulli–Euler beam. 5. Conclusions
Fig. 15. Comparison of distributions of respective critical velocity with respect to the gap D0 in various failure patterns.
This paper presents the development of a novel true meshless numerical technique, termed the HermiteCloud method. It employs the Hermite interpolation theorem to construct the interpolation functions and uses the point collocation technique for discretization of partial differential equations. The necessary auxiliary conditions are derived for the complete set of partial differential equations with mixed Dirichlet and Neumann boundary conditions. Various 2-D numerical studies have been conducted to validate the convergence characteristics and numerical accuracy of the Hermite-Cloud method and it is concluded that the developed Hermite-Cloud method is a very efficient
542
H. Li et al. / Engineering Structures 26 (2004) 531–542
numerical technique for solving general PDVPs. By the Hermite-Cloud method, the analysis of the nonlinear fluid-structure interaction of near-bed submarine pipelines under a current is completed. The numerical results illustrate the relation between the critical current velocities and the distance from pipelines to the seabed in different pipeline failure pattern. They are in good agreements with previous work.
References [1] Liu GR. Mesh free methods: moving beyond the finite element method. CRC Press LLC; 2003. [2] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. International Journal for Numerical Methods in Engineering 1994;37:229–56. [3] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. International Journal for Numerical Methods in Engineering 1995;20:1081–106. [4] Liu WK, Jun S, Li S, Adde J, Belytschko T. Reproducing kernel particle methods for structural dynamics. International Journal for Numerical Methods in Fluids 1995;38:1665–79. [5] Mukherjee YX, Mukherjee S. On boundary conditions in the element free Galerkin method. Computational Mechanics 1997;19:267–70. [6] Hegen D. Element-free Galerkin methods in combination with finite element approaches. Computer Methods in Applied Mechanics and Engineering 1996;135:143–66. [7] Zhu T, Atluri SN. Modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method. Computational Mechanics 1998;21:211–22. [8] Liew KM, Wu HY, Ng TY. Meshless method for modeling of human proximal femur: trestment of nonconvex boundaries and stress analysis. Computational Mechanics 2002;28:390–400. [9] Liew KM, Ng TY, Wu YC. Meshfree method for large deformation analysis—a reproducing kernel particle approach. Engineering Structures 2002;24:543–51. [10] Liew KM, Wu YC, Zou GP, Ng TY. Elasto-plasticity revisited: numerical analysis via reproducing kernel particle method and parametric quadratic programming. International Journal for Numerical Methods in Engineering 2002;55:669–83. [11] Liew KM, Huang YQ, Reddy JN. Moving least squares differential quadrature method and its application to the analysis of shear deformable plates. International Journal for Numerical Methods in Engineering 2003;56:2331–51.
[12] Ren J, Liew KM, Meguid SA. Modelling and simulation of the superelastic behaviour of shape memory alloys using the element-free Galerkin method. International Journal of Mechanical Sciences 2002;44:2393–413. [13] Gingold RA, Moraghan JJ. Smoothed particle hydrodynamics: theory and applications to non-spherical stars. Monthly Notices of the Astronomical Society 1977;181:375–89. [14] Onate E, Idelsohn S, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics: applications to convective transport and fluid flow. International Journal for Numerical Methods in Engineering 1990;39:3839–66. [15] Liszka TJ, Duarte CAM. Tworzydlo WW. hp-meshless cloud method. Computer Methods in Applied Mechanics and Engineering 1996;139:263–88. [16] Duarte CA, Oden JT. An h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 1996;139:237–62. [17] Aluru NR, Li G. Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation. International Journal for Numerical Methods in Engineering 2001;50:2373–410. [18] Ng TY, Li H, Cheng JQ, Lam KY. A new hybrid meshless-differential order reduction (hM-DOR) method with applications to shape control of smart structures via distributed sensors/ actuators. Engineering Structures 2003;25:141–54. [19] Liew KM, Huang YQ, Reddy JN. A hybrid moving least squares and differential quadrature (MLSDQ) meshfree method. International Journal of Computational Engineering Science 2002;3:1–12. [20] Liew KM, Lim HK, Tan MJ, He XQ. Analysis of laminated composite beams and plates with piezoelectric patches using the element-free Galerkin method. Computational Mechanics 2002;29:486–97. [21] Liew KM, Zhao X, Ng TY. The element-free kp-Ritz method for vibration of laminated rotating cylindrical panels. International Journal of Structural Stability and Dynamics 2002;2:523–58. [22] Ren J, Liew KM. Mesh-free method revisited: two new approaches for the treatment of essential boundary conditions. International Journal of Computational Engineering Science 2002;3:219–33. [23] Zhao X, Liew KM, Ng TY. Vibration analysis of laminated composite cylindrical panels via a meshfree approach. International Journal of Solids and Structures 2003;40:161–80. [24] Lam KY, Wang QX, Zong Z. A nonlinear fluid–structure interaction analysis of a near-bed submarine pipeline in a current. Journal of Fluids and Structures 2002;16:1177–91. [25] Reddy JN. An introduction to the finite element method. New York: McGraw-Hill; 1993.
A meshless Hermite-Cloud method for nonlinear fluid-structure analysis of near-bed submarine pipelines under current Hua Li a,, J.Q. Cheng a, T.Y. Ng a,b, Jun Chen a, K.Y. Lam a a
b
Institute of High Performance Computing, National University of Singapore, 1 Science Park Road, #01-01, The Capricorn, Singapore Science Park II, Singapore 117528, Singapore School of Mechanical & Production Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore Received 28 June 2003; received in revised form 1 December 2003; accepted 1 December 2003
Abstract Through higher-order Hermite construction of the interpolation functions, a novel true meshless numerical technique—the Hermite-Cloud method is developed for the analysis of nonlinear fluid-structure interaction of near-bed submarine pipelines under a current. Following the development of the discrete point collocation technique throughout the computational domains, the present method constructs the approximations of both the unknown functions and their first-order derivatives based on the classical reproducing kernel particle method, except that a fixed reproducing kernel approximation is employed instead. For partial differential equations with boundary conditions, certain auxiliary conditions are required due to the use of the Hermite theorem. After validation of the presently developed Hermite-Cloud method by several two-dimensional elasticity examples, the method is used for the failure analysis of a near-bed submarine pipeline under a horizontal current with consideration of a nonlinear fluid–structure interaction effect. The numerical result indicates that the distance from the pipelines to seabed and the current velocity have remarkable influences on the deformation behavior of the pipelines. It also shows that various critical current velocities corresponding to different failure patterns exist, and the presently simulated relation of the critical velocities agrees with the results in the previous publications. # 2003 Elsevier Ltd. All rights reserved. Keywords: Meshless method; Hermite-Cloud; RKPM; Point collocation; Nonlinear fluid–structure interaction; Submarine pipelines
1. Introduction The development of the finite element method (FEM) has formed the basis for the development of commercial computational mechanics based software for the modeling and simulation of various mechanical computer aided engineering (MCAE) problems. In this endeavor however, the requirement of mesh generation in the FEM results in certain disadvantageous features, such as iterative remeshing for the tracking of dynamic processes in large deformation problems and the requirement for large storage and memory due to the large number of element nodes involved in FEM discretization. To overcome these FEM associated deficiencies, various meshless numerical techniques [1] have
Corresponding author. Tel.: +65-64191561; fax: +65-64191380. E-mail address: [email protected] (H. Li).
0141-0296/$ - see front matter # 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2003.12.005
been recently developed, for example, the backgroundmesh type of approaches including the Galerkin-based techniques [2–12] and the direct domain discretization such as point collocation-based techniques [13–18]. Some interesting work should also be noted [19–23]. To explore more efficient collocation-based techniques, this paper develops a novel true meshless technique — the Hermite-Cloud method, which is based on the Hermite interpolation theorem and point collocation as well as reproducing kernels, and represents a significant development in reproducing kernel techniques. It employs the Hermite theorem for the construction of the interpolation functions, where the shape functions are constructed to correspond respectively to the unknown functions and the first-order derivatives of these unknown functions. It is formulated on the basis of the classical reproducing kernel particle method (RKPM) but fixed kernels are used here
532
H. Li et al. / Engineering Structures 26 (2004) 531–542
instead [17]. Furthermore, for a given set of the partial differential equations (PDEs) with the Dirichlet and/or Neumann boundary conditions, certain differential-type auxiliary conditions are derived based on the Hermite theorem to generate a complete set of equations for the governing partial differential boundary value (PDBV) problem. Following the scattering of a set of points in the computational domain including the edges, the point collocation technique is applied for the PDBV discretization. The approximate solutions of both the unknown functions and their first-order derivatives are thus expressed in terms of both the shape functions and unknown point values leading to a complete set of discrete algebraic equations being described. Finally, they are solved with respect to the unknown point values and the numerical solutions of the PDBV problem can be computed in a straightforward manner. To validate the developed Hermite-Cloud method, several two-dimensional (2-D) elasticity examples, including a higher-order patch subjected to a uniform unity-intensity load, cantilever beams subjected to a linearly varying axial and a shear end loads respectively, and a coupled thermo-elasticity problem, are analyzed numerically and the results are in good agreements with either exact results or existing numerical meshless solutions. For the application of the present Hermite-Cloud method in the submarine engineering, a nonlinear fluid–structure interaction of near-bed submarine pipelines under a current is studied here. The behavior of near-seabed submarine pipelines is greatly affected by the environment factors, such as ocean currents, earthquakes and underwater explosions. This paper investigates the influence of current on the submarine pipelines. When the distance between a supported submarine pipeline and the seabed is very large, the current around the pipeline is hardly affected by the seabed and the flow is symmetrical to the pipeline. In the vertical direction, there is no net force applied on the pipeline. However, when the distance comes to be small, the velocity of current between the pipeline and seabed becomes higher, compared with that above the pipeline. Then based on Bernoulli’s equation, the current pressure above the pipeline is higher than that between the pipeline and seabed. A downward net force caused by the pressure difference is thus produced and the pipeline is expected to deform under the external load. The deflection of pipelines enlarges with the increase of current velocity. When the velocity reaches a critical value, the pipeline will fail to work. In this paper, the two types of failure patterns are considered. One is caused by instability, in which the pipeline will rest on the seabed completely, and the other results from the limit of material properties, in which the stress/deflection caused by the external load reaches the allowable stress/deformation of the constitutive
materials. The variation of critical velocities with the distance from pipelines to seabed in various failure patterns is obtained. The distribution of the pipeline stress and the relation between the pipeline deflection and current velocity are simulated, respectively. The discussion on these computed results is also made in detail, which indicates a good agreement with the published work [19].
2. Development of meshless Hermite-Cloud method 2.1. Theoretical formulation As mentioned earlier, as one in the class of mathematical operators that can reproduce a function by integration transform over a computational domain, the reproducing kernel method for a 2-D real function f ðx;yÞ in a domain X is expressed as follows: ð f ðx;yÞ ¼ Uðx p;y qÞf ðp;qÞdpdq; ð1Þ X
where Uðx;yÞ is a real window function. According to the definition of the reproducing kernel method, the ideal window functions should be orthogonal, and its integration over X should be unity, so as to exactly reproduce f ðx;yÞ. The selection of window functions differentiates the various reproducing kernel methods. However, it is not often easy to select suitable functions as the window function must satisfy both the conditions to exactly reproduce the unknown function f ðxÞ. In the classical RKPM [3,4], an approximate solution is constructed by introducing a window function Uðx;yÞ as the product of a correction function Cðx;y;p;qÞ and a kernel function Kðx;yÞ, that is, Uðx p;y qÞ ¼ Cðx;y;p;qÞKðx p;y qÞ:
ð2Þ
Substituting Eq. (2) into Eq. (1), the approximation ~f ðx;yÞ of the unknown function f ðx;yÞ can be represented as follows ð ~f ðx;yÞ ¼ Cðx;y;p;qÞKðx p;y qÞf ðp;qÞdpdq: ð3Þ X
Further to that mentioned above, the selection of kernel functions also differentiates the various reproducing kernel methods. If the fixed kernel is used as the present kernel function Kðx;yÞ, Eq. (3) can be rewritten as ð ~f ðx;yÞ ¼ Cðx;y;p;qÞKðxk p;yk qÞf ðp;qÞdpdq ð4Þ X
where the point ðxk ;yk Þ is the center point of the fixed kernel corresponding to the kernel function Kðxk p;yk qÞ. Further, depending on different PDBV problems, the kernel function may be constructed by
H. Li et al. / Engineering Structures 26 (2004) 531–542
different forms of weighted window functions, such as the Gaussian functions, spline functions or radial basis functions. In the present Hermite-Cloud method, we employ a cubic spline function to construct the kernel function Kðxk p;yk qÞ ¼ W ððxk pÞ=DxÞW ððyk qÞ=DyÞ=ðDxDyÞ;
conditions, Eq. (9), is rewritten in discrete form as NT X Cðx;y;pn ;qn ÞKðxk pn ;yk qn Þbi ðpn ;qn ÞDSn
bi ðx;yÞ ¼
n¼1 NT X ¼ Bðpn ;qn ÞC ðx;yÞKðxk pn ;yk qn Þ n¼1
ð5Þ
where W ðzÞis a cubic-spline window function and is given as 8 0 j zj 2 > > > > < ð2 jzjÞ3 1 j zj 2 ð6Þ W ðzÞ ¼ 6 > > 2 j z j > > j zj 1 : z2 1 3 2 in which z ¼ ðxk pÞ=Dx or z ¼ ðyk qÞ=Dy. Dx and Dy denote the cloud size with the center point ðxk ;yk Þ. These are adjusted according to the point distribution and accuracy requirement, due to the consistency conditions, of the reproducing kernel approach. It is well-known mathematically that a continuous function can be represented as a sum of the linearly independent functions. As a result, the correction function Cðx;y;p;qÞ in Eq. (4) may be expressed as
bi ðpn ;qn ÞDSn
BT ðx;yÞ ¼ Aðxk ;yk ÞC ðx; yÞ
C ðx; yÞ ¼ A1 ðxk ; yk ÞBT ðx; yÞ;
ð12Þ
where Aðxk ;yk Þ is a symmetric constant matrix at the fixed-cloud center point ðxk ; yk Þ, Aij ðxk ;yk Þ ¼
NT X
bi ðpn ;qn ÞKðxk pn ;yk qn Þ
bj ðpn ;qn ÞDSn
ði; j ¼ 1; 2; . . . ; bÞ:
ð13Þ
By substituting Eqs. (5), (7) and (12) into Eq. (4), the approximate expression f~ðx;yÞ of the unknown real function f ðx;yÞ can be obtained as ð ~ f ðx;yÞ ¼ Bðp;qÞC ðx;yÞKðxk p;yk qÞf ðp;qÞdpdq ðX ¼ Bðp;qÞA1 ðxk ;yk ÞBT ðx;yÞ X
Kðxk p;yk qÞf ðp;qÞdpdq;
ð14Þ
which can further be discretized as ð8Þ
X
ði ¼ 1; 2; . . . ; bÞ:
ð11Þ
we obtain
ð7Þ
In the coefficient vector CT ðx;yÞ=fc1 ;c2 ; . . . ;cb g of Eq. (7), the coefficients ci ði ¼ 1; 2; . . . ; bÞ are unknown and can be determined through the following consistency conditions of the reproducing kernel techniques, ð bi ðx;yÞ ¼ Cðx;y;p;qÞKðxk p;yk qÞ bi ðp;qÞdpdq
ð10Þ
n¼1
where C ðx;yÞ is a bth-order column coefficient vector and Bðp;qÞ a bth-order row basis-function vector, in which bi ðp;qÞ ði ¼ 1; 2; . . . ; bÞ are linearly independent basis functions (b denotes the degree of the polynomials of the basis function). Usually their definitions depend on the PDBV problem to be solved. For example, for a 1-D PDE system, we may take BðpÞ ¼ f1; p; p2 ; . . . ; pb1 g. For a 2-D second-order PDE system, a quadratic basis function Bðp;qÞ is required Bðp;qÞ ¼ fb1 ðp;qÞ;b2 ðp;qÞ; . . . ;bb ðp;qÞg ¼ f1;p;q;p2 ;pq;q2 g ðb ¼ 6Þ:
ði ¼ 1;2; . . . ;bÞ;
where NT is the total number of scattered points covering the interior domain X and the surrounding edges. The subscript n represents the nth scattered point. DSn is defined as the cloud area of the nth point. It is obvious that Eq. (10) is a set of linear algebraic equations with respect to the coefficients ci ði ¼ 1; 2; . . . ; bÞ. Thus, if it is rewritten in the following matrix form with respect to the coefficient vector CT ðx;yÞ ¼ fc1 ; c2 ; . . . ; cb g,
Cðx;y;p;qÞ ¼ Bðp;qÞC ðx;yÞ ¼ fb1 ðp;qÞ;b2 ðp;qÞ; . . . ;bb ðp;qÞgfc1 ;c2 ; . . . ;cb gT ;
533
ð9Þ
If a set of points is scattered arbitrarily in the integral domain X and along its edges, the point collocation technique can be applied for PDBV discretization. By substituting Eq. (7) into Eq. (9), the consistency
f~ðx;yÞ ¼
NT X ðBðpn ;qn ÞA1 ðxk ;yk ÞBT ðx;yÞ n¼1
Kðxk pn ;yk qn ÞDSn Þfn ¼
NT X Nn ðx;yÞfn ;
ð15Þ
n¼1
where fn is defined as the point value of the nth point. Nn ðx;yÞ are the shape functions that consist of the basis functions and are simply polynomials in x and y. As such, any derivatives of the shape functions can be obtained easily by differentiation of the basis function vector Bðx;yÞ. Furthermore, it is evident that the present shape functions Nn ðx;yÞ satisfy the consistency
534
H. Li et al. / Engineering Structures 26 (2004) 531–542
conditions, Eqs. (9) or (10), for all independent basis functions bi ðx;yÞ ði¼ 1; 2; . . . ; bÞ. In particular, when b1 ðx;yÞ ¼ 1:0 (i¼ 1), b2 ðx;yÞ ¼ x (i¼ 2) and b3 ðx;yÞ ¼ y (i¼ 3) are taken for the discretized consistency conditions Eq. (10), we have 1:0 ¼
NT X
Nn ðx;yÞ; x ¼
n¼1
NT X Nn ðx;yÞxn ; n¼1
NT X Nn ðx;yÞyn : y¼
ð16Þ
n¼1
According to the Hermite theorem, the first-order derivatives of the unknown real function f ðx;yÞ with respect to the variables x and y are gx ðx;yÞ ¼
@f ðx;yÞ @f ðx;yÞ ; gy ðx;yÞ ¼ : @x @y
ð17Þ
n¼1
These result in additional unknown functions. If they are also discretized by imposing Eq. (15) in a similar manner, their approximate expressions can be written as ~ gx ðx;yÞ ¼
NS X
Mm ðx;yÞgxm
differential governing equations, consists of all three unknown point-value vectors ffn g1NT , fgxm g1NS and fgym g1NS , thus making it possible to directly compute the approximate solutions of both the unknown functions and their first-order derivatives. In addition, the Hermite-based interpolation approximation, Eq. (20), is also used for the construction of the following necessary auxiliary conditions. Due to the inclusion of additional unknown functions gx ðx;yÞ and gy ðx;yÞ, certain auxiliary conditions are required to generate a complete set of PDBV equations. Based on Eq. (17) for the defined gx ðx;yÞ and gy ðx;yÞ, the following auxiliary conditions are developed naturally by considering Eqs. (16), (18) and (20) ! NS NT NT X X X Nn;x ðx;yÞfn ðNn;x ðx;yÞxn Þ Mm ðx;yÞgxm
ð18Þ
m¼1
ð19Þ
where NS ð NT Þ is the total number of scattered points, and gxm and gym are the unknown point values for the mth point. Mm ðx;yÞ are the shape functions corresponding to the unknown first-order differential functions gx ðx;yÞ and gy ðx;yÞ, but their basis functions are now Bg ðx;yÞ 2 Rb1 . Based on the Hermite interpolation theorem, a true meshless approximation f~ðx;yÞ of the unknown real function f ðx;yÞ can now be finally constructed in the following form f~ðx;yÞ ¼
NT X Nn ðx;yÞfn n¼1
þ
NS X m¼1
þ
NS X m¼1
x
NT X
ðNn;x ðx;yÞyn Þ Mm ðx;yÞgym ¼ 0
ð21Þ
m¼1 NS NT X X m¼1
n¼1
!
ðNn;y ðx;yÞxn Þ Mm ðx;yÞgxm ¼ 0;
ð22Þ
n¼1
where the variable in the subscript after a comma indicates partial differentiation with respect to that variable. These auxiliary conditions are a necessary requirement in the implementation of the proposed Hermite-Cloud method. The above formulation thus defines the HermiteCloud method. In summary, it is based on the Hermite interpolation theorem, and consists of the approximation f~ðx;yÞ of the unknown function f ðx;yÞ, the approximations ~gx ðx;yÞ and ~gy ðx;yÞ of the first-order derivatives gx ðx;yÞ and gy ðx;yÞ, and further combines the necessary auxiliary conditions. 2.2. Numerical implementation
!
In general, engineering PDBV problems can be written as
Nn ðx;yÞxn Mm ðx;yÞgxm
n¼1
! NT X y Nn ðx;yÞyn Mm ðx;yÞgym :
!
n¼1
n¼1
m¼1
n¼1
! NS NT NT X X X Nn;y ðx;yÞfn ðNn;y ðx;yÞyn Þ Mm ðx;yÞgym
m¼1 NS X ~ gy ðx;yÞ ¼ Mm ðx;yÞgym ;
m¼1 NS NT X X
ð20Þ
n¼1
The presently constructed Hermite-based interpolation approximation has many computational advantages. Most notably, the computational accuracy at scattered discrete points in the domain is much refined not only for the approximate functions, but also for their firstorder derivatives. This can be easily understood from the subsequent formulation, Eqs. (29) and (32), where the final unknown vector of the set of the algebraic governing equations, derived by the discretized partial
Lf ðx;yÞ ¼ Pðx;yÞ PDEs in interior domain X
ð23Þ
f ðx;yÞ ¼ Qðx;yÞ Dirichlet boundary condition on CD ð24Þ @f ðx;yÞ ¼ Rðx;yÞ Neumann boundary condition on CN @n ð25Þ where L is a differential operator and f ðx;yÞ an unknown real function. By the point collocation technique and taking f~ðx;yÞ as the approximation of f ðx;yÞ, the problem is discretized and expressed
H. Li et al. / Engineering Structures 26 (2004) 531–542
535
approximately by Lf~ðxi ;yi Þ ¼ Pðxi ;yi Þ i ¼ 1;2; . . . ;NX
ð26Þ
f~ðxi ;yi Þ ¼ Qðxi ;yi Þ i ¼ 1;2; . . . ;ND
ð27Þ
@ f~ðxi ;yi Þ ¼ Rðxi ;yi Þ i ¼ 1;2; . . . ;NN ; @n
ð28Þ
where NX, ND and NN are the numbers of scattered points in the interior domain, and along the Dirichlet and Neumann edges, respectively, and the total number of scattered points is thus NT ¼ ðNX þND þNN Þ. Substituting the approximations, Eqs. (18) and (20), into Eqs. (26) and (28), and further combining the auxiliary conditions of Eqs. (21) and (22), followed by rearrangement of the resulting equations, a set of discrete algebraic governing equations with respect to the unknown point values fi , gxi and gyi , is obtained and can be written in matrix form as ½Hij ðNT þ2NS ÞðNT þ2NS Þ fFi gðNT þ2NS Þ1 ¼ fdi gðNT þ2NS Þ1 ; ð29Þ where fdi g and fFi g are ðNT þ 2NS Þ-order column vectors, with fFi gðNT þ2NS Þ1 ¼ fffi g1NT ; fgxi g1NS ; fgyi g1NS gT
ð30Þ
fdi gðNT þ2NS Þ1 ¼ ffPðxi ;yi Þg1NX ; fQðxi ;yi Þg1ND ; fRðxi ;yi Þg1NN ; f0g12NS gT
ð31Þ
and [Hij] is a ðNT þ2NS Þ ðNT þ2NS Þ coefficient square matrix 2
6 LNj ðxi ;yi Þ NX NT 6 6 6 Nj ðxi ;yi Þ ND NT 6 6 ½0 NN NT ½Hij ¼ 6 6 6 6 Nj;x ðxi ;yi Þ NS NT 6 6 4 Nj;y ðxi ;yi Þ NS NT
Fig. 1. Geometry and point distribution for the higher-order patch subjected to a uniform unidirectional stress of unity magnitude.
3. Validation study of Hermite-Cloud method To validate the numerical stability and accuracy of the present Hermite-Cloud method, numerical comparisons are made in this section for several 2-D planestress elasticity problems. These include a higher-order patch subjected to a uniform unity-intensity unidirectional stress, cantilever beams subjected to linearly varying axial loads or shear end loads, and a 2-D thermo-elasticity analysis. In order to measure the numerical accuracy of the present Hermite-Cloud method, a refined global error n is defined [5] for the numerical comparison: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u NT X 1 u t 1 ð33Þ ðf~i fi Þ2 : n¼ jfmax j NT i¼1 3.1. Higher-order patch For an isotropic higher-order patch subjected to uniform unidirectional stress of unity magnitude, as
NT P L xi Nn ðxi ;yi Þxn Mj ðxi ;yi Þ n¼1
n¼1
n¼1
NX NS
½0ND NS Mj ðxi ;yi Þ NN NS N PT Nn;x ðxi ;yi Þxn Mj ðxi ;yi Þ NS NS n¼1 NT P Nn;x ðxi ;yi Þxn Mj ðxi ;yi Þ NS NS
3
NT P L yi Nn ðxi ;yi Þyn Mj ðxi ;yi Þ NX NS
½0ND NS Mj ðxi ;yi Þ NN NS N PT Nn;x ðxi ;yi Þyn Mj ðxi ;yi Þ NS NS n¼1 NT P Nn;x ðxi ;yi Þyn Mj ðxi ;yi Þ NS NS
7 7 7 7 7 7 7: 7 7 7 7 7 5
n¼1
ð32Þ
Solving numerically the complete set of linear algebraic equations, Eq. (29), (NT þ NS ) point values {Fi}, consisting of the NT point values {Fi} and 2NS point values {gxi} and {gyi}, are obtained and the approximate solution f~ðx;yÞ and its first-order derivatives ~ gx ðx;yÞ and ~ gy ðx;yÞ of the PDBV problem can be computed through Eqs. (20), (18) and (19).
shown in Fig. 1, the plane-stress equilibrium governing equations is given in terms of the displacements u and v as follows: 2 E @ u 1 l @2u 1 þ l @2v þ þ ð34Þ þ Xx ¼ 0 1 l2 @x2 2 @y2 2 @x@y 2 E @ v 1 l @2v 1 þ l @2u þ þ ð35Þ þ Xy ¼ 0; 1 l2 @y2 2 @x2 2 @x@y
536
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 2. Numerical comparison for the higher-order patch subjected to a uniform unidirectional stress of unity magnitude: (a) displacement u, (b) displacement v.
where E is the elastic modulus and l the Poisson’s ratio. Xx and Xy denote the body forces in the x and y directions, respectively. If the material constants E ¼ 1 and l ¼ 0:25 are taken, the exact solution of this problem is uðx;yÞ ¼ x;
vðx;yÞ ¼ 0:25y:
ð36Þ
After implementation of the presently developed Hermite-Cloud method, the displacements u and v are obtained numerically and their profiles are plotted in Fig. 2(a,b). It is observed that, for a point distribution 53, comparisons with the exaction solution reveal global errors of n ¼ 3:19 107 for the displacement u, and n ¼ 3:32 107 for the displacement v. It is clearly evident that the present Hermite-Cloud method can achieve almost exact results at the scattered discrete points even with a coarse 53 point distribution for the present patch test problem. This example clearly demonstrates the efficiency and ease of implementation of the present Hermite-Cloud method are much enhanced over that of existing generic meshless techniques, such as the EFG method, since the present Hermite-Cloud method does not require a background mesh and has better consistency characteristics in the construction of the shape functions.
The exact solution for the case shown in Fig. 3 is given by uðx;yÞ ¼ xy; vðx;yÞ ¼
ð4x2 þ y2 Þ : 8
ð37Þ
Applying the Hermite-Cloud method for this elastic analysis, the numerical displacements u and v are computed. The displacement profiles are shown in Fig. 4(a,b). It is again observed that, for the same coarse 53 point distribution, comparison with the exact solution reveals global errors of n¼ 1:61 106 for the displacement u and n¼ 3:86 106 for the displacement v. In addition, the exact solution for the cantilever beam subjected to a shear end load, as shown in Fig. 5, is given by uðx;yÞ ¼
Pðy 0:5DÞ ð6L 3xÞx þ ð2 þ lÞðy2 2DyÞ 6EI ð38Þ
P 3lðy2 2Dy þ 0:5D2 ÞðL xÞ 6EI þð1 þ 1:25lÞD2 x þ ð3L xÞx2 ;
vðx;yÞ ¼
ð39Þ
3.2. Cantilever beam subjected to various loadings Here, a 2-D plane-stress cantilever beam is considered with the same equilibrium governing equations as Eqs. (34) and (35). It is subjected to a linearly varying axial end load as shown in Fig. 3, and subjected to a shear end load as shown in Fig. 5, respectively.
Fig. 3. Geometry and point distribution for a cantilever beam subjected to a linearly varying axial end load.
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 4.
537
Numerical comparison for the cantilever beam subjected to a linearly varying axial end load: (a) displacement u, (b) displacement v.
where I ¼ D3 =12 and D and L are respectively the height and length of the cantilever beam. The displacement fields are solved numerically using the Hermite-Cloud method, and the computed displacement profiles are shown in Fig. 6(a,b). Further, in order to examine the convergence of the Hermite-Cloud method, Table 1 compares the convergence characteristics between the present Hermite-Cloud method and
the meshless hM-DOR method [18]. It is observed from the table that global errors of both the displacement u and first-order derivative u,x in the present HermiteCloud method decrease monotonically with increasing number of scattered points, and the global-errors are lower than those of the hM-DOR method [18]. These numerical results thus show that the computational accuracy of the Hermite-Cloud method at scattered discrete points in the domain is much refined not only for the approximate solutions, but also for the first-order derivatives of these solutions. 3.3. Coupled thermo-elasticity
Fig. 5. Geometry and point distribution for a cantilever beam subjected to a shear end load.
Fig. 6.
As shown in Fig. 7, a plane-stress analysis of coupled thermo-elasticity is carried out by the developed Hermite-Cloud method. The partial differential governing equations that couple the heat conduction equation with the generic thermo-stress elastic equations, given
Numerical comparison for the cantilever beam subjected to a shear end load: (a) displacement u, (b) displacement v.
538
H. Li et al. / Engineering Structures 26 (2004) 531–542
Table 1 Numerical comparisons between the present Hermite-Cloud method and the hM-DOR [18] method for a cantilever beam subjected to a shear end load (E ¼ 3:0 107 , D ¼ 1, L ¼ 8 and l ¼ 0:25) Points distribution(NxNy)
511 1111 2111 6511 8111
Global error (n) for u
Global error (n) for u,x
hM-DOR [18] (%)
Hermite-Cloud (%)
hM-DOR [18] (%)
Hermite-Cloud (%)
35.9 17.9 6.16 4.28 4.01
35.9 17.7 5.92 1.95 0.384
30.9 15.0 5.38 3.92 3.90
30.9 14.9 5.19 2.13 1.92
in terms of the displacements u and v and temperature T, are @2T @2T þ 2 ¼h @x2 @y
ð40Þ
@2u 1 l @2u 1 þ l @2v @ðT T0 Þ ð1 þ lÞa ¼0 þ þ 2 2 @x 2 @y 2 @x@y @x ð41Þ @2v 1 l @2v 1 þ l @2u @ðT T0 Þ ð1 þ lÞa ¼ 0; þ þ @y2 2 @x2 2 @x@y @y ð42Þ where h is the heat source, T0 the initial reference temperature and a the coefficient of thermal expansion. The considered boundary conditions are uðx;yÞ ¼ 0; rxy ðx;yÞ ¼ 0; Tx ðx;yÞ ¼ 0; at x ¼ 0 and 1;
ð43Þ
vðx;yÞ ¼ 0; rxy ðx;yÞ ¼ 0; Tðx;yÞ ¼ 0; at y ¼ 0
ð44Þ
rxy ðx;yÞ ¼ 0; ryy ðx;yÞ ¼ 0; Tðx;yÞ ¼ 0; at y ¼ 1:
ð45Þ
7
If we take T0 ¼ 0, h ¼ 20, E ¼ 3:0 10 , a ¼ 0:001 and l ¼ 0:25, the exact solutions of this coupled thermo-elasticity problem are Tðx;yÞ ¼ 10y2 ; uðx;yÞ ¼ 0; vðx;yÞ ¼ 10að1 þ lÞy3 =3 ¼ 0:0125y3 =3:
Fig. 7. Geometry of the 2-D thermoelasticity problem.
Again by using the Hermite-Cloud method, the variable fields Tðx;yÞ, uðx;yÞ and vðx;yÞ and their firstorder derivatives are obtained numerically, with a regular 55 point distribution. The numerical results show that the standard error of the displacement u is of the order of 1011 and for the temperature Tðx;yÞ, the global error is less than 1.0107. A probable reason for the excellent results is that the exact solutions of u and T can be found in the basis functions. However, the same cannot be said for the numerical results of the displacement v since its exact cubic solution is not found in any of the basis functions. Compared with exact solution, Fig. 8 presents the accuracy variation of the displacement v with the increase of the regularly scattered points, for NT ðtotal point numberÞ ¼ 25 ð5 5Þ; 121 ð11 11Þ and 441 ð21 21Þ. It is evident that the computed accuracy of displacement v is significantly refined with the increase of the number of scattered points. Fig. 9 compares the convergence characteristics of the present displacement vðx;yÞ and the first-derivatives v;y ðx;yÞ, with that obtained by the meshless Finite Cloud method [17]. Present convergence rates are 2.12 for vðx;yÞ and 1.89 for v;y ðx;yÞwhile that for the Finite Cloud method are 2.2 for vðx;yÞ and 1.95 for v;y ðx;yÞ. It is therefore clear that the present Hermite-Could method attains better
ð46Þ
Fig. 8. Variation of the numerical displacement v with the point distribution density for the thermo-elasticity case.
H. Li et al. / Engineering Structures 26 (2004) 531–542
539
ditions for the above problem are given as
Fig. 9. Convergence comparison between the present HermiteCloud method and Finite- Cloud method [17] for the thermo-elasticity case (n, global error; h, point distance).
computational accuracy for both the field-variable solution and the first-order derivative than the Finite Cloud method.
; on CW ; wðx0 Þ ¼ w
ð49Þ
hðx0 Þ ¼ h; on Ch @h ; on CM Mðx0 Þ ¼ EI ¼M @x x¼x0 @w ; on CV V ðx0 Þ ¼ GAks h þ ¼V @x x¼x0
ð50Þ ð51Þ ð52Þ
in which Cw , Ch , CM and CV are the boundaries where w, h, M and V satisfy, respectively. In Eq. (47), f is the fluid force caused by the current pressure difference. Although it has an analytical solution, it converges slowly. In order to speed up the computation, Lam et al. used the boundary element method (BEM) to construct the approximation of the exact solution. They are in a very good agreement with each other [24]. The expression of fluid force f(x) is given as f ðxÞ ¼
1 qAU02 cðdÞ; 2
ð53Þ
where qis the mass density, U0 the current velocity, c and d are dimensionless coefficient and have the relation as 4. Numerical study of a near-bed submarine pipeline under current When a submarine pipeline is placed near the seabed, the horizontal current is not symmetrical to the pipeline due to the existence of seabed. As a result, a downward external force is produced and applied on the pipeline, resulting in the deformation of pipeline. This problem can be simplified as a beam with different boundary support conditions. In the previously published work [24], Lam et al. studied the steady behavior of the pipeline under a steady nonlinear external force, in which the pipeline was simplified as a Bernoulli– Euler beam by the assumption that there is no transverse shear strain. However it is not applicable in the case that the ratio of length to diameter of the pipeline is not very large, where the effect of shear deformation should be included. In order to break through the limit, the Timoshenko beam theory [25] is introduced here to include the transverse shear deformation. Based on the Timoshenko beam theory, the general governing equations of the pipeline are written as @ @w GAks þh þf ¼0 ð47Þ @x @x @ @h @w EI þ h ¼ 0; ð48Þ GAks @x @x @x where w is the deflection of the pipeline, h the rotation, G the shear modulus, A the cross section area, ks the shear correction coefficient, E the elasticity modulus, I the moment of inertia. The general boundary con-
cðdÞ ¼
2:23d 2 þ 12:54d þ 0:02 0:77d 3 þ 0:44d 2 þ 0:02d
ð54Þ
in which d is defined as d ¼ ðD0 wðxÞ Rs Þ=ð2Rs Þ, and D0 is the distance between the central line of pipeline at initial status and the seabed, Rs the pipe outer radius. It is noted from the Eqs. (53) and (54) that the fluid force f(x) is a nonlinear term related to the deflection w(x). The governing Eqs. (47) and (48) for the pipeline deformation thus are the coupled nonlinear functions. It is also observed that, when the material and dimension of pipelines are given, the fluid force f(x) is a function of both the current velocity U0 and the gap D0. Namely these two parameters have significant influences on the deformation behavior of the near-bed submarine pipelines. In order to demonstrate the influences in detail, a practical example is designed here. A circular steel pipeline is placed near the seabed under a steady current, as shown in Fig. 10. Several physical and material parameters are given, the pipeline length L¼ 10 m, Rs ¼ 0:5 m, the thickness of steel pipeline ts ¼ 0:02 m, D0 ¼ 0:7 m, Es ¼ 2:11 1011 N=m2 , Poisson ratio m ¼ 0:3, qs ¼ 7800 kg=m3 , steel yield stress Ys ¼ 2:5 108 N=m2 , and allowable deflection/span ratio Yw ¼ 0:004. It is noted that the ratio L/Rs is not very large. As such, the pipeline can be simplified as a thick beam fixed at two ends. By the Timoshenko beam theory and HermiteCloud meshless method, numerical analysis is carried out and discussions are also presented in succession.
540
H. Li et al. / Engineering Structures 26 (2004) 531–542
Fig. 10. Schematic diagram of a submarine pipeline’s deformation under a current.
The variation of mid-point deflection wL=2 with the current velocity U0 is illustrated in Fig. 11. It is observed that, with the increase of current velocity, the deflection enlarges until a critical value Ucb. When the velocity is smaller than the critical value, the deflection converges with the iterations and the pipeline reaches the stable equilibrium. However, if the velocity is over the critical one, the deflection does not converge, which is defined as the instability, i.e. the pipeline has fallen into the seabed. It is the first failure pattern of pipelines—the instability failure, as mentioned above. The influence of the distance D0 between the pipeline and seabed under the critical current velocities caused by instability Ucb is revealed in Fig. 12. It is shown that Ucb is a monotonically increasing function with D0. As indicated in Eqs. (53) and (54), if other variables are fixed, the fluid force f(x) decreases with the enlargement of D0 and increases with the increasing U0. Therefore, in order to obtain the same critical fluid force causing the pipeline to lose the stability, when the gap D0 becomes large, it is reasonable to require a higher critical velocity Ucb. It is also expected that the critical velocity is infinite when the gap becomes infi-
nite, since the effect of seabed on the current is then zero and the flow symmetry is not changed at all. As well known, sometimes the failure due to material properties happens before the instability failure. Thus it is necessary to discuss material failure. The distribution of stress along the pipeline is plotted in Fig. 13,
Fig. 11. Variation of the deflection of mid-point of the pipeline with respect to the current velocity U0 (when D0 ¼ 0:7 m).
Fig. 13. Distribution of the stress along the pipeline (when U 0 ¼ 10 m=s and D0 ¼ 0:7 m).
Fig. 12. Effect of the gap D0 on the critical velocities Ucb of the instability failure.
H. Li et al. / Engineering Structures 26 (2004) 531–542
541
Fig. 14. Critical velocities of (a) strength failure, (b) deflection failure (when D0 ¼ 0:7 m).
in which the curve is found to be similar to that of a two-end fixed beam under uniform load, but actually they are different. For the beam under uniform load, the maximal absolute value of stresses (moments) along the beam occurs at the mid-point, yet on the contrary, the present maximal absolute value of stresses (moments) appears at the end points, as shown in Fig. 13. This phenomenon results from the different external load distribution. In this problem, the fluid force f(x) is a decreasing function of pipeline deflection w(x), and it is clear that the deflection of two ends is smaller than that of mid-point. This implies that the fluid force at two ends is larger than that of mid-point. Obviously it is different from the case of uniform load distribution. Fig. 14 presents the critical velocities caused by the two types of material failures—the strength failure and the deflection failure. As mentioned above, the endpoint stress and mid-point deflection firstly reach the
respective critical values along the whole pipeline. Thus their corresponding critical current velocities are defined as the critical values of the whole system. In the Fig. 14(a,b), both the end-point stress and midpoint deflection enlarge with the increase of current velocity U0, and when they reach the yield stress and allowable deflection, the corresponding critical velocity Ucs and Ucw are obtained. Similar to the instability failure, the distance between the pipeline and seabed D0 has an influence on the critical velocities of material failure. In order to compare these two failure patterns, the distributions of respective critical velocities (Ucb due to instability failure, Ucs due to strength failure and Ucw due to deflection failure) with respect to the gap D0 are synthetically plotted in Fig. 15. It is found that there are bifurcations at point A and B. In other words, beyond point A, the pipeline will not break down due to the stability loss, because the strength failure always appears before the instability failure. The results agree with Lam’s work [24], in which the pipeline is simplified as a Bernoulli–Euler beam. 5. Conclusions
Fig. 15. Comparison of distributions of respective critical velocity with respect to the gap D0 in various failure patterns.
This paper presents the development of a novel true meshless numerical technique, termed the HermiteCloud method. It employs the Hermite interpolation theorem to construct the interpolation functions and uses the point collocation technique for discretization of partial differential equations. The necessary auxiliary conditions are derived for the complete set of partial differential equations with mixed Dirichlet and Neumann boundary conditions. Various 2-D numerical studies have been conducted to validate the convergence characteristics and numerical accuracy of the Hermite-Cloud method and it is concluded that the developed Hermite-Cloud method is a very efficient
542
H. Li et al. / Engineering Structures 26 (2004) 531–542
numerical technique for solving general PDVPs. By the Hermite-Cloud method, the analysis of the nonlinear fluid-structure interaction of near-bed submarine pipelines under a current is completed. The numerical results illustrate the relation between the critical current velocities and the distance from pipelines to the seabed in different pipeline failure pattern. They are in good agreements with previous work.
References [1] Liu GR. Mesh free methods: moving beyond the finite element method. CRC Press LLC; 2003. [2] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. International Journal for Numerical Methods in Engineering 1994;37:229–56. [3] Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. International Journal for Numerical Methods in Engineering 1995;20:1081–106. [4] Liu WK, Jun S, Li S, Adde J, Belytschko T. Reproducing kernel particle methods for structural dynamics. International Journal for Numerical Methods in Fluids 1995;38:1665–79. [5] Mukherjee YX, Mukherjee S. On boundary conditions in the element free Galerkin method. Computational Mechanics 1997;19:267–70. [6] Hegen D. Element-free Galerkin methods in combination with finite element approaches. Computer Methods in Applied Mechanics and Engineering 1996;135:143–66. [7] Zhu T, Atluri SN. Modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method. Computational Mechanics 1998;21:211–22. [8] Liew KM, Wu HY, Ng TY. Meshless method for modeling of human proximal femur: trestment of nonconvex boundaries and stress analysis. Computational Mechanics 2002;28:390–400. [9] Liew KM, Ng TY, Wu YC. Meshfree method for large deformation analysis—a reproducing kernel particle approach. Engineering Structures 2002;24:543–51. [10] Liew KM, Wu YC, Zou GP, Ng TY. Elasto-plasticity revisited: numerical analysis via reproducing kernel particle method and parametric quadratic programming. International Journal for Numerical Methods in Engineering 2002;55:669–83. [11] Liew KM, Huang YQ, Reddy JN. Moving least squares differential quadrature method and its application to the analysis of shear deformable plates. International Journal for Numerical Methods in Engineering 2003;56:2331–51.
[12] Ren J, Liew KM, Meguid SA. Modelling and simulation of the superelastic behaviour of shape memory alloys using the element-free Galerkin method. International Journal of Mechanical Sciences 2002;44:2393–413. [13] Gingold RA, Moraghan JJ. Smoothed particle hydrodynamics: theory and applications to non-spherical stars. Monthly Notices of the Astronomical Society 1977;181:375–89. [14] Onate E, Idelsohn S, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics: applications to convective transport and fluid flow. International Journal for Numerical Methods in Engineering 1990;39:3839–66. [15] Liszka TJ, Duarte CAM. Tworzydlo WW. hp-meshless cloud method. Computer Methods in Applied Mechanics and Engineering 1996;139:263–88. [16] Duarte CA, Oden JT. An h-p adaptive method using clouds. Computer Methods in Applied Mechanics and Engineering 1996;139:237–62. [17] Aluru NR, Li G. Finite cloud method: a true meshless technique based on a fixed reproducing kernel approximation. International Journal for Numerical Methods in Engineering 2001;50:2373–410. [18] Ng TY, Li H, Cheng JQ, Lam KY. A new hybrid meshless-differential order reduction (hM-DOR) method with applications to shape control of smart structures via distributed sensors/ actuators. Engineering Structures 2003;25:141–54. [19] Liew KM, Huang YQ, Reddy JN. A hybrid moving least squares and differential quadrature (MLSDQ) meshfree method. International Journal of Computational Engineering Science 2002;3:1–12. [20] Liew KM, Lim HK, Tan MJ, He XQ. Analysis of laminated composite beams and plates with piezoelectric patches using the element-free Galerkin method. Computational Mechanics 2002;29:486–97. [21] Liew KM, Zhao X, Ng TY. The element-free kp-Ritz method for vibration of laminated rotating cylindrical panels. International Journal of Structural Stability and Dynamics 2002;2:523–58. [22] Ren J, Liew KM. Mesh-free method revisited: two new approaches for the treatment of essential boundary conditions. International Journal of Computational Engineering Science 2002;3:219–33. [23] Zhao X, Liew KM, Ng TY. Vibration analysis of laminated composite cylindrical panels via a meshfree approach. International Journal of Solids and Structures 2003;40:161–80. [24] Lam KY, Wang QX, Zong Z. A nonlinear fluid–structure interaction analysis of a near-bed submarine pipeline in a current. Journal of Fluids and Structures 2002;16:1177–91. [25] Reddy JN. An introduction to the finite element method. New York: McGraw-Hill; 1993.