Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390 www.elsevier.com/locate/cma
A finite volume method for incompressible linear elasticity I. Bijelonja a, I. Demirdzˇic´ a
b,* ,
S. Muzaferija
b
Masˇinski fakultet Sarajevo, Vilsonovo sˇetalisˇte 9, 71000 Sarajevo, Bosnia and Herzegovina b CD-adapco, Du¨rrenhofstrasse 4, D-90402 Nu¨rnberg, Germany
Received 6 June 2005; received in revised form 29 November 2005; accepted 9 January 2006
Abstract This paper describes development of a finite volume based method for modelling of incompressible linear elastic body problems. The method is based on the solution of the integral form of conservation equations governing momentum balance and the introduction of pressure as an additional variable. The incompressibility is enforced by employing the volumetric constraint equation. A collocated variable arrangement is used and the spatial domain is discretised using finite volumes of arbitrary polyhedral shapes. A segregated approach is employed to solve resulting set of coupled linear algebraic equations, embedding a SIMPLE based algorithm for displacement–pressure coupling. Numerical experiments showed that the developed method appears to be locking free. 2006 Elsevier B.V. All rights reserved. Keywords: Finite volume method; Incompressible body; Elasticity
1. Introduction Many problems of physical importance involve analysis of deformation of incompressible or almost incompressible solids. For example, some elastomers, rubber-like materials and materials in inelastic conditions, may exhibit almost incompressible response. Sometime the compressibility effects may be so small that they could be neglected, and the material would be idealized as incompressible. It is well known that numerical analysis of solids in incompressible limit could lead to difficulties. For example, fully integrated displacement based lower-order finite elements suffer from volumetric locking, which usually accompanies pressure oscillation in incompressible limit. The volumetric locking phenomenon and pressure oscillation have been studied extensively and a number of approaches have been proposed to overcome the problem. Methods that have been proposed under the finite element framework include the mixed formulation [1–3], the selective reduced integration method [4], the hourglass control method on under-integrated elements [5], the B-bar typed formulation [6,7], and the pressure projection method [8]. Several post-processing pressure smoothing techniques [9] have been also developed to produce a smooth pressure distribution. These techniques seem to be able to smooth out non-physical pressure oscillations, but require an additional computation effort. At present time, contrary to what was originally claimed, it is also clear that meshless methods could experience volumetric locking [10]. The interest for finite volume (FV) application to the structural analysis problems involving incompressible materials has grown recently. Henry and Collins [11] introduced the linear elastic incompressible material model undergoing small axisymmetric deformation using finite volume modelling. To account for volumetric invariance, a condition between the radial and axial displacements and the hydrostatic pressure as the additional dependent variable were used. The rectangu*
Corresponding author. Tel.: +49 911 946 4362; fax: +49 911 946 4366. E-mail address:
[email protected] (I. Demirdzˇic´).
0045-7825/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.01.005
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
6379
lar cells and SIMPLEC [12] algorithm were employed in discretisation. Later, the model was incorporated into a commercially available FLOW3D code, in order to model the blood flow in artery [13]. In these works no locking phenomenon is registered. A small strain incompressible material model was introduced to a geometrically unrestricted finite volume model in the form of the mixed formulation [14,15]. The auxiliary pressure equation was used to impose incompressibility constraints. From the results of several benchmark solutions, the finite volume method appeared to offer a number of advantages over equivalent finite element models as the finite volume solution was conservative and incompressibility was satisfied exactly for each discretisation element of the solution domain, rather than in average sense. Although certain restrictions on mesh configuration had to be imposed to avoid locking, these restrictions were less severe than those of the equivalent finite element meshes. Numerical calculation with meshes consisting of triangular cells showed excellent agreement with analytical results. Meshes consisting of quadrilateral FV cells displayed too stiff behavior, indicating a locking phenomenon. Employing the fact that the equation of Stokes flow are identical in form to the equations of isotropic incompressible linear elasticity (only the physical interpretation of the variables is different), the model developed for the analysis of incompressible fluids [16] is adopted in this paper to describe, in a displacement–pressure based formulation, deformation of an incompressible elastic solid. The method is based on the solution of the integral form of momentum balance equation. Incompressibility constraint is enforced using an integral form of the volumetric constraint equation and the hydrostatic pressure as an additional dependent variable. The method employs numerical meshes consisting of contiguous control volumes of completely arbitrary topology. A collocated variable arrangement is used, and a segregated approach is employed to solve resulting set of coupled linear algebraic equations, embedding a SIMPLE [12] based algorithm for calculation of hydrostatic pressure featuring in the constitutive equation. The method is applied to several test cases for which an analytical and finite element solution exists. Comparisons of the numerical and analytical results show a very good performance of the method. Numerical experiments show that developed discretisation scheme is free of locking for finite volume meshes of arbitrary topology. In the next section the governing equations together with constitutive relations are given. This is followed by a brief description of the finite volume discretisation procedure and the solution algorithm. Finally, the method’s capabilities are demonstrated by applying it to a number of test cases. 2. Mathematical formulation and problem description The behavior of a continuum is governed by the following momentum transport equation: Z Z Z o qv dV ¼ r ds þ qf b dV ; ot V S V
ð1Þ
which is valid for an arbitrary part of continuum of the volume V bounded by the outward pointing surface vector s. In Eq. (1), q is the mass density, v is the velocity, r is the Cauchy stress tensor, and fb is the resulting body force. It is assumed here that angular momentum balance equation is satisfied exactly due the shear stresses conjugate principle. The constitutive equations in an isotropic linear elastic solid may be written as r ¼ 2l þ k div uI;
ð2Þ
where 1 ¼ ½grad u þ ðgrad uÞT 2
ð3Þ
is the linear strain tensor, u is the displacement vector, I is the identity tensor, k and l are the Lame´ constants, related to the is Young’s modulus E and the Poisson’s ratio m by the following relations: k¼
Em ; ð1 þ mÞð1 2mÞ
l¼
E . 2ð1 þ mÞ
ð4Þ
In the incompressible limit m ! 1/2, the Lame´ constant k becomes unbounded, the pressure cannot be determined from the strain field, and an alternative formulation of the constitutive law is therefore necessary. In this formulation the constitutive equation is written as r ¼ 2l pI;
ð5Þ
where the pressure p must be determined as a part of the solution to the boundary-value problem and thus represents an additional unknown. The additional equation that needs to be introduced is the kinematic condition of incompressibility: div u ¼ 0.
ð6Þ
6380
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
2.1. Resulting set of equations Introducing constitutive relation (5) for incompressible solid into momentum balance equation (1), the linear elastic solid body problem reduces to the solution of the following equation: Z Z Z Z o ou T q dV ¼ lðgrad uÞ ds þ ½l grad u pI ds þ qf b dV . ð7Þ ot V ot S S V Integrating the kinematic condition of incompressibility (6) and using the Gauss divergence theorem the following equation is obtained: Z Z div u dV ¼ u ds ¼ 0. ð8Þ V
S
Eqs. (7) and (8), make a closed set of four equations with four unknown functions (three displacement vector components ui and pressure p) of spatial coordinates and time. 2.2. Initial and boundary condition To complete the mathematical model, initial and boundary conditions have to be specified. As initial conditions, the pressure, displacements, and in transient cases velocities, have to be specified at all points of the solution domain. Boundary conditions have to be specified at all times at all solution domain boundaries. They can be either of Dirichlet type (displacement boundary condition) or of Neumann type (traction boundary condition). 3. Numerical method In this section the equation governing momentum balance (7) and the kinematic condition of incompressibility (8) are discretised by employing the finite volume method described in detail in [16,17]. Note that system of Eqs. (7) and (8) can be viewed as the equations of Stokes flow of incompressible fluid, with different physical interpretation of variables. From numerical point of view it is important to note that diffusivity coefficient l, which presents viscosity in fluid flow problems, is much smaller than shear modulus which play a role of diffusivity coefficient in solid body analysis. Having all this in mind, the numerical procedure developed in [16] for analysis of incompressible fluid flow is adopted here for discretisation of incompressible solid body equations (7) and (8). 3.1. Space and time discretisation In order to obtain discrete counterpart of Eqs. (7) and (8), the time interval of interest is divided into an arbitrary number of time increments dt, and the space is discretised by a number of contiguous, non-overlapping control volumes (CV), with computational points at their centers (Fig. 1). The boundary nodes, needed for the specification of boundary conditions, reside at the center of boundary cell faces. 3.2. Discretisation of momentum equation After the space discretisation, the momentum equation (7) is written for each displacement component ui (i = 1, 2, 3) and for each control volume as follows:
Fig. 1. A control volume of an arbitrary shape.
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
o ot
Z
f X oui dV ¼ q ot j¼1
n
V
Rate of change
Z
l grad ui ds þ
Sj
nf Z X j¼1
½lðgrad uÞ ii pii ds þ
Sj
Diffusion
6381
Z qfbi dV V
ð9Þ
Source terms
where ii are the Cartesian base vectors and nf is the number of faces enclosing the cell P0. In order to evaluate integrals in the above equation, a distribution of dependent variables and physical properties of material in space and time have to be approximated. 3.2.1. Spatial distribution The following linear spatial distribution is used: wðrÞ ¼ wP 0 þ ðgrad wÞP 0 ðr rP 0 Þ;
ð10Þ
where w stands for dependent variables ui or p, or a physical material property, rP 0 is the position vector of control volume center P0 (Fig. 1). The unknown vector ðgrad wÞP 0 is calculated by applying a least square fit to Eq. (10): ðgrad wÞP 0 ¼ D1 c;
ð11Þ
where D¼
nf X T ðrP j rP 0 Þ ðrP j rP 0 Þ; j¼1
c¼
nf X T ðrP j rP 0 Þ ðwP j wP 0 Þ.
ð12Þ
j¼1
Eq. (10) is used to calculate values of w at the cell-faces j necessary for evaluation of the surfaces integrals in (9), leading to a second-order symmetric formula: 1 1 wj ¼ ðwP 0 þ wP j Þ þ ½ðgrad wÞP 0 ðrj rP 0 Þ þ ðgrad wÞP j ðrj rP j Þ; 2 2
ð13Þ
where rj is the position vector of the cell-face center j. 3.2.2. Temporal distribution According to the adopted fully implicit time discretisation scheme, diffusion and source terms are evaluated at the current instant of time tm, where m is time step counter. Henceforth, the time step counter will normally be omitted, and all values of ui and p will refer to the current time tm, unless indicated otherwise. 3.2.3. Rate of change Using the backward differencing in time and the midpoint rule the transient term in Eq. (9) is approximated as Z ðqV ÞP 0 o oui dV q ðui 2um1 þ uim2 ÞP 0 ; ð14Þ i ot V ot dt2 where the volume integral is approximated using the midpoint rule, and the time increment dt = tm tm1. 3.2.4. Diffusion Using the midpoint rule for approximation of the surface integral, the diffusion flux of ui through a cell-face j can be approximated by Z l grad ui ds lj ðgrad ui Þj sj ; ð15Þ Sj
where lj stands for the cell-face mean value of the diffusivity, obtained by using Eq. (13), and the face value of (grad ui) is constructed in the following manner: " # ðui ÞP j ðui ÞP 0 grad ui dj jdj jsj ðgrad ui Þj ¼ ðgrad ui Þj þ ; ð16Þ jdj j jdj j dj s j where the over bar denotes the arithmetic average of the values calculated at nodes P0 and Pj. The second term on the righthand side (term in [ ] brackets) represents the difference between the second-order central difference approximation of the derivative in the direction of vector dj and the value obtained by interpolating cell-center gradients. It vanishes if the spatial variation of ui is linear or quadratic. Otherwise, its magnitude is proportional to the second-order truncation error of the scheme and reduces accordingly with grid refinement. This correction term detects and smoothes out any unphysical
6382
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
oscillations that might occur in the iteration process when solving resulting set of algebraic equations. The first part of this term is treated implicitly, and the rest of the diffusion flux explicitly. 3.2.5. Source terms The surface and volume integrals in Eq. (9) are calculated using again the midpoint rule: Z ½lðgrad uÞii pii ds ½lðgrad uÞii pii j sj
ð17Þ
Sj
and
Z V
qfbi dV ðqfbi ÞP 0 V P 0 .
ð18Þ
3.2.6. Initial and boundary conditions To start the calculation, the initial values u0i and p0, and in case of a transient run velocities v0i are required. On the faces coinciding with the boundary of the solution domain, boundary conditions have to be applied. In case of Dirichlet boundary condition, the expressions for diffusion fluxes (15) and sources (17) and (18) remain valid, except that ðui ÞP j is replaced by the boundary value (ui)B. On the boundary regions where Neumann boundary conditions are specified, the boundary fluxes are added to the source term, while the variable values at the boundary are obtained by extrapolating displacement values from the interior of the solution domain. 3.2.7. Resulting algebraic equations After assembling all terms featuring in Eq. (9), for each CV a following algebraic equation, which links the value of dependent variable ui at the control volume center with its values at the points in the neighborhood, is obtained: aui0 ðui ÞP 0
ni X
auij ðui ÞP j ¼ bui ;
ð19Þ
j¼1
where ni is the number of internal faces of the cell P0, sj sj aui j ¼ lj ; dj s j nf X ðqV ÞP 0 auij þ ; aui 0 ¼ dt2 j¼1 X nf nf X sj sj bui ¼ lj ðgrad ui Þj sj grad ui dj ½lðgrad uÞ ii pii j sj þ dj s j j¼1 j¼1 nB X ðqV ÞP 0 m1 m2 þ ðqfbi V ÞP 0 þ 2u u þ auiB ðui ÞB i i P0 dt2 B¼1
ð20Þ
and nB = nf ni is the number of boundary faces surrounding cell P0. 3.3. Discretisation of the volumetric constraint equation and calculation of pressure It can be noted that the pressure which appears in the source term of the discretised momentum equation is unknown, while at the same time Eq. (8) describing kinematic condition of incompressibility has not been used yet. After the space discretisation, Eq. (8) is written for each control volume as follows: nf Z nf Z X X u ds uj sj ¼ 0. ð21Þ j¼1
Sj
j¼1
Sj
In order to find displacement vector components ui and pressure p from the discretised momentum equation (9) and volumetric constraint equation (21) by an iterative solver a spatial solution strategy is adopted, which incorporates a pressurecorrection method of the SIMPLE-type [12]. According to the collocated version of the SIMPLE algorithm [18], the displacement at the cell face uj is constructed in the following manner: V P 0 pP j pP 0 grad p dj jdj jsj uj ¼ uj ; ð22Þ au0 jdj j jdj j dj s j
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
6383
where uj is the spatially interpolated displacement, the over bar denotes the arithmetic average, and au0 is the corresponding momentum equation central coefficient. The second term in Eq. (22) is a diffusion term, analogous to the term introduced by Eq. (16), when the diffusive transport of variable ui was discussed. Its role is to smooth out the oscillatory pressure profile, and at the same time introduces the pressure into the volumetric constraint equation in such a manner that a pressurecorrection equation can easily be constructed. This is achieved by employing the predictor–corrector procedure which will be briefly outlined here. The so-called predictor stage values of u and p (featuring in expression (22) for uj ), which satisfy the momentum equation, do not necessarily satisfy the volumetric constraint equation (21). By correcting the cell-face velocities, e.g., 0 0 V P 0 pP j pP 0 jdj jsj uj ¼ uj ð23Þ au0 jdj j dj sj and requiring that the double starred cell-face displacement satisfy the volumetric constraint equation (21), an equation for the pressure-correction p 0 is obtained: nf X ap00 p0P 0 ap0j p0P j ¼ bp0 ; ð24Þ j¼1
with coefficients: V P 0 sj sj ap0j ¼ ; au 0 dj s j
ap00 ¼
nf X
ap0j ;
j¼1
bp 0 ¼
nf X
uj sj ;
ð25Þ
j¼1
where all variables have their predictor stage values. After the field of pressure-correction p 0 is obtained, the displacement and pressure are corrected: uP 0 ¼ uP 0 ;pred þ uP 00 ¼ uP 0 ;pred pP 0 ¼ pP 0 ;pred þ
nf 1 X p 0 sj ; au0 j¼1 j
ð26Þ
bp p0P 0 ;
where bp is an under-relaxation factor (typically, bp = 0.2–0.3). The boundary conditions for the pressure-correction equation (24) depend on the boundary conditions for the momentum equations. On those portions of the boundary where displacement is prescribed, displacement correction is zero, which implies a zero-gradient Neumann boundary condition on the pressure correction (see Eq. (23)). If the traction is prescribed at the boundary, the pressure is calculated via constitutive relation and its correction is zero, leading to a Dirichlet boundary condition for the pressure correction. 3.4. Solution algorithm After assembling Eqs. (19) and (24) for all CVs, four sets of N coupled linear algebraic equations is obtained, where N is the number of CVs. The solution of this system of algebraic equations is obtained by employing a segregated iterative algorithm. It consists of temporary decoupling of Eqs. (19) and (24) by assuming that coefficients are known (calculated by using dependent variable values from the previous iteration or the previous time step). As a result, a set of linear algebraic equations for each dependent variable is obtained: A/ / ¼ b / ;
ð27Þ 0
where / stands for ui or p , A/ is an N · N matrix, vector / contains values of dependent variable / at N nodal points and b/ is the source vector. The matrices A/ resulting from the discretisation method described above are sparse, with the number of non-zero elements in each row equal to the number of nearest neighbors plus one. They are also diagonally dominant, which makes the equation system (27) easily solvable by an iterative method which retains the sparsity of the matrix A. Eq. (27) are solved using the conjugate gradient method with an incomplete Cholesky preconditioning. There is no need to solve them to a tight tolerance, since the coefficients and sources are only an approximation (based on the values of dependent variables from the previous iteration or time step) and reduction of the sum of absolute residuals by an order of magnitude normally suffices. Note that in order to promote stability of the solution method, an under-relaxation is often necessary. 4. Test cases In this section a set of test cases chosen to demonstrate capabilities of the method is presented. The results of calculations are compared with analytical solutions and finite element predictions.
6384
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
4.1. Inflation of an infinitely long tube This is a standard test problem designed to identify the volumetric locking and existence of the stress oscillation in a numerical solution since the elasticity of this problem is dominated by pressure, and the constraints in axial direction further intensifies the numerical difficulty in the incompressibility limit. An infinitely long tube, with inner radius of 8 m, an outer radius of 16 m is subjected to an internal pressure p = 10 MPa. Material of the tube is assumed incompressible and linear elastic with Young’s modulus E = 1000 MPa. The analytical solution for the radial and hoop stress, and radial displacement field is [19] pi r2i r2o pi r2i r2o 1 pi r2i r2o rrr ¼ 2 ; ð28Þ 1 2 ; rhh ¼ 2 1 þ 2 ; ur ¼ 2 2 2l r2o r2i r ro ri r ro ri r where pi is internal pressure, l is the shear modulus and ri and ro are the inner and the outer radius of the tube, respectively. The problem is considered as a plane strain, with a space domain made of 10 segment of the tube discretised with two different grids made of quadrilateral and triangular cells, as shown in Fig. 2. Fig. 3 (left) shows numerical and analytical results of the radial displacement along the radius of the tube. Numerical results are shown for two grids, one made of quadrilateral and the other of triangular cells, as shown in Fig. 2. It can be seen that the both finite volume profiles follow exactly the analytical ones. In Fig. 3 (right) the relative displacement error norm: Relative error ¼
ku uh k2 ; kuk2
ð29Þ
where u and uh are the exact and numerical displacement vector solutions, is plotted together with an ideal slope of the second-order scheme as a function of the (quadrilateral) cell number in the radial direction of the tube, demonstrating the second-order accuracy of the present method. Fig. 4 shows numerical and analytical results of the radial and hoop stress distribution along the radius of the tube. The numerical results correspond to the numerical mesh made of quadrilateral cells shown in Fig. 2. It can be seen that the finite volume results show a good agreement with the analytical ones, and that there is no indication of any stress oscillation. 4.2. Cook’s membrane problem To demonstrate the performance of the present discretisation scheme in a bending dominated response, a tapered panel clamped on one and subjected to a constant in-plane shearing stress s = 0.0625 MPa on the free end, as shown in Fig. 5, is
0.16
10.00
0.12
1.00
Relative error
Radial displacement (m)
Fig. 2. Computational mesh made of quadrilateral (left) and triangular cells (right).
0.08
0.10
analytical FVM (quadrilateral cells) FVM (triangular cells)
0.04 8
10
12 Radius (m)
14
16
0.01
1
10 Number of cells in radial direction
100
Fig. 3. Radial displacement vs. radius of the tube (left) and relative displacement error norm vs. number of cells in radial direction (right). The dashed line represents the theoretical slope for a second-order method.
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390 0
6385
20
Hoop stress (MPa)
Radial stress (MPa)
16 -4
-8 analytical FVM
12 8 analytical FVM
4
-12
0 8
10
12 Radius (m)
14
16
8
10
12 Radius (m)
14
16
Fig. 4. Radial (left) and hoop stress distributions (right) vs. radius of the tube.
τ 16 m
44 m
44 m
48 m Fig. 5. Cook’s membrane problem sketch (left) and a finite volume mesh (right).
100
10 Tip displacement error (%)
Vertical top corner displacement x 102 (mm)
considered. In the context of linear elasticity, this simulation is commonly referred to as Cook’s membrane problem. The problem is complicated as both shear and bending are present in combination with a skew numerical mesh. The material is considered incompressible and linear elastic with Young’s modulus E = 206.9 MPa. Plane strain deformation of the body is assumed. The space domain is discretised using a number of systematically refined numerical grids made of quadrilateral cells. The coarsest mesh is shown in Fig. 5. Since the analytical solution for this problem is not known, in Fig. 6 the numerical results for the vertical top corner displacement of the membrane as a function of the number of cells/elements per side are compared with the finite element calculations presented in [20] obtained with a Q1/ME2 element which is two-dimensional counterpart of the three-dimensional mixed-enhanced formulation with nine enhanced modes, standard tri-linear interpolation functions, and standard 8-pt quadrature rule. It can be seen a very close agreement between the two solutions.
8 6 4 FVM
2
FEM (ν = 0.4999)
0
10
1
0.1 0
10
20
30
40
50
Number of cells per side
60
70
1
10 100 Number of cells per side
1000
Fig. 6. Vertical top corner displacement (left) and corresponding error (right) vs. number of cells for the Cook’s membrane problem. The dashed line represents the theoretical slope for a second-order method.
6386
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390 y
D
C
4m
Solution domain
r E
Θ
1
x A
σ0
B 4m
σ0
Fig. 7. Infinite plate with a circular hole subjected to unidirectional tension (left) and computational mesh (right).
To check the order of accuracy of the present discretisation method, the error of the vertical top corner displacements together with an ideal slope of the second-order scheme is also plotted as a function of number of cells per side of the membrane in Fig. 6. It can be seen that the convergence towards the analytical solution is asymptotically of the second order. The analytical solution of the problem is evaluated from the numerical solution using the Richardson’s extrapolation. 4.3. Infinite plate with a circular hole An infinite plate with a circular hole under uniform tension as depicted in Fig. 7 is analysed. The problem is considered as a plane strain and fully incompressible material with Young’s modulus E = 1000 MPa is assumed. Normal stresses r0 = 1 MPa are uniformly distributed along two edges of the plate. According to the analytical solution [19] the stress field corresponding to this problem is a2 3 3 a4 cos 2h þ cos 4h þ rxx ¼ r0 1 2 cos h ; 2 r4 r 2 2 a 1 3 a4 ð30Þ cos 2h cos 4h cos 4h ; ryy ¼ r0 2 2 r4 r 2 2 a 1 3 a4 sin 2h þ sin 4h þ rxy ¼ r0 2 sin 4h ; 2 r4 r 2 and the displacement field is r0 a r a a3 uðr; hÞ ¼ cos h þ ð2 cos h þ cos 3hÞ 3 cos 3h ; 4l a r r r0 a r a a3 sin h þ sin 3h 3 sin 3h ; vðr; hÞ ¼ 4l a r r
ð31Þ
where a is the hole radius and l is the shear modulus. Taking into account the symmetry of the problem, the solution domain shown in Fig. 7 is used for the numerical analysis. In order to eliminate the influence of the finite plate dimensions the traction calculated from the analytical solution is imposed at the boundaries BC and CD. The symmetry boundary conditions are applied at the boundaries AB and DE, while zero traction is prescribed at the hole boundary EA. The polyhedral grid employed for numerical calculations is presented in Fig. 7. In Fig. 8 analytical and numerical results of displacement profile along the boundary AB and stress rxx profile along the boundary ED are presented. The numerical values of stresses are calculated in the computational points at the geometrical centers of the cells next to the boundary ED. It can be observed that numerical results show a very good performance of the method. 4.4. Force acting on a semi-infinite solid To test the method on a three-dimensional problem the well-known Boussinesq’s problem of an half-space with a point load is considered. The problem is approximated by analysing a cubic volume shown in Fig. 9. According to the analytical solution [22], in case of incompressible media the displacement field and normal stress in the direction of the load corresponding to this problem are
4
4
3
3
6387
analytical FVM
σxx (MPa)
Radial displacement x 103 (m)
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
2
analytical
1
2
1
FVM
0
0 1
2
3
4
1
Radial distance (m)
2
3
4
Radial distance (m)
Fig. 8. Displacement profile along the boundary AB (Fig. 7) and stress profile along the boundary ED of the plate.
z
1
1
1
y
x
Fig. 9. Problem sketch for the Boussinesq’s problem (left) and computational grid (right).
P z2 1þ 2 ; 4plq q
3Pz3 ; ð32Þ 2pq5 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where P is the force magnitude, l is the shear modulus, and q ¼ x2 þ y 2 þ z2 . In a numerical analysis of the problem material of the body is assumed incompressible with Young’s modulus E = 1000 MPa. At planes x = 1 m, y = 1 m and z = 0 displacements are prescribed according to the analytical solution, and at planes x = 0 and y = 0 symmetry plane boundary conditions are prescribed. The point load with a magnitude P/4 = 1000 kN is approximated with an equivalent pressure load distributed over the face of the cell next to the point load, while the rest of the plane z = 1 is traction free. The numerical mesh used in calculation is shown in Fig. 9 (right). Fig. 10 shows numerical and analytical results of the displacement and normal stress component in the direction of the load along the z axis. It can be seen that the numerical calculations show a good agreement with the analytical results. u¼
P xz ; 4pl q3
v¼
P yz ; 4pl q3
w¼
rzz ¼
4.5. Epoxy disc enclosed in thin steel ring To show a robustness of the method, an epoxy disc enclosed in a thin steel ring as depicted in Fig. 11 is considered. The star shaped hole at the center has six symmetrically located leaflets whose dimensions are shown in Fig. 11. A uniform pressure of 6.894 MPa inside the star shaped hole is applied and 25.4 mm long cracks in the radial direction at the tip of the leaflets are inserted. A plane strain deformation of this composite cylindrical body is assumed and, because of the symmetry of the geometry and the loading conditions, deformation of only a 30 segment of the body is investigated. In planes h = 0 and h = 30 the symmetry plane boundary conditions are imposed, and the outer surface of the steel ring is taken to be traction free. The epoxy disc is assumed to be perfectly bonded to the steel ring so that the displacement and the surface traction are continuous across their common interface. Material of the epoxy is modelled as a totally incompressible Hookean with the
6388
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390 0.0E+00
0.0E+00
Displacement w (m)
-6.0E+02 -4.0E-02 -1.2E+03 -8.0E-02 -1.8E+03 analytical
analytical
-1.2E-01
-2.4E+03
FVM
-1.6E-01 0.00
0.25
0.50
0.75
1.00
FVM
-3.0E+03 0.00
0.25
z-coordinate (m)
0.50
0.75
1.00
z-coordinate (m)
Fig. 10. Displacement (left) and normal stress component (right) profile in the direction of the load along the z axis.
12.7 816.4 30o 116.6 389
Epoxy
Steel
Fig. 11. A cross-section of the epoxy disc enclosed in a thin steel ring (dimensions in millimeters) (left) and computational mesh (right).
y-coordinate (mm)
shear modulus equal to 10.19 MPa, while the steel ring is assumed to obey the Hooke’s law with Young’s modulus E = 200,000 MPa and Poisson’s ratio m = 0.3. The finite volume calculations are performed using the mesh locally refined around the leaflet corner and the crack, consisting of 2844 control volumes, as shown in Fig. 11 and results are compared with finite element results obtained with the commercial code ABAQUS 6.11 in Ref. [21], where the space domain was divided into 4824 CPE8H (8-node element with biquadratic interpolation for displacement and linear interpolation for the pressure field), and the epoxy was modelled as a Hookean with the shear modulus equal to 10.19 MPa and the Poisson’s ratio was varied from 0.49 to 0.4999. Deformed shapes of the crack face calculated by finite volume and finite element method (FEM) are plotted in Fig. 12. The finite element results are shown for two different Poisson’s ratios, 0.49 and 0.4999. It can be seen that crack opening displacements decrease with an increase in Poisson’s ratio and that curvature of the deformed crack at the tip increases with an increase in Poisson’s ratio. Finite volume calculations agree very well with finite element results for Poisson’s ratio of 0.4999. Finite volume and finite element calculations of strain components xx and yy vs. radial distance from the crack tip at points ahead of the crack front are presented in Fig. 13. It seems that the FEM solution exhibits non-physical oscillations 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 390
FVM
crack tip
395
400
405 410 415 x-coordinate (mm)
420
Fig. 12. Deformed shapes of the crack face.
425
0.10
0.10
0.08
0.08
εxx εyy
0.06
Strain components
Strain components
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
0.04 0.02 0.00 -0.02 -0.04
6389
εxx εyy
0.06 0.04 0.02 0.00 -0.02 -0.04
-0.06 0
1
2
3
4
5
Distance r (mm)
6
7
-0.06 0
1
2
3
4
5
6
7
Distance r (mm)
Fig. 13. Finite element (left) and finite volume (right) calculation of strain components vs. radial distance from the crack tip at points ahead of the crack front.
in the vicinity of the crack tip up to about 5 mm from the crack tip. Besides, it is evident that displacement field is not divergence free, thus violating incompressibility constrain equation in this region. Namely, for an incompressible material and plane strain deformation, from incompressibility constrain equation follows that xx must equal to yy. On the other hand, the finite volume solution does not show oscillatory behavior near the crack tip. It can be also concluded that finite volume solution produces a divergence free displacement field with an exception of the first computational point in the cell next to the crack tip. 5. Conclusions In this paper a finite volume based numerical method for predicting the linear elastic behavior of incompressible elastic solid is successfully applied. The method solves the momentum balance and volumetric constraint equations written in an integral form. The incompressibility constraint is enforced by employing the hydrostatic pressure as an additional variable. A segregated approach employed to solve resulting set of coupled algebraic equations, embedding a SIMPLE based algorithm for calculation of hydrostatic pressure featuring in the constitutive equation shows a good convergence behavior. The presented test cases demonstrate very good accuracy of the method. The volumetric locking phenomenon, common in numerical description of an incompressible material behavior, is not registered. This is due to the conservative nature of the method which ensures that the volumetric incompressibility constraint is enforced exactly for each cell as well as for the whole body. Another phenomenon, non-physical pressure oscillations, which are also inherent in numerical analysis of incompressible body, are effectively eliminated embedding a third-order pressure diffusion term into the SIMPLE based algorithm for displacement–pressure coupling. The present method has been extended to the large strain analysis of incompressible hyperelastic materials [23], and work is in progress to apply it to compressible solids. References [1] L.R. Herrman, Elasticity equations for incompressible and nearly incompressible materials by a variational theorem, AIAA J. 3 (1965) 1896–1900. [2] T. Belytschko, W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy by a variational theorem, Comput. Methods Appl. Mech. Engrg. 54 (1986) 279–301. [3] T.S. Sussman, K.J. Bathe, A finite element formulation for incompressible elastic and inelastic analysis, Comput. Struct. 26 (1987) 357–409. [4] D.S. Malkus, T.J.R. Hughes, Mixed finite element methods—reduced and selective integration techniques: a unification of concepts, Comput. Methods Appl. Mech. Engrg. 15 (1978) 63–81. [5] T. Belytschko, J.S.J. Ong, W.K. Liu, J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comput. Methods Appl. Mech. Engrg. 43 (1984) 251–276. [6] T.J.R. Hughes, Generalization of selective integration procedures to anisotropic and nonlinear media, Int. J. Numer. Methods Engrg. 15 (1980) 1413– 1418. [7] J.C. Simo, R.L. Taylor, K.S. Pister, Variational and projection methods for volume constraint in finite deformation elasto-plasticity, Comput. Methods Appl. Mech. Engrg. 51 (1985) 177–208. [8] J.S. Chen, C. Pan, C.T. Wu, A pressure projection method for nearly incompressible rubber hyperelasticity, J. Appl. Mech. 63 (1996) 862–876. [9] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, New Jersey, 1987. [10] J. Dolbow, T. Belytschko, Volumetric locking in the element free Galerkin method, Int. J. Numer. Methods Engrg. 46 (1999) 925–942. [11] F.S. Henry, M.W. Collins, Prediction of transient wall movement of an incompressible elastic tube using a finite volume procedure, in: Proceedings of the Second International Conference on Computers in Biomedicine, BIOMED 93, UK, 1993. [12] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980.
6390
I. Bijelonja et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 6378–6390
[13] F.S. Henry, M.W. Collins, A novel predictive model with compliance for arterial flows, in: Proceedings of the 1993 ASME Winter Annual Meeting, New Orleans, 1993. [14] M.A. Wheel, Modelling incompressible materials using a mixed structural finite volume approach, in: Computational Mechanics in UK 1997, ACME, vol. 5, London, 1997. [15] M.A. Wheel, A mixed finite volume formulation for determining the small strain deformation of incompressible materials, Int. J. Numer. Methods Engrg. 44 (1999) 1843–1861. [16] I. Demirdzˇic´, S. Muzaferija, Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology, Comput. Methods Appl. Mech. Engrg. 125 (1995) 235–255. [17] I. Demirdzˇic´, S. Muzaferija, Finite volume method for stress analysis in complex domains, Int. J. Numer. Methods Engrg. 37 (1994) 3751–3766. [18] C.M. Rhiw, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (1983) 1525–1532. [19] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1970. [20] E.P. Kasper, R.L. Taylor, A mixed-enhanced strain method, part i: geometrically linear problems, Comput. Struct. 75 (2000) 237–250. [21] R.C. Batra, H.K. Ching, Energy release rates in a constrained epoxy disc with Hookean and Mooney-Rivlin materials, Theor. Appl. Fracture Mech. 38 (2002) 165–175. [22] A.S. Saada, Elasticity: Theory and Applications, Pergamon Press, New York, 1974. [23] I. Bijelonja, I. Demirdzˇic´, S. Muzaferija, A finite volume method for large strain analysis of incompressible hyperelastic materials, Int. J. Numer. Methods Engrg. 64 (2005) 1594–1609.
