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

    *                                                                  * 

    *  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.

    * The divergence-of-gradient operator: 

                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 

    coordinate directions lead to:

    

      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.

        

    Re-arrangement so as to give prominence to the div_grad term, leads 

    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.