```

********************************************************************

*                                                                  *

*  The equations of displacements, stresses and strains in solids  *

*                                                                  *

********************************************************************

by: D B Spalding

Contents

--------

1. Nomenclature

2. Strain related to stress

3. Normal stress related to strain

4. Balance of forces

5. Differential equations for the displacements

6. The differential equations for the displacements; 2D

7. The differential equations for the displacements; 1D

8. Cylindrical geometry

9. The representation of the differential equations in PHOENICS

10. Concluding remarks

--------------------------------------------------------------------

1. Nomenclature

---------------

* Cartesian coordinate directions:                  x, y, z

* Radius, axial distance and angle in

cylindrical-polar cordinates :  r, a, z

* Strains (extensions) in x, y, and z directions:   ex, ey, ez

* Linear thermal expansion:                         et

* Normal stresses in x, y, z directions, divided by Young's modulus,

presumed uniform:                                 nx, ny, nz

* Shear stresses in xy, xz, yz planes, divided by Young's modulus:

sxy, sxz, syz

* Displacements in x, y and z directions:           u, v, w

* Poisson's ratio:                                  R

* Various constant functions of Poisson's ratio:  R1, R2, R3, R4 ...

namely:

R1  =  (1 - R) / [(1  + R) (1 - 2 R)]

R2  =  R / [(1 + R) (1 - 2 R)]

R3  =  1 / (1 - 2 R)

R4  =  1 / [ 2 (1 + R)]

* derived functions:

R5  =  R1  -  R2 R2 / R1          = 1 / (1 + R)

R6  =  R2  -  R2 R2 / R1          = R / [(1 + R) (1 - R)]

R7  =  R3  -  R3 R2 / R1          = 1 / (1 - R)

R8  =  R1 - 2 R2 R2 / ( R1 + R2) ]

R9  =  R3 [ 1 - 2 R2 R2 / ( R1 + R2) ]

R10 =  1 / [ 2 ( 1 + R ) ( 1 - 2 R) ]

R11 =  1 / R4                     = 2 (1 + R)

R12 =  2 ( 1 + R ) / ( 1 - R )

R13 =  2 ( 1 + R ) (1 - R )

R14 =  2 ( 1 - R )

Note that, if R = 0, then R1 = R3 = R5 = R8 = R11 = 1.0,

R2 = R6 = 0.0 ,

R4 = R10 = 0.5  and

R12= 2.0

* First-order differentiation operators

with respect to x, y, z, a:  ddx., ddy., etc

* Second-order differentiation operators with respect to

x & x   y & y   z & z   x & y   x & z  and  y & z:

d2dxx.  d2dyy.  d2dzz.  d2dxy.  d2dxz.      d2dyz.

div_grad.= d2dxx. +  d2dyy.  + d2dzz

* The multiplication operator is implied by any space in an equation

which is not adjacent to one of the operators: + - / = or a

differential operator.

* External forces per unit volume, including gravity, aceleration,

etc, in the x, y and z directions:  fx, fy , fz

2. Strain related to stress

---------------------------

The following are expressions of Hooke's and Poisson's

experimentally-based laws:

2.1 Normal strains related to normal stresses

---------------------------------------------

ex =  nx   -   R ( ny + nz )  +  et                (2.1-1)

ey =  ny   -   R ( nx + nz )  +  et                (2.1-2)

ez =  nz   -   R ( nx + ny )  +  et                (2.1-3)

[ These correspond to the top part of matrix D on page 58 of Huebner &

Thornton (H&T form now on), except that et (thermal strain) is not

mentioned by them.

In the absence of any stresses, they correspond to ex=ey=ez=et, which

seems reasonable. ]

2.2 Shear stresses related to shear strains

-------------------------------------------

sxy = R4 ( ddy.u  + ddx.v )               = syx    (2.2-1)

sxz = R4 ( ddz.u  + ddx.w )               = sxz    (2.2 2)

szy = R4 ( ddy.w  + ddz.v )               = szy    (2.2-3)

[ The origin of these equations is not clear. They appear to correspond

to the bottom part of the C matrix of H&T. But the coefficient s there

are (1-2R)/2, whereas R4 = 1/{2(1+R)}  ]

3. Normal stress related to strain normal strain

------------------------------------------------

3.1 General 3D equations

------------------------

The following are derived from the above by algebraic

re-arrangement:

nx = R1 ex  +  R2 ( ey + ez )  -  R3 et        (3.1-1)

ny = R1 ey  +  R2 ( ex + ez )  -  R3 et        (3.1-2)

nz = R1 ez  +  R2 ( ex + ey )  -  R3 et        (3.1-3)

The details of the derivation, for nx only, are as follows:

Addition of R times equations 2.1-2 and 2.1-3 to

(1 - R) times equation 2.1-1 leads to:

(1 - R) ex   +   R ey   +   R ez

=

(1 - R) nx  -  (1 - R) R (ny + nz)  +  (1 - R) et

+     R ny  -        R R (nx + nz)  +        R et

+     R nz  -        R R (nx + ny)  +        R et    (3.1-4)

in which the multipliers of both ny and nz sum to zero.

Hence,

(1 - R) ex   +   R ey   +   R ez  =

(1 - R - 2 R R) nx      + (1 + R) et                 (3.1-5)

ie

(1 - R) ex   +   R ey   +   R ez  =

(1 + R) (1 - 2 R) nx    + (1 + R) et                 (3.1-6)

ie

nx  = [ (1 - R)/ {(1  + R) (1 - 2 R)} ] ex  +

[ R / {(1 + R) (1 - 2 R)}] ( ey + ez) -

[ 1 / (1 - 2 R)} ] et                          (3.1-7)

ie

nx = R1 ex  +  R2 ( ey + ez )  -  R3 et             (3.1-1)

[ Note that the shear-stress equations (1.2-1.2.3) have not been used

here.]

3.2 2D equations

----------------

If ez (for example) is zero everywhere, which implies that the

process is two-dimensional. movement in the third dimension being

prevented, then the equations for nx and ny become:

nx = R1 ex  +  R2 ey  -  R3 et            (3.2-1)

ny = R1 ey  +  R2 ex  -  R3 et            (3.2-2)

nz = R2 ( ex + ey )   -  R3 et            (3.2-3)

while the equation for nz remains as above.

If however there is no constraint in the  z direction, ez may be

non-zero, while nz obeys:

nz = 0  .                                 (3.2-3)

Then ez may be eliminated from the equations of section 3.1, with

the result:

nx = R5 ex  +  R6 ey  -  R7 et            (3.2-4)

ny = R5 ey  +  R6 ex  -  R7 et            (3.2-5)

The derivation, for nx only, is as follows:

From 3.1-3, with nz = 0, leads to:

ez =  - (R2/R1) ( ex + ey )  +  (R3/R1) et    (3.2-6)

Substitution of ez, from 3.2-6 into 3.1-1, yields:

nx = R1 ex  +  R2 ( ey - (R2/R1) (ex + ey))

+  (R2 R3 /R1)  et - R3 et

ie

nx = (R1 - R2 R2/R1) ex  + R2 (1 - R2/R1) ey

- R3 (1 - R2/R1) et   (3.2-7)

ie, from the definitions of R5, R6 and R7:

nx = R5 ex  +  R6 ey  -  R7 et                 (3.2-4)

Expressed in terms of R, this is:

nx = [ 1 / (1 + R) ] ex  +

{ R / [(1 + R) (1 - R)] } ey  -

[ 1 / (1 - R) ]  et                 (3.2-4)

3.3 1D equations

----------------

If ey and ez (for example) are zero everywhere, the equations for

nx, ny and nz become:

nx = R1 ex  -  R3 et                        (3.3-1)

ny = R2 ex  -  R3 et                        (3.3-2)

nz = R2 ex  -  R3 et                        (3.3-3)

If however there are no constraints in the y and z direction, ey and

ez may be non-zero, while ny and nz obey:

ny = 0  and

nz = 0  ;

then ey and ez may be eliminated from the equations of section 3.2,

with the result:

nx = R8 ex  -  R9 et           (3.3-4)

The derivation is as follows. Equations 3.1-1. 3.1-2 and 3.1-3

become:

nx = R1 ex  +  R2 ( ey + ez )  -  R3 et        (3.3-5)

0 = R1 ey  +  R2 ( ex + ez )  -  R3 et        (3.3-6)

0 = R1 ez  +  R2 ( ex + ey )  -  R3 et        (3.3-7)

Since there is nothing to distinguish ey from ez, they may be taken

as equal, with the results:

nx = R1 ex  +  2 R2 ey  -  R3 et               (3.3-5)

0 = ( R1 + R2 ) ey  +  R2 ex  -  R3 et        (3.3-6)

Substitution for ey from 3.3-6 in 3.3-5 yields:

nx = R1 ex  +  [ 2 R2 / (R1 + R2) ] [ - R2 ex + R3 et]  -

R3 et                                     (3.3-7)

ie

nx = [ R1 - 2 R2 R2 / ( R1 + R2) ] ex  -

R3 [ 1 - 2 R2 R2 / ( R1 + R2) ] et      (3.3-8)

In terms of R, this is:

4. Balance of forces

--------------------

The balances of the forces acting on unit volume in the three

ddx. nx  +  ddy. sxy  +  ddz.sxz  +  fx  =  0           (4.-1)

ddy. ny  +  ddx. sxy  +  ddz.syz  +  fy  =  0           (4.-2)

ddz. nz  +  ddx. sxz  +  ddy.syz  +  fz  =  0           (4.-3)

5. The differential equations for the displacements; 3D

-------------------------------------------------------

5.1 The x-direction displacement u

----------------------------------

Differentiation of the equation 4.-1 leads,

after substitution from the equations of section 3.1, to:

ddx. [ R1 ex  +  R2 ( ey + ez )  -  R3 et ]  +

ddy. [ R4 ( ddy.u  + ddx.v) ]  +

ddz. [ R4 ( ddz.u  + ddx.w) ]  +

fx  =                                            0.0

Since, by definition:

ex  =  ddx. u

ey  =  ddy. v

ez  =  ddz. w

the strains ex, ey, ez can be eliminated from the above differential

equation, with the result:

ddx. [ R1 ddx. u  +  R2 ( ddy. v + ddz. w )  -  R3 et ]  +

ddy. [ R4 ( ddy. u  + ddx. v) ]  +

ddz. [ R4 ( ddz. u  + ddx. w) ]  +

fx                                     =  0.0

If Poisson's ratio R is taken as constant, implying that all Rs are

also constant, this equation can be written in terms of second-order

operators, as follows:

ddx. [ R1 ddx. u  +  R2 ( ddy. v + ddz. w )  -  R3 et ]  +

R4 d2dyy. u  +  R4 ddx.  ddy. v  +

R4 d2dzz. u  +  R4 ddx.  ddz. w  +

fx  =  0.0

with further rearrangement to:

ddx. [ ( R1 - R4 ) ddx. u + ( R2 + R4 ) ( ddy. v + ddz. w ) ] -

ddx. [ R3 et ]  +

R4 ( d2dxx. u  + d2dyy. u  + d2dzz. u ) +

fx  =  0.0

But, as a consequence of the above definitions,

R1 - R4  = R2 + R4  =  1 / [2 (1 + R) (1 - 2 R)]  =  R10

Therefore, the equation can be reduced to:

R10 ddx. dil  -  R3 ddx et  +  R4 div_grad. u +  fx  = 0

where

dil = ddx. u  +  ddy. v  +  ddz. w    , the dilatation.

to:

div_grad. u  +   ddx. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]

+   fx / R4  =  0.0                            (5.1-1)

[ This should be compared with equation C.15 of H&T. In that equation,

the multiplier of dil is 1/(1-2R).

R10/R4 = 2(1+R)/{2(1+R)(1-2R)} i.e. 1/(1-2R) . So they equations agree.

The multiplier of fx in H&T is 2(1+R)/E, which is again in agreement.

However, H&T are silent about  et .]

5.2 The y- and z-direction displacements v and w

------------------------------------------------

Replacing x successively by y and z, and u by v and w, leads to the

corresponding equations:

div_grad. v  +  ddy. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]

+  fy / R4  =  0.0

div_grad. w  +  ddz. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]

+  fz / R4  =  0.0

Note that, as a consequence of the definitions,

R10 / R4  =  1 /  ( 1.0  - 2 R )          =  R3

and

R3 / R4   =  2 ( 1 + R ) / ( 1 - 2 R )

The three equations can thus be written as:

--------------------------------------------------------------------

|u|           |x|                            |x|

div_grad. |v|  +   R3 dd|y|. ( dil - R11 et )  +  R11 f|y|  =  0.0

|w|           |z|                            |z|

(5.2-1)

--------------------------------------------------------------------

where it is to be noted that:

R3  = 1 /( 1 - 2 R)

R11 = 1 / R4         = 2 ( 1 + R )

6. The differential equations for the displacements; 2D

-------------------------------------------------------

6.1 When displacements in the third direction are zero

------------------------------------------------------

Let the zero-displacement direction be z. Then the equation of

section 5.2 still holds, for directions x and y, with the proviso

that ddz. terms are absent from both the div_grad and dil

definitions.

6.2 When stressess in the third direction are zero

--------------------------------------------------

When nz, say, is zero, the analysis must start from

nx = R5 ex  +  R6 ey  -  R7 et

rather than from

nx = R1 ex  +  R2 ey  -  R3 et   .

Thereafter, apart from the substitutions of Rn+4 for Rn, and from

the absence of z-direction derivatives, the analysis proceeds as in

section 6. The resulting equations are:

--------------------------------------------------------------------

|u|          |x|                            |x|                               |

div_grad. |v|  +  R3 dd|y|. ( dil - R11 et )  +  R11 f|y|  =  0.0

(6.2-1)

--------------------------------------------------------------------

whereby it should be noted that R11 remains in place, but R13 has

replaced R12.

7. The differential equations for the displacements; 1D

-------------------------------------------------------

7.1 When displacements in the second and third directions are zero

------------------------------------------------------------------

Let the zero-displacement directions be y and z. Then the equation of

section 5.2 still holds, for directions x and y, with the proviso

that ddz. terms are absent from both the div_grad and dil

definitions.

7.2 When stressess in the second and third directions are zero

--------------------------------------------------------------

When ny and nz, say, are zero, the analysis must start from:

nx = R8 ex  +  R9 et  .

Thereafter, the analysis proceeds as in section 6. The resulting

equation is:

------------------------------------------------------------------

|

|  div_grad. u  +  R11 ddx. ( dil + R14 et )  +  fx / R4 =  0.0

|                                                            (7.2-1)

------------------------------------------------------------------

8. Cylindrical geometry

-----------------------

8.1 Nomenclature

----------------

When the grid is cylindrical polar, and x, y and z are replaced by

a (ie angle), r (ie radial distance) and z (now axial distance),

* the strains become:                                 ea, er, ez

* the normal stresses become:                         na, nr, nz

* the displacements become:                            u,  v,  w

8.2 Strain related to stress

----------------------------

The following relations replace those of section 2 above:

ea =  na   -   R ( nr + nz)  +  et

er =  nr   -   R ( na + nz)  +  et

ez =  nz   -   R ( na + nr)  +  et

sar = R4 ( ddr.u  -  u / r  +  dda.v / r) = sra

saz = R4 ( ddz.u  +  dda.w / r)           = sza

szr = R4 ( ddr.w  +  ddz.v )              = srz

8.3 Normal stress related to strain, general 3D

-----------------------------------------------

Because the e-to-n relations are the same as before (ie in section

2), so are the n-to-e relations (ie as in section 3).

8.4 Balance of forces

---------------------

In the angular direction, the balance of forces becomes:

dda. na / r  +  ddr. sar  +  ddz. saz  +  fa  +

( sra + sar ) / r  =  0.0      .

In the radial direction, the balance of forces becomes:

ddr. ( nr r ) / r  +  dda. sar / r  +  ddz. szr  -  na / r

+  fr   =  0.0        ,

which may also be expressed as:

ddr. nr  +  dda. sar / r  +  ddr. szr + ( nr - na ) / r

+  fr   =  0.0      .

In the axial direction, the balance of forces becomes:

ddz. nz  +  ddr. ( srz r) / r  +  dda. saz / r

+  fz  =  0.0       ,

which may also be written as:

ddz. nz  +  ddr. srz   +  dda. saz / r  +  srz / r

+  fz  =  0.0       .

Thus, as compared with the corresponding equations appearing section

3, the following extra terms appear:

-  in the equation for the angular (formerly x) direction:    none;

-  in the equation for the radial  (formerly y) direction:

+ ( nr - na ) / r   ;

-  in the equation for the axial (formerly, and still, z) direction:

+ srz / r           .

8.5 The differential equations for the displacements

----------------------------------------------------

Since the differences in the differential equations from their

cartesian-coordinate counterparts all derive from the above-

mentioned extra terms, it can be deduced that the differential

equations in cylindrical-polar coordinates can be expressed as:

div_grad. u  +  R11 dda. ( dil + R12 et ) / r +  fa / R4  =  0.0

div_grad. v  +  R11 ddr. ( dil + R12 et )

+  [fr +  ( nr - na ) / r ] / R4  =  0.0

div_grad. w  +  R11 ddz. ( dil + R12 et )

+  (fz + srz / r ) / R4  =  0.0

9. The representation of the differential equations in PHOENICS

--------------------------------------------------------------

9.1 Comparison of the equations for displacements and velocities

----------------------------------------------------------------

The momentum equations which are solved by PHOENICS, when the

viscosity is set to unity and the convection terms are neglected,

can be represented as:

----------------------------------------------------------

|u|           |x|                   |x|

div_grad. |v|    -    dd|y|. ( p )   +  mom_so|y|    =  0.0

|w|           |z|                   |z|            (9.1-1)

----------------------------------------------------------

where p stands for pressure, and mom_so stands for momentum source

per unit volume.

Comparison of these equations for those for the displacements shows

that they can be made identical if the pressure is related to the

dilatation and to the thermal expansion by:

-  p  =  R3 ( dil - R11 et )

ie     -  p  =  [1 / (1 - 2 R)] [ dil  -  2 ( 1 + R ) et ]

(9.1-2)

and if:

|x|  =       |x|

mom_so|y|  =  R11 f|y|

|z|  =       |z|

ie:

|x|  =               |x|

mom_so|y|  =  2 ( 1 + R ) f|y|                      (9.1-3)

|z|  =               |z|

9.2 The dilatation in PHOENICS

------------------------------

In PHOENICS, the dilatation, ie   ddx. u  +  ddy. v  +  ddz. w  ,

appears as a volumetric source in the mass-conservation equation.

Therefore the required relationship can be contrived by providing,

for the part of the domain which is occupied by solids, a volumetric

source term calculated as:

vol_sor = dil = - (1 / R3) p  +  R11 et

ie:

vol_sor = dil = - (1 - 2 R) p  +  2 (1 + R ) et     (9.2-1)

9.3 Normal-stress boundary conditions

-------------------------------------

The expression for nx in equation 3.1-1 can be re-written in terms

of the dilatation as follows:

nx = ( R1 - R2 ) ex  +  R2 ( ex + ey + ez )  -  R3 et

ie

nx = ( R1 - R2 ) ex  +  R2 dil  -  R3 et          (9.3-1)

Replacement of dil by p, from equation 9.1-2, then leads to:

nx = ( R1 - R2 ) ex  +

R2 [- (1 - 2 R) p  +  2 (1 + R ) et ]  -

R3 et                                        (9.3-2)

This can be rearranged with ex on the left-hand side and replaced by

its equivalent, ddx.u, as:

ddx. u = [ 1 / (R1 - R2) ] nx +

[ R2 (1 - 2 R) / (R1 - R2 )] p +

[ R3 - 2 R2 (1 + R) / ( R1 - R2 ) ] et   (9.3-3)

which becomes, after simplification:

ddx. u = ( 1 + R) nx  +  R p + ( 1 + R ) et       (9.3-4)

Since  ddx.u ie the gradient of the displacement u, this equation

can be used as the normal-stress boundary condition of the

differential equation for u.

It can be represented by way of a source term in the computational

cell within the solid and adjacent to the fluid, as illustrated

below.

.      |      .      |     .   /|

.      |      .      |     .   /|

.      |      .      |     .   /|

-----*------>------*------>-----*---->=====> normal stress

.      |      .      |     .   /|

.      |   solid     |  solid  /|   fluid

.      |      .      |     .   /|

Corresponding equations can be derived for the v and w boundary

conditions, namely:

ddy. v = ( 1 + R) ny  +  R p + ( 1 + R ) et       (9.3-5)

ddz. w = ( 1 + R) nz  +  R p + ( 1 + R ) et       (9.3-6)

A consequence is that, when a direct stress S is to be exerted at a

boundary, the RATCH and COVAL statements in the RHOENICS Q1 file

take the form:

PATCH(name,north/east/high,,,,,,)

COVAL(name,v1/u1/w1,fixflu,S*(1+R))

A further consequence is that, within EARTH, sources per unit area

are applied to cells adjacent to solid-fluid boundaries which are

equal to:

R p + ( 1 + R ) et     .

10. Concluding remarks

----------------------

Further topics to be dealt with in the continuation of the current

document include:-

* allowance for non-uniform material properties;

* the relationship of the dilatation to the PHOENICS pressure

variable;

* the solution algorithm;

* a series of exact solutions of the equations, which can be used as

tests;

* a series of practically-interesting problems, which can be used so

as to interest engineers in the use of the technique.

It should also be remarked that all the above derivations require to

be checked carefully and independently.

```