*                                                                  * 

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

    *                                                                  * 



    by: D B Spalding  (latest update 04.04.06)





    1. Nomenclature

    2. The differential equations for displacements

    3. Boundary conditions for displacements

      3.1 Forces and stresses at solid-fluid surfaces

      3.2 Boundary conditions involving displacement gradients

      3.3 Discussion

    4. The rotation equations

      4.1 Differential equations

      4.2 Boundary conditions

    5. Differential equations for displacements in terms of rotations

      5.1 Mathematical manipulations

      5.2 The resulting new equations

      5.2.1 Equations in terms of rotation gradients 

      5.2.2 Expressing the rotation-gradients via displacement gradients 

      5.2.3 Introducing the shear stress 

      5.2.4 The complete differential equation 

      5.2.5 Discussion 

      5.2.6 Finite-volume expression

      5.4 One- and two-dimensional situations

        5.4.1 The displacement gradients at constrained boundaries

        5.4.2 Gradients of ex, ey and ez at boundaries

    6. Implementation in PHOENICS

      6.1 Comparison of the differential equations

      6.2 The material-dependent properties

      6.3 The dilatation

      6.4 Applied-force boundary conditions

      6.5 Thermal-expansion boundary conditions

      6.6 Rotation-gradient sources

    7. Special cases

      7.1 Uniform direct stress

      7.2 Uniform thermal stress

      7.3 Non-uniform direct stress

      7.4 Non-uniform thermal stress

      7.5 Uniform bending stress

      7.6 Non-uniform bending stress

    8. Solutions in terms of 'stress-functions'

      8.1 Introduction

      8.2 The torsion of straight bars of non-uniform cross-section

      8.3 The bending of beams of non-uniform cross-section  



    1. Nomenclature


    * Literature reference: 

         SP Timoshenko and JN Goodier,

         Theory of Elasticity (3rd edition)

         McGraw Hill, 1970:                             T&G


         BA Boley and JH Weiner

         Theory of Thermal Stresses

         John Wiley, 1960                               B&W


    * 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

    * Shear-strain components in the three directions:  exy, eyz, ezx

    * Linear thermal expansion:                         et

    * Normal stresses in x, y, z directions             sx, sy, sz 

    * Shear stresses in xy, yz, zx planes               sxy, syz, szx 

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

    * Poisson's ratio:                                  P

    * Dilatation                                        D

    * Young's modulus                                   E

    * Lame's constant lamda = E*P/[(1+P)*(1-2*P)]       L                                         

    * Lame's constant mu = E/[2*(1+P)]                  G


      These definitions imply that the expressions

        L+G, 2*G and L+2*G  , which appear in equations 

      below, have the significances:

        L+G   =  E / [ 2 * (1+P) * (1-2P) ]

        2*G   =  E / (1 + P )

        L+2*G =  E * (1-P) / [ (1+P) * (1-2P) ]


    * Constant needed for expressing  thermal-expansion                                     

      effects        alpha * E / (1-2*P)                H

      where alpha is the longitudinal 

      thermal-expansion coefficient


    * Externally-applied forces in x, y and Z dirctns   X, Y, Z

    * The temperature excess above the initial state    T

      which entails that      et = alpha*T


    * Differentiation with respect to x, y or z is indicated by

      C:\WINDOWS\DESKTOP\ONLINE~1a next-line x, y or z


      so that du/dx is written as                       u


      likwise du/dy is written as                       u


      likewise du/dz is written as                      u


    * Thus: the dilatation, e, becomes            e = u + v + w 

                                                      x   y   z

    * the gradients    u  ,  v   and  w  , respectively,

                       x     y        z                           

      are equal also to:  ex ,  ey  and  ez by definition;

      therefore the dilatation can also be written as:

                                              e = ex + ey + ez              


    * the definitions of the shear-strain components (also called 

      'angular distortions') can likewise be written as:

        exy = u + v   ;  eyz = v + w   ;  ezx = u + w

              y   x            z   y            z   x


    * the definitions of the rotations and torques become:


        x-wise rotation                         i  = v  -  w

                                                     z     y


        y-wise rotation                         j  = w  -  u

                                                     x     z


        z-wise rotation                         k  = u  -  v

                                                     y     x



    *                                                                  *

    * Important note:                                                  *

    * Strictly speaking, i, j and k represent TWICE the rotations, as  *

    * properly derived by T&Goodier on page 233 and given              *

    * by them the symbols omega_x, omega_y and omega_z .               *

    *                                                                  *



        x-wise torque                           I  = Y  -  Z

                                                     z     y 


        y-wise torque                           J  = Z  -  X   

                                                     x     z


        z-wise torque                           K  = X  -  Y

                                                     y     x


        the torques being per unit volume    


    * div_grad is represented by a next-line ^, 

      so that     div_grad u becomes                    u




    2. The differential equations for displacements



    From T&G, pages 241 and 457,


       G*u  +  (L+G)*e  -  H*T   +    X    =  0          [2.1]

         ^           x       x 


       G*v  +  (L+G)*e  -  H*T   +    Y    =  0          [2.2]

         ^           y       y 


       G*w  +  (L+G)*e  -  H*T   +    Z    =  0          [2.3]

         ^           z       z 



    The corresponding normal stresses, from page 456 of the same 

    book, are:    


          sx =  L*e  +  2*G*u  -  H*T                    [2.4]



          sy =  L*e  +  2*G*v  -  H*T                    [2.5]



          sz =  L*e  +  2*G*w  -  H*T                    [2.6]



    The corresponding shear stress, from pages 7 and 10 of the same 

    book, are:


          sxy = G*exy = G*( u  +  v )                    [2.7]

                            y     x


          syz = G*eyz = G*( v  +  w )                    [2.8]

                            z     y


          szx = G*ezx = G*( w  +  u )                    [2.9]

                            x     z


    with temperature gradients making no appearance because they

    can affect rotations, e.g. u  -  v ,

                               y     x


    but can not afffect 'shear strains' (also called 'angular distortions') 

    such as:              e.g. u  +  v .

                               y     x 



    3. Boundary conditions for displacements



    3.1 Forces and stresses at solid-fluid surfaces



    At surfaces separating solid from fluids, the forces per unit area 

    in the cartesian directions x, y and z, namely X", Y" and Z", are 

    related to the stresses in the solid by the following equations, 

    given by T&G, p.236 :


      X" = sx * l + sxy * m + szx * n                    [3.1-1]


      Y" = sy * m + syz * n + sxy * l                    [3.1-2]


      Z" = sz * n + szx * l + syz * m                    [3.1-3]



    Here l, m and n are the direction cosines of the outward normal to 

    the surface, defined by:


      l = x                                              [3.1-4]



      m = y                                              [3.1-5]



      n = z                                              [3.1-6]



    where s stands for the distance along the normal.


    For surfaces aligned with the cell faces of a cartesian grid, the 

    values of the direction cosines are:


    *  for East faces:  l = 1 ; m = 0 ; n = 0     

    *  for North faces: l = 0 ; m = 1 ; n = 0     

    *  for High faces:  l = 0 ; m = 0 ; n = 1     



    3.2 Boundary conditions involving displacement gradients



    Substitution in 3.1-1, 3.1-2 and 3.1-3  of the expressions for 

    sx, sy, sz, sxy, syz, and szx of section 2 enables boundary 

    conditions for the displacements u, v and w to be written in terms 

    of their gradients and the external forces per unit area. This is 

    will be done for faces aligned with grid-cell surfaces, as follows:


    * For East cell faces (l=1, m=0, n=0):


      u = ( X" - L * e + H * T ) / (2 * G )              [3.2-1]



      v = - u + Y" / G                                   [3.2-2]

      x     y


      w = - u + Z" / G                                   [3.2-3]

      x     z


    * For North cell faces (l=0, m=1, n=0):


      v = ( Y" - L * e + H * T ) / (2 * G )              [3.2-4]



      u = - v + X" / G                                   [3.2-5]

      y     x


      w = - v + Z" / G                                   [3.2-6]

      y     z


    * For High cell faces (l=0, m=0, n=1):


      w = ( Z" - L * e + H * T ) / (2 * G )              [3.2-7]



      u = - w + X" / G                                   [3.2-8]

      z     x


      v = - w + Y" / G                                   [3.2-9]

      z     y



    Each of the gradients on the left-hand side of these equations has, 

    when multiplied by (L+G), the significance of a source in the 

    displacement-balance equation.


    3.3 Discussion



    It is worth pointing out to fluid-dynamicists that the 

    just-mentioned sources are partly similar to and partly different 

    from those which apply to velocities at solid-fluid boundaries. 

    Thus, if the equations for u are considered:


    (1) The X"s in each of the equations have their counterparts in the 

    corresponding equations for x-direction velocity; for x-directed 

    forces, whether applied normally or tangentially to the surface, 

    surely act as sources of x-directed momentum.


    (2) The H * T terms do not appear in the equations for velocity; 

    nevertheless, it is understandable that they appear in precisely the 

    same form for each displacement, because thermal expansion has no 

    preferred direction.


    (3) However, the L*e term represents a first big surprise for the 

    fluid-dynamicist; for the dilatation D contains all three 

    displacements; therefore gradients of v and w can provide sources of 

    u. Fluid-dynamics has not counterpart to this effect.


    (4) The surprise may be somewhat allayed by the recognition that the 

    material property L is zero when the Poisson's Ratio is zero; so the 

    L*e term can be regarded as a Poisson's-Ratio effect, which 

    is a peculiarity of solids from which fluids are mercifully free.


    (5) No such excuse alleviates the surprise occasioned by the 

    appearance of gradients v and y in the equations 3.2-5 and 3.2.8; 

    for they express an effect of the resistance to shear of solid 

    bodies which a fluid dynamicist might have thought as already 

    handled by the appearance of viscosity-like terms in the 



    (6) Finally, surprising or not, the unexpected terms MUST be 

    included if correct solutiions are to be obtained.



    4. The rotation equations



    4.1 Differential equations



    It will be shown to be useful, for the solution of solid-stress 

    problems by means of PHOENICS, to make use of differential equations 

    for the rotations i, j and k, defined in section 1.


    Differentiation of equation 2.1 with respect to y, and of equation 

    2.2 with respect to x, followed by subtraction, leads to:


              G*k  +   K   =  0.0                           [4.1-1]



    where k and K are respectively the rotation and the externally-

    applied torque defined above.


    Corresponding operations on the pairs of equations 2.2 with 2.3,

    and 2.3 with 2.1, lead to two further equations, namely:



              G*i  +   I   =  0.0                           [4.1-2]




              G*j  +   J   =  0.0                           [4.1-3]



    Here it is notable that the terms containing D and T make no 

    appearance, by reason of the identities:


          D  =  D     and    T  =  T                     .

          xy    yx           xy    yx  


    These equations do not appear in T&G; but it is hardly possible that 

    they have escaped attention by specialists in the theory of 




    4.2 Boundary conditions



    There being no guidance from T&G, the boundary conditions will be 

    derived from physics-based reasoning, attention being given to the 

    boundary conditions the k on north and east surfaces of a body of 

    which the cross-section in the x and y directions is indepedent of 

    the third dimension, z. 


    (1) As a first case, let the  body be a beam with uniform thickness in 

    the y direction much smaller than its length in the x-direction; and 

    let the only external forces be equal and opposite couples at the 

    low-x and high-x ends, only the latter being free to rotate.


    In these circumstances, the beam will bend in such a manner as to

    have uniform curvature; therefore rotation k will be linear in x, 


          k = k1 * x/x1

    where the length of the beam is x1 and the rotation at x=x1 equals 



    A linear k distribution satisfies the differential equation with 

    no sources at the north and south. Since there are no body forces 

    there either, the not unreasonable idea that "no forces means no

    sources" gains some support.


    It should be noted that the y-direction displacement is parabolic in 

    x, e.g. :

          v = a*x**2



    (2) Let us next consider the case in which the beam is loaded by a 

    tranverse y-direction force at the x=x1 end, which is known from 

    elementary beam theory to have as solution:


          v = a*x**2 + b*x**3


    This corresponds to a parabolic k distribution; and this satisfies 

    the differential equation for k only if there IS a finite source of 

    k at the north and south boundaries, on which however there are 

    still no external forces.


    This contradicts the "no force means so source" hypothesis; and it 

    presents a dilemma: by calculating the moment exerted by the force 

    at x=x1 we can deduce what is the bending moment a the x=x 

    cross-section; but this implies some "action at a distance" notion 

    which is not in keeping with the approach via differential 

    equations; and it gives no guidance as to what is to be done when 

    the body is not beam-like, i.e. long and thin.  


    A way out of the dilemma may be the following:

    - admit that there ARE force-like actions at the north and south 

    boundaries in BOTH cases;

    - recognise these as the "pseudo forces" which enter the 

    u-displacement equation via equation 3.2-5 above;



    ????????? to be continued ???????



    5. Differential equations for displacements in terms of rotations



    5.1 Mathematical manipulations



    Gradients of the dilatation D appear in the  equations 2.1, 2.2 

    and 2.3 . These can be written out in full as:


         D  =   u    +    v     +   w                       [5.1-1]    

         x      xx        xy        xz                           


         D  =   u    +    v     +   w                       [5.1-2]    

         y      xy        yy        yz                           


         D  =   u    +    v     +   w                       [5.1-3]   

         z      xz        yz        zz



    It is now useful to differentiate the definitions of the rotation 


         i with respect to x and y, 

         j with respect to y and z, and 

         k with respect to z and x, with the following results:


            k   =  u     -    v                         [5.1-4]    

            x      xy         xx


            k   =  u     -    v                         [5.1-5]    

            y      yy         xy 


            i   =  v     -    w                         [5.1-6]    

            z      zz         yz                            


            i   =  v     -    w                         [5.1-7]    

            y      yz         yy                            


            j   =  w     -    u                         [5.1-8]    

            x      xx         xz                            


            j   =  w     -    u                         [5.1-9]    

            z      xz         zz                             


    The following substitutions effect the desired simplifications:


    (a) into  5.1 from 5.5 and 5.9, to yield:     


         D  =   u    +    u     +   u   -  k   -  j          

         x      xx        yy        zz     y      z




         D  =   u  +  j   -  k                [5.1-10]

         x      ^     z      y



    (b) into  5.2  from  5.4  and  5.6, to yield:     


         D  =   v  +  k   -  i                [5.1-11] 

         y      ^     x      z



    (c) into  5.3  from  5.7  and  5.8, to yield:     


         D  =   w  +  i   -  j                [5.1-12] 

         z      ^     y      x




    5.2 The resulting new differential equations




    5.2.1 Equations in terms of rotation gradients


    Substitution for dilatation gradients from 5.10, 5.11 and 5.2 into

    the differential equations for the displacements yields three new 

    equations, namely:



     |                                                                |

     | (L+2*G)*u  +  (L+G)*(j - k)  -  H*T   +    X    =  0  [5.2-1]  |

     |         ^            z   y        x                            |

     |                                                                |

     | (L+2*G)*v  +  (L+G)*(k - i)  -  H*T   +    Y    =  0  [5.2-2]  |

     |         ^            x   z        y                            |

     |                                                                |

     | (L+2*G)*w  +  (L+G)*(i - j)  -  H*T   +    Z    =  0  [5.2-3]  |

     |         ^            y   x        z                            |

     |                                                                |



    These equations are easier to solve than these of section 2, from

    which they were derived, because:


    * the terms involving the rotation gradients can be evaluated 

      separately; and indeed they can often be represented by

      analytical expressions in terms of x, y and z; and


    * the linkages between u, v and w, formerly expressed by the 

      dilatation gradients, have disappeared.


    It is true that some linkages remain because D appears in the 

    boundary-conditions listed in section 3; but these are appreciably 




    5.2.2 Expressing the rotation-gradients via displacement gradients



    The definitions of j and k imply that:


      (j - k) =  w  - u    -    u  +  v

       z   y    zx   zz        yy    yx


    which can be written as:


      (j - k) =  w  - u    -    u  +  v

       z   y     xz   zz        yy    xy


    i.e. as:



      (j - k) = ez  - u    -    u  +  ey

       z   y     x    zz        yy     x


    i.e. as:                                                                   



      (j - k) = ez + ey - u  -  u         [5.2-4]                               

       z   y     x    x   zz    yy    


    Similar substitutions lead to:


      (k - i) = ex + ez - v  -  v         [5.2-5]

       x   z     y    y   xx    zz   




      (i - j) = ex + ey - w  -  w         [5.2-6]

       y   x     z    z   xx    yy                




    5.2.3 Introducing the shear stress



    The shear stresses sxy, syz and sxz are related to terms in equations

    2.2-4, 2.2-5 and 2.2-6 by the following relations:


          sxy/G = u  +  v                     [2.7]

                  y     x


          syz/G = v  +  w                     [2.8]

                  z     y


          sxz/G = w  +  u                     [2.9]

                  x     z


    With the aid of these equations the terms  u,  u,  v,  v,  w,  w

                                              yy  zz  xx  zz   yy xx


    can be eliminated, with the results:


    (j - k) = 2*(ez + ey) - sxy/G - sxz/G        [5.2-7]

     z   y        x    x     y       z       


    (k - i) = 2*(ex + ez) - sxy/G - syz/G        [5.2-8]

     x   z        y    y     x       z                 


    (i - j) = 2*(ex + ey) - sxz/g - syz/G        [5.2-9]

     y   x        z    z     x       y                 


    These equations are useful only where the shear stresses are known, 

    for example at solid-fluid boundaries. 




    5.2.4 The complete differential equation



    It will be seen that equation 5.2-4 contains the terms: - u  - u .

                                                             yy   xx

    Since those two terms appear with positive signs in the  u  term of


    the differential equation 5.2-1, it apppears appropriate it cancel 

    them out. The complete differential equation can now be written as:



    (L+2*G)*u + G*( u + u) + (L+G)*(ey + ez) - H*T + X  =  0  [5.2-10]   

           xx      yy  zz            x    x      x                   


    with corresponding equations for v and w, namely:       


    (L+2*G)*v + G*( v + v) + (L+G)*(ex + ez) - H*T + Y  =  0  [5.2-11]   

           yy      xx  zz            y    y      y                   


    (L+2*G)*w + G*( w + w) + (L+G)*(ey + ez) - H*T + Z  =  0  [5.2-12]   

           zz      xx  zz            z    z      z                   



    Let it now be observed that, the definitions of the relevant 

    quantities imply that H, defined above as: alpha * E / (1-2*P) ,

    can also be expressed as:  alpha * ( 3*L + 2*G).


    With alpha abbreviated to A, the roles of L and G in the equations 

    can be better clarified by re-writing them (omitting the *, for 

    brevity, except between A and T) as:


    (L+2G)u + G( u + u) + (L+G)(ey + ez) - (3L+2G)A*T + X  =  0  [5.2-13]

         xx     yy  zz           x    x             x                    


    (L+2G)v + G( v + v) + (L+G)(ex + ez) - (3L+2G)A*T + Y  =  0  [5.2-14]

         yy     xx  zz           y    y             y                    


    (L+2G)w + G( w + w) + (L+G)(ey + ez) - (3L+2G)A*T + Z  =  0  [5.2-15]

         zz     xx  zz           z    z             z                    




    5.2.5 Discussion



    It must now be admitted that these equations could have been derived 

    directly from those section 2, and that the introduction of rotation 

    has proved to be a needless digression.


    Thus, the first term of 5.2-13 can be split into two, thus:


           (L+2G)u   = (L+G)u + Gu

                xx         xx   xx


    Therefore the first three terms can be written as:


           (L+G)u + Gu + G( u + u) + (L+G)(ey + ez)            

               xx   xx     yy  zz           x    x 

    i.e. as:

               G( u +  u + u) + (L+G)( u + ey + ez) 

                 xx   yy  zz          xx    x    x  

    i.e. as:

               G u  + (L+G)e

                 ^         x


    These are precisely the first two terms of equation 2.1 .


    Which form of equation is likely to prove most convenient for 

    numerical solution? Probably that of section 5.2.4 because, in a 

    segregated-variables solution procedure, it is possible to express 

    all the terms containing dependent variable (u, v or w) by way of 

    the finite-volume 'molecule', so that only those containing the 

    strain gradients need to be inserted explicitly as source terms.


    Let it however be remembered that the introduction of the strain 

    gradients into the above expresions is not an obligatory move; for

      ex  can also be expressed as   u   .

       y                            xy


     In the following section it will appear that the second of these 

     has some advantages.  


    5.2.6 Finite-volume expression



    Whatever the form of the differential equation and boundary 

    conditions, it is how they are expressed in finite-voume (i.e. FV)

    form which is decisive for the satisfactory working of the computer 

    program. Attention will now be given to the FV equation for the 

    displacement u, the grid being cartesian.


    The following diagram represents the grid in two dimensions; but the 

    existence of the z (High-Low) direction must also be borne in mind.


    !         !         !         !         !         !

    !         !         !    N    !         !         !

    !         !         !   -->   !         !         !

    !   ^     !         !/////////!         !         !

    !   !     !         !/////////!         !         !


    !   !     !        NW/////////NE        !   !     !

    !   !     !    W    !/// P ///!    E    !   !     !

    !   v     !   -->   !   -->   !   -->   !  dyp    !

    !   ^     !         !\\\\\\\\\!         !   !     !

    !   !     !         !\\\\\\\\\!         !   !     !


    !   !     !        SW\\\\\\\\\SE        !         !

    !   !     !         !\\\ S \\\!         !         !

    !   v     !         !   -->   !         !         !

    !         !         !<--dxp-->!         !         !

    !         !     <--dxw--><-- dxe-->     !         !


    Consider the terms of equation 5.2-1 which concern displacements,


              (L+2*G)*u  +  (L+G)*(j - k)

                      ^            z   y 

    Integration of the first term over the control volume surrounding P

    leads ( with obvious notation) to:


         (L+2G)*{  dyp*dzp*[( uE -uP) / dxe  - ( uP - uW) / dxw]

                 + dxp*dzp*[( uN -uP) / dyn  - ( uP - uS) / dys]

                 + dxp*dyp*[( uH -uP) / dyh  - ( uP - uL) / dzl] }



    This is a typical diffusion-type expression, requiring no comment.


    Consider now the second part of the second term, involving the 

    gradient of the z-direction rotation component, k, the integration 

    of which over the control volume surrounding P yields:


         - (L+G)*dxp*dzp*(kn - ks) 


    where kn and ks are the values of k in the control volumes shaded

    ////// and \\\\\\ respectively.   



    These quantities can be represented in discretized form as:


         kn = ( uN - uP ) / dyn  - ( vNE - vNW) / dxp , and


         ks = ( uP - uS ) / dys  - ( vSE - vSW) / dxp 


    Corresponding expressions can be derived for the j-gradient term, 

    involving uH, uL, wHE, wHW, wLE, wLW, dzh and dzl. However the 

    necessary conclusions can be drawn by confining attention to the 



    All the above terms involve displacement differences. Their combined 

    influence can therefore be represented as:


        ( uE - uP )  * [ (L+2G)* dyp*dzp / dxe ]

      + ( uW - uP )  * [ (L+2G)* dyp*dzp / dxw ]

      + ( uN - uP )  * [  G* dxp*dzp / dyn ]  

      + ( uS - uP )  * [  G dyp*dzp / dys ]  

      + ( vNE - vNW) * [ (L+G) * dzp ] 

      - ( vSE - vSW) * [ (L+G) * dzp ]


    with corresponding z-direction terms.


    Consider now the situation in which the point P is located at the

    solid-fluid interface, as shown below:




           !         !        

           !         !    N   

           !         !   -->  

     ^     !         !/////.  

     !     !         !/////.  

    dyn----!----k----^----kn    ^

     !    NNW       NW/////.    !

     !     !    W    !/// P.    !

     v     !   -->   !   -->   dyp

     ^     !         !\\\\\.    ! 

     !     !         !\\\\\.    ! 

    dys----!----k----^----ks-   v-

     !    SSW       SW\\\\\.   

     !     !         !\\\ S.   

     v     !         !   -->   

           !         !<--dxp-->!       


     The rotations kn and ks cannot be computed as before, because 

     there is no direct way of calculating dv/dx. This quantity could

     be computed indirectly from the shear-stress boundary condition 2.7



          sxy = G*exy = G*( u  +  v )                    [2.7]

                            y     x

          kn =  2*dy/dy - sxy  




          kn =  2*(uN - uP) / dyn  - sxyn/G


          ks =  2*(uP - uS) / dys  - sxys/G


    However, there does not appear to be any advantage in so doing.


    Instead, it appears preferable to presume that dv/dx varies little

    with x, so that its value can be extrapolated from the next cells to 

    the west. The variation of dv/dx with y is what is important.



    5.3 Cylindrical polar coordinates



    5.3.1 Equations for displacements in terms of polar co-ordinates



    According to B&W, p. 257, the equations for the displacements in 

    terms of gradients of rotation can be written in polar 

    co-ordinates as follows:



    (L+2G)D - (L+G)(j - k)  - (3L+2G)A*T + X  =  0             [5.3-1]

         rx         z   y            rx                    


    (L+2G)D - (L+G)( k - i) - (3L+2G)A*T + Y  =  0             [5.2-2]

          y         rx   z             y                    


    (L+2G)D - (L+G)(ri - j) - (3L+2G)A*T + Z  =  0             [5.2-3]

          z         ry  rx             z                    




    * the dilatation, D =  rv  +   u  +  w                     [5.2-4] 

                           ry     rx     z                          

    * the rotation about the z-axis, k = ru  -  v              [5.2-5]

                                         ry    rx

    * the rotation about the x-axis, i = v  -  w               [5.2-6]

                                         z     y

    * the rotation about the y-axis, j =  w  -  y              [5.2-7]

                                         rx     y

    and, x is now measured in radians;


         r stands for radial distance, which is identical to y,

           so that dr=dy; 


         where r appears it has the significance of a multiplier

         whereas the appearance of y on the second line signifies

         differentiation, so that:  


         rv  stands for [d(rv)/dy]/r                              



          u  stands for [d(u)/dx]/r



         ru  stands for [d(ru)/dy]/r   



          v  stands for [d(v)/dx]/r



         and so on.



    5.4 One- and two-dimensional situations                          



    5.4.0 Definitions and preliminary considerations                 



    (a) Definitions and examples                                     



    It frequently arises that the shape of the solid in question, and 

    the directions and distributions of displacement-producing 

    influences, whether mechanical or thermal, result in distributions 

    of displacement and temperature which vary in one direction only.


    [ For brevity, the words 'width', 'thickness' and length will be 

    associated with the cartesian coordinates x, y and z respectively.


    In axi-symmetrical situations, for which polar coordinates may be

    used, the radial direction will always be y; and x will represent 

    angle around the z-axis. ]


    Examples of one-dimensional distributions are:


    A. A long wire of uniform diameter is pulled by equal and opposite 

    forces, one at each end. Its longitudinal-direction strain and 

    stress are uniform along its length; but the longitudinal 

    displacement varies linearly with longitudinal position (z).


    B. A uniform-sectioned block is compressed by equal and opposite 

    forces applied uniformly to high-x and low-x faces. All 

    displacements, stresses and strains depend on the x coordinate alone; 

    and this is true whether its thickness and/or length are large or 

    small compared with its width.


    C. A similar block is placed suddenly into contact at its high-x and 

    low-x faces with objects of uniform but unequal temperature. The 

    temperature of the material within the block changes with time, as a 

    consequence of heat conduction; but the temperature at any particular 

    time depends only on x.



    Examples of two-dimensional distributions are:


    d. The forces on the high-x and low-x faces of a block such as that 

    in example B. vary with its position in the width dimension, y. Then 

    the stresses and displacements vary with both x and y, but not with 

    the third (length) coordinate, z.


    e. In examples B. and C., the width varies with x, but not the 

    length. Then the stresses, displacements and temperatures vary 

    with both x and y, but not with  z.


    f. The block of example A. is clamped firmly at the low-x face and 

    subjected at the other to a force which is uniform over that face 

    but which is directed in the y direction. The block therefore acts 

    as a cantilever beam, subject to bending stresses.



    (b) A difference between the solid-stress and fluid-flow equations



    Fluid-flow specialists are familiar with the fact that, although all 

    real flow  phenomena are three-dimensional, much can be gained from 

    the study of one- and two-dimensional idealisations; however they 

    enjoy an advantage which the solid-stress specialist does not, 



    * the dimensions which are NOT the 'one', or one of the 'two', need 



    Thus, if a fluid flows in such a manner that its properties vary 

    with x and y but not with z, the terms which are connected with z, 

    for example:   u      or     u    , are all zero.

                   z             zz


    Where displacements, strains and stresses in solids are concerned, 

    however, this is not so. Thus, as will be explained below (see 

    section 5.4. ), u    may be zero without  u   being zero also.

                    v                         vv


    Why is this? In part it is due to the phenomenon expressed by the 

    non-zero value of Poisson's ratio, namely: stresses in the 

    x-direction normally cause either strains or stresses in the y- 

    and z-directions as well as strains in the y-direction. 


    However, this is not the whole of the matter; for inter-connexions 

    between the three directions are implied by the fact that gradients 

    of the dilatation D, which involve gradients of u, v and w, appear 

    in the differential equations for equations each one of these 



    It seems therefore wisest to warn the CFD specialist that, although 

    the equations governing displacements are superficially similar to 

    those governing velocities, there exist subtle differences which 

    those seeking to solve them must respect.



    (c) Consequences for the equation-solving procedures



    (i) The solid-stress-analysis literature pays little attention to 

    one-dimensional situations; and, when it does so, it is usually 

    postulated that the material is unconstrained in the other two 

    directions, as in the wire under tension of case A. above.


    (ii) By contrast, special nomenclature is employed for two-dimensional

    situations, attention being focussed on two extremes, namely, if the 

    x-y plane is that in which variations are to be determined:


    * In extreme (1), the material is free to expand or contract in the 

    z-direction. This is called the 'PLANE-STRESS' situation, because 

    direct stresses in the z-direction are negligible. 


    It arises, typically, when the z-direction dimensions of the solid 

    in question are small compared with those in the other two 

    directions. However, this is far from being sufficient; for it is 

    possible that such a body, however thin, could be confined between 

    two totally inflexible planes, which prevented their z-direction 



    For the 'plane-stress' condition to prevail, such confining surfaces 

    must NOT exist.


    * In extreme (2), the material is NOT free to expand or contract in 

    the z-direction, being constrained, it is supposed by rigid surfaces 

    to which it must adhere, while still being free to slide past it in 

    the x or y directions. 


    This is called the 'PLANE-STRAIN' situation, because strains in the 

    z-direction are negligible.


    This situation arises when inflexible confining surfaces of the 

    just-discussed kind DO exist; but it arises also when the 

    z-direction dimensions of the solid are large compared with those in 

    the x- and y-directions; for then, if significant variations with z 

    did exist, they would be opposed and nullified with the 

    consequential shear stresses on planes of constant x and y.


    (iii) For one-dimensional situations, the fact that the solid-stress 

    literature appears to have neglected them does not mean that they 

    SHOULD be ignored. It appears therefore appropriate to define two 

    additional situations, namely:


    (3) LINEAR STRESS, in which (if z is the longitudinal direction) 

    stresses in the x- and y-directions are neglible. The 'wire' of 

    example A. illustrates this.


    (4) LINEAR STRAIN, in which the solid is enclosed in a totally rigid 

    casing which prevents its boundaries from moving in (say) the y- 

    and z-directions.



    (d) Particular consequences for the numerical analyst



    The fluid-dynamicist, when dealing with a one-dimensional situation, 

    will discretize space in the appropriate (say, x) direction; and the 

    only solved-for velocity will then be u. Velocities in the other 

    directions, namely v and w, will make no appearance.


    Similarly, if the problem is two-dimensional (with variations in x 

    and y, say), the u and v velocities will be solved for; but not the 

    z-direction velocity, w.


    If a computer code designed for computing velocities is to be employed 

    (as PHOENICS is) for computing displacements, the differences beteen 

    the equation systems, mentioned in (b) above, must be addressed; and 

    actions may be needed in three areas, namely:


    (1) the calculation of the dilatation;

    (2) the calculation of the forces at boundaries; and

    (3) (when thermal effects are significant) by provision of sources 

        of displacement in the solved-for directions.








    5.4.1 The displacement gradients at constrained boundaries



    Let it be supposed that a solid rectangular block is subjected to 

    uniform pressure on its low-x and high-x faces. The problem might be 

    regarded as being one-dimensional, with only equation 5.13 requiring 

    to be solved. 


    This is true; but the y and z directions cannot be completely 

    disregarded for the reason that the x-direction behaviour depends on 

    whether the low-y and high-y, and/or the low-z and high-z, faces are 

    held fixed, allowed to move freely, or partially constrained.


    The dependence of the x-direction displacements on the y- and 

    z-direction constraints is expressed by way of the dilatation, i.e.

    the sum of the strains ex, ey and ez; therefore ey and ez must be 

    given values.


    Two extreme cases are common, namely:

    (a) the strains are zero because the surfaces are held fixed; and

    (b) the stresses are zero, because the surfaces are unconstrained.


    However, although neglected by most textbooks, a third situation may 

    arise, namely: 

    (c) that of linear constraint, whereby the displacement may be 

    proportional to the force exerted at the boundary.


    Since case(c) contains the other two as special cases, with 

    proportionality constant zero for case (a) and infinite for case(b), 

    it is this which will be considered.


    Let the forces exerted on the high-x, high-y and high-z boundaries 

    be respectively:    - bx*ex , - by*ey  and  - bz*ez , where ex, ey 

    and ez are the displacements at the constrained boundaries. Then it 

    follows from equations 3.1, 3.2 and 3.3, that:


          ex =  ( -bx*ex -  L*e  +   H*T ) / (2*G) 


          ex =  ( - L*e  + H*T) / ( 2*G + bx )          [5.16]



          ey =  ( -by*ey -  L*e  +   H*T ) / (2*G)


          ez =  ( - L*e  + H*T) / ( 2*G + by )          [5.17]



          ez =  ( -bz*ez -  L*e  +   H*T ) / (2*G)


          ez =  ( - L*e  + H*T) / ( 2*G + bz )          [5.18]



    These are the equations to be used when (and ONLY when) the strain 

    in question has NOT been deduced from the differential equation for 

    the corresponding displacement.    


    5.4.2 Gradients of ex, ey and ez at such boundaries



    (a) The problem


    Equation 5.13 contains the term u , i.e. u   +   u   +   u   ;

                                    ^        xx      yy      zz

    and, in a one-dimensional situation, with x as the direction in

    which u is allowed to vary, it might be supposed that the 

    u and u                                

    yy    zz

    terms might be neglected. This is not the case however, at 

    least when thermal expansion is present.


    Consider, for example, the case in which the temperature T varies 

    with x, and the proportional linear thermal expansion, et;

    the proportionality coefficient is the material constant, alpha. 


    Let it also be supposed that the y- and z-direction surfaces are not 

    constrained in any way, so that the three strains at any point obey

    the equality relation:


              ex  =  ey  =  ez  =  et                         [5.19]


    Then initially-plane (constant-x) surfaces will become spheres of 

    radius equal to:   alpha * T  .



    Of course, if movement in x is allowed, but not in z, the 

    initially-plane surfaces will become cylinders of the same radius 

    about an axis normal to the xy plane; and, if movement in z is 

    allowed, but not in x, the initially-plane surfaces will become 

    cylinders of the same radius about an axis normal to the xz plane.


    Finally, if movement is allowed in neither the y nor z directions, 

    plane surfaces will remain plane.                          


    (b) The distribution of u


    Consider the simple case in which  et  is linear in x, for example:


      et =  a * ( - 1 + 2*x/l)                                [5.20]


    which corresponds to the case of a block of thickness l, which is 

    cooled to a lower-than-average temperature where x=0 and heated to a 

    correspondingly higher-than-average temperature at the opposite 

    face, the linear temperature being brought abour, in the steady 

    state, by heat conduction.


    In this case, alpha * T   = 2*a/l                         [5.21]


    The uncontrained nature of the process entails that the linear

    extensions in each direction, ex, ey and ez are all exactly equal 

    to et, from which the following expressions concerning first and 

    second gradients of displacement and rotation can be derived:-



    (i) u  :  Differentiation of  u  = ex (by definition) = et

        xx                        x

     leads, via 5.20,  to:   u     =  2*a/l                   [5.22] 



    (ii) u  :  Integration of  v  = ey (by definition) = et

         yy                    y


     leads, via 5.20,  to:   v  =  a*(-1+2*x/l)*y + const     [5.23] 


     whereafter differentiation with respect to x leads to:

                             v  =  2*a*y/l                    [5.24]


     This relation is useful because, since there can be no shear 

     stresses, equation 2.7 leads to:

                u  +  v   =  0,

                y     x

     from which it follows, via differentiation of 5.24, that:

                u   =  -  v     =  - 2*a/l                    [5.25]

                yy        xy   


    (iii) u  :  A similar manipulation of the corresponding terms

          zz    involving w and z lead to:


                u   =  - 2*a/l                                [5.26]



    (c) Consequences for one- and two-dimensional analyses



    It follows that the supposition mooted at the start of section 

    5.4.2, namely that the u   and   u   tems might be neglected is

                           yy        zz

    far from the truth. If they were neglected, u would equal 2*a/l;                      


    but, when they are taken into account, u equals - 2*a/l !



    (d) The rotation-gradient terms



    Equation 5.13 also contains the term:  (L+G)*(j - k)  .

                                                  z   y  

    In one- or two-dimensional analyses one might be tempted to set

    z- and y-direction gradients to zero here also; but that we be 

    equally wrong, as will now be explained.


    The term j , by reason of the definition of j, can be expressed as:


       j  =   w    -   u                                    

       z      xz       zz                              


    Since w = a*(-1 + 2*x/l)z + constant, derived in a simlar manner to 

    5.23, differentiation leads to:

       w  =   2*a/l   ; and so, by reason of 5.26, we deduce that



       j  =   4*a/l                                       [5.27]


    A similar argument shows that k  has the same value; and the upshot 


    is that, to account for the effect of the unconsidered y and z 

    directions, it is necessary to add to the terms in the equation:


     -(L+2*G)*(4*a/l) to account for the second differential 

                       coefficients in the first term, and

      (L+G)*(8*a/l)    to account for the rotation gradients, wherein


    it will be recalled that   2*a/l = alpha * T/x  .      




    6. Implementation in PHOENICS


    6.1 Comparison of the differential equations



    PHOENICS treats the displacements as though they were velocities; 

    therefore what is to be done can be deduced by comparing equation 

    5.13, for example, with the corresponding equation for the 

    x-direction velocity. 


    This is, for the steady state, with uniform effective viscosity emu, 

    and with the convection terms omitted:



      emu*u  +  p  + source_terms  + X  = 0                   [6.1]

          ^     x


    Comparison of this equation with 5.13 dictates that:


    (a) the emu store should contain  L+2*G ;


    (b) the pressure-gradient term should be omitted;


    (c) (L+G)*(k - i )  should appear as a source term; 

               z   y     


    (d) so also should    -  H*T .



    Similar conclusions can be drawn about the equations for the y and x 



    The equations are of course solved in finite-volume form, as far as 

    possible in the same way for both displacements and velocities. How 

    this is done will now be discussed.                          



    6.2 The material-dependent properties



    (a) The PHOENICS PROPS file


    PHOENICS provides information about several properties of fluid and

    solid materials by way of the so-called PROPS file, which is arranged

    in columns, thus:

  IMAT density   viscosity  specific  thermal thermal  compress

  i.e.           or Poisson   heat     conduct expnsn.  -ibility or

  PRPS          's ratio    capacity  -ivity  coefff.  1/Young's Modulus

    The 'or' in columns 3 and 6 refers to whether the material is fluid     

    or solid. Thus, in the case of solid aluminium, where the entry is:


     <100>   ALUMINIUM at 27 deg c

   100    2700.0     0.35      896.0     204.0 2.35e-5  1.47e-11


    the 0.35 is the  Poisson's ratio and the 1.47e-11 is the reciprocal

    of the Young's Modulus.

    (b) The three-dimensionally-stored quantities



    When so instructed, PHOENICS creates three-dimensional storage for 

    each of the above six properties; therefore, unless other 

    arrangements are made, it will place the viscosity or the Poisson's 

    ratio in the corresponding storage location according to whether 


    However, within the solid it is COMPOSITE properties which are 

    needed, as expressed by the terms L+G, L+2G, H, etc which appear in 

    the differential equations and boundary conditions. Therefore, in 

    order to simplify the coding when PHOENICS reads the above-described

    PROPS ile, it places, for solid-containing cells:


         L+2*G  in the viscosity store,

           L+G  in the compressibility store, and

             H  in the thermal-expansion coefficient store,


    these being computed, in accordance with definitions, from:         


       L+2*G  =  E * (1-P) / [ (1+P) * ( 1-2*P)]             [6.2]


         L+G  =  E / [ 2 * ( 1+P ) * (1-2*P) ]               [6.3]


           H  =  alpha * E / (1-2*P)                         [6.4]


    where alpha is the thermal-expansion coefficient.


    These three composite properties are those appearing in the 

    differential equation of section 5.2. In the boundary conditions of 

    section 3, the L and G which are needed on their own are easily 

    deduced from:


          G = (L+2*G) - (L+G)                                [6.5]

          L = 2*(L+G) - (L+s*G)                              [6.6]


    When desired, the Poisson's ratio, P, and the thermal-expansion

    coefficient, alpha, can be recovered from the stored quantities via:


          P     = 1.0 - 0.5 * (L+2G) / (L+G)                 [6.7]

          alpha = H / [4.0*(L+G) - (L+2G) ]                  [6.8]


    (c) A device to prevent round-off error. 



    As is exemplified by the above data for aluminium, the reciprocal of 

    the Young's Modulus is usually of the order of 1.e-11. For this 

    reason, PHOENICS internally multiplies it by 1.e11, in order that 

    excessively small coefficients shall not be generated , so inducing 

    round-off error. Applied forces are correspondingly divided by 1.e11 

    so as to ensure that displacements, stresses and strains are 

    correctly reported.      


    6.3 The dilatation



    Although the dilatation, D, does not appear in the differential 

    equations which are solved, it does so in the boundary conditions. 

    Therefore D is calculated within the solver; and its values are 

    stored in the locations which, for fluid regions, are filled by 

    values of the pressure.


    For reduced-dimensionality situations, of the kind discussed 

    in section 5.2 above, the following practice is adopted:


    * If NX = 1, so that no equation for u is solved, the strain ex is

      computed from equation 5.16, the constant bx being supplied by way 

      of a spedat statement in the q1 file thus:

         SPEDAT(BOUNDARY,XCONST,R,value of constant bx)


    * If NY = 1, so that no equation for u is solved, the strain ex is

      computed from equation 5.17, the constant bx being supplied by way 

      of a spedat statement in the q1 file thus:

         SPEDAT(BOUNDARY,YCONST,R,value of constant by)


    * If NX = 1, so that no equation for u is solved, the strain ex is

      computed from equation 5.18, the constant bx being supplied by way 

      of a spedat statement in the q1 file thus:

         SPEDAT(BOUNDARY,ZCONST,R,value of constant bz)


    Often the values of the constants are either:

    * 1.e20, i.e. a large number allowing no displacement, 


    * 0.0, implying no constraint whatever.    


    6.4 Applied-force boundary conditions



    As shown in section 3, displacement-gradient boundary conditions 

    depend on D, T and G, as well as on the applied force. However, if 

    a force is defined by way of a patch of FIXFLU type, having FORC as 

    the first four characters of its name, the fourth argument of the 

    corresponding COVAL can be the force itself; for then the D, T and 

    G influences will be introduced automatically.


    Forces exerted by the fluid on the solid are represented 

    internally, i.e. without user intervention.



    6.6 Rotation-gradient sources



     PHOENICS provides, when necessary, cell-by-cell storage space for

     three rotation-gradient sources, namely:


          RGRX  =  (k - i )                                  [6.9]

                    z   y  


          RGRY  =  (i - j )                                  [6.10]

                    x   z  


          RGRZ  =  (j - k )                                  [6.11]

                    y   x  


      Their cell-by-cell distribution may often be calculated, once for 

      all, at the start of the flow simulation. Alternately they may be 

      computed  by solving the differential equations for i, j and k.




    7. Special cases


    7.1 Uniform direct stress


    a. The problem



    A rectangular parallelepiped, of dimensions lx, ly and lz, is held 

    at its low-x, low-y and low-z faces and subjected to an applied 

    tensile force of X" per unit area at its high-x face.


    Three cases are considered:

    (i) both high-y and high-z faces are free to move

    (ii) both high-y and high-z faces are prevented from moving

    (iii) the high-y face is free to move but the high-z face is not.


    Case(i): motion allowed in y and z directions.


    * Application of Hooke's law, and the definitions of Young's modulus

      E and the Poisson's ratio P lead directly to:

        ex = X"/E                                     [7.1]

        ey = P*X"/E                                   [7.2]

        ez = P*X"/E                                   [7.3]

      from which the dilatation may be deduced as

        D  = (X"/E) * ( 1 - 2P )                      [7.4]


    * As a consequence of 7.1, the maximum displacement will equal

      ex * lx, i.e. (X"/E)*lx.                        [7.1a]


    * It is a useful check on their derivations to ensure that the above 

    differential and boundary-condition equations are satisfied by 

    equations 7.1, 2, 3 and 4.


    The differential equations 5.1, 2 and 3 reduce to

        u  =  0 ,   v = 0   and   w  =  0  

        ^           ^             ^

    because there are no sources from rotation gradients, temperature 

    gradients or applied forces within the body.


    These equations are all satisfied by uniformity of ex, ey and ez, 

    which are the gradients of u, v and w, expressed by 7.1, 2  and 3.   


    The boundary conditions reduce, in the absence of temperature 

    effects and of Y" and Z", to


      ex =  ( X" -  L*e ) / (2*G)            [3.1a]



      ey =  (    -  L*e ) / (2*G)            [3.2a]



      ez =  (    -  L*e ) / (2*G)            [3.3a]


    Summation leads to:


       D = ( X" - 3*L*e ) / (2*G)


    and so to 


       D * ( 2*G + 3*L) = X" 


    The definitions of G and L imply that:     

       2*G + 3*L  =  E / ( 1 - 2P )



       D * E / ( 1-2P)  = X"  


    which is consistent with 7.4 .


    Therefore both the differential and boundary-condition equations 

    have passed their tests.   



    Case (ii): motion disallowed in y and z directions.


    * The solution cannot for this case be obtained by direct 

      application of Hooke's law, and the definitions of Young's 

      modulus and the Poisson's ratio; but simultaneous solution of

      the general equation leads to the following results:


    ** The differential equations 5.1, 2 and 3 dictate constancy of 

       u, v and w  gradients because of the absence of source terms;

       therefore ex. ey and ez are all constants.


    ** ey and ez are both zero, by reason of the prevention of movement 

       in the y and z directions. 


    ** D and ex are therefore identical, because D = ex + ey + ez  by



    ** Boundary condition 3.1 thus reduces to:


          ex = ( X" -  L*ex ) / (2*G)            [3.1b]

       with the consequences:

          ex =  X" / ( L + 2*G )

       and so, from the definitions of L and G, 


          ex = ( X"/E )*(1+P)*(1-2P)/(1-P)   [7.5]


    ** Boundary condition 3.2 enables the y-direction stress to be

       deduced, thus:


         ey =  0  = ( Y" -  L*e ) / (2*G)            [3.2b]




         Y" = L * D


       and so, since D equals ex, which is given by equation 7.5 ,


         Y" = X" * P / (1-P)                   [7.6]


     ** Similarly, boundary condition 3.3 leads to:


         Z" = X" * P / (1-P)                   [7.7]



    Case (iii): motion disallowed in y but not in z directions.


    ** this case is characterised by ey=0 and sz=0. Therefore 3.3a and 

       3.2b are valid.


    ** Then the dilatation is given by:


         D = ex + ey + ez

           = ex      - L*e/2*G


       so that 


         D = ex * 2G / (L + 2*G) 



    ** Substitution in 3.1 now leads to:


         ex = [ X" - L * ex * 2*G / ( L + 2*G)] / 2*G


       and so to


         ex = (X"/E) * (1+P) * (1-P)            [7.8]


    ** Finally, 


         sy = L * D

            = L * (X"/E) * (1+P) * (1-P)

            = (X"/E) * P / (1-2*P)              [7.9]



    7.2 Uniform temperature


    a. The problem


    A rectangular parallelepiped is held at its low-x, low-y and low-z 

    faces and is heated to a temperature T above its reference value.


    Four cases are considered:

    (i) the high-x, high-y and high-z faces are free to move;

    (ii) the high-x face is fixed but the high-y and high-z faces 

        remain free;

    (iii) the high-x and high-y faces are fixed but the high-z 

        remains free;

    (ix)  all faces are fixed;



    Case(i): motion allowed in all three directions


    Since there are no stresses in any direction, addition of equations 

    3.1, 3.2 and 3.3, together with the substitution of D for ex+ey+ez,

    leads to:

         ( 2*G + 3*L )*e = H*T

    which, when G, L and H are expressed via their definitions in terms

    of E, P and alpha in the nomenclature, reduces to:


         D = 3 * alpha*T                 [7.10]

    and so

         ex = ey = ez = et = alpha*t     [7.11]


    Neither poisson's ratio P, nor Young's modulus E, plays any part in 

    these formulae, which correspond to stress-free isotropic thermal 



    Case(ii): motion prevented only in the x-direction


    In this case, ex must be zero; therefore the relationship between the 

    dilatation D and the thermal strain et can be derived by adding 

    equations 3.1 and 3.2, and setting ey+ez equal to D.


    The result, after substitution from the definitional relationships, 



         D =        alpha*T * 2 * (1 + P)     [7.12]

         ey = ez =  alpha*T * (1 + P)         [7.13]


    Understandably, when P is zero, the dilatation is alpha*T times 2,

    rather than times 3 as in case (i), because the material can expand 

    in only two directions. With finite P however, the x-direction 

    compressive stress causes additional strains in the y- and 

    z-directions; so the dilatation becomes larger.


    The magnitude of the x-direction stress is obtained by substitution 

    from 7.12, and the definitional relationships, into 3.1. The result 


        sx = - alpha*E*T                      [7.14]



    Case(iii): motion prevented in both the x- and y-directions


    When ez is the only non-zero strain, the dilatation is deduced by 

    setting D=ez in equation 3.3, with the result, after the usual 



        D = ez = alpha*T * ( 1+P ) / ( 1-P )  [7.15]


    Once again, the dilatation increases with increase in Poisson's ratio, 

    and indeed twice as rapidly as before, because there are two 

    compressive stresses (i.e. those in the x- and y-directions) 

    contributing to z-direction strain.    


    The corresponding x- and y-direction stresses are:


        sx = - (alpha*E)*T / ( 1-P )          [7.16]


    Case(iv): Motion prevented in all three directions.


    In this case, ex, ey, ez and therefore D also, are all zero; and the 

    stresses in the three directions are directly deducible from 3.1, 3.2 

    or 3.3 as:


        sx = sy = sz = - H*T = - (alpha*E) * T / ( 1-2*P)  [7.17]


    This shows that the stress increases with Poisson's ratio, 

    corresponding to the fact that compression in one direction tends to 

    cause expansion in a direction at right angles which, if it is to be 

    prevented, requires a further compressive stress.    


    7.3 Non-uniform direct stress



    Let the force applied to the high-x face of the block in section 7.1 

    be removed; but let the same total force be applied uniformly through 

    the body, as though by a gravitational body force. We may then 

    expect the displacement at the end to be exactly one half of what 

    was before.


    Let the value of this force, per unit volume, be X, which will be 

    equal to the previous X" divided by lx.


    Once again, three cases will be considered, namely:

    (i) strain is permitted in the y- and z-directions;

    (ii) strain is not permitted is not permitted in either y- or 


    (iii) strain is permitted in the y- but not the z-direction .


    In all cases, the direct stress per unit area which is needed to 

    balance the force is: X*(lx - x), having its maximum value at the 

    fixeb end where x=0 and its minimum at the free end where x=lx . 


    Case (i) Y- and z-direction movement permitted.


    The distribution of x-direction displacement will be deducible by 

    integrating the differential equation:


       ex = u =  (X/E)*( lx - x )


    with the result:


            u = (X/E)*( lx*x - 0.5*x**2 )                     [7.18]


    As expected, the maximum displacement, namely 0.5*(X/E)*lx**2 where 

    x-lx, is exactly half that given by equation 7.1a (with X"= X*lx).     


    The y- and z-direction strains will be - P*ex, so that the 

    dilatation will be:


            D = ( 1.0 - P ) * (X/E)*( lx - x )                [7.19]


    which has the same value at x=0 as in section 7.1, but which 

    diminishes to zero at the free end (x=lx))        


    The stresses in the y- and z-directions will of course be zero. 


    Case(ii) No movement permitted in y- or z-directions



    The analysis of this case follows the same pattern was established 

    for case (ii) of section 7.1; indeed all the equations there 

    continue to be valid if X" is replaced by X*(lx - x).


    The stresses, strains and dilatation of course now vary linearly 

    with x.


    Case(iii) No movement permitted in y- or z-directions



    What has just been written about case (ii) applies to case(iii) also.



    7.4 Non-uniform temperature


    (a) Linear temperature profile in an unconstrained block



    (1) The problem



    * Let it now be supposed that the low-x face of the rectangular block 

    is held at the temperature -T0 and the high-x face at temperature 

    +T0. If the thermal conductivity is uniform and there are no other 

    heat sources, the steady-state temperature distribution will obey:


       T = T0 * ( -1 + 2x/l ) ,                             [7.20]


    where l is the x-direction dimension of the block


    * Let the block be free to expand in all three directions, so that 

    the corresponding stresses are all zero.


    It is easy to see that the local strains are represented by:


       ex = ey = ez = et = e/3 = alpha*T;                   [7.21]


    so they are also linear in x; and the x-direction displacements u 

    can be derived by integration of the definitional relation:


       u  =  alpha*T0*( -1 + 2x/l)                          [7.22]



    with the result:


       u  =  alpha*T0 * ( -x + x**2/l )                     [7.23]


    and the consequence that u has a minimum = (1/4)*alpha*T0*l at 

    x=l/2 and the value zero at the high-x face.


    * The block does not change in overall length because the shrinkage 

    near its colder surface is exactly compensated by its expansion near 

    the hotter surface.


    * The only material property which influences the solution is alpha; 

    for, there being no stress, neither Young's modulus nor Poisson's 

    ratio have any role tp play.



    (2) Comparison with the equations of sections 3 and 5



    * The above solution has been arrived at by direct physical 

    arguments; but it must also be consistent with, and derivable from, 

    the above-derived general equations. This will now be checked.


    * Equation 3.1 enables the stress in the x-direction, X", to be 

    deduced as:


      X" = ex*2*G  +  L*e  -   H*T                   [3.1c]


    Substitution from the solutions for ex, D and T leads to:


      X" = 2*G*alpha*T  + L*3*alpha*T  - H*T


    whereafter introduction of the definitions of G, L and H leads to:


      X" = alpha*T*E * ( 1/(1+P) + 3*P/((1+P)*(1-2P)) - 1/(1-2P) )  

         = alpha*T*E * ( 1-2P + 3P - (1+P) ) / ((1+P)*(1-2P)) )  

         = alpha*T*E * ( 0 ) / ((1+P)*(1-2P)) )  

         = 0 


    as expected.  


    * Equation 5.13 can be checked as follows:


      (L+2*G)*u  +  (L+G)*(j - k )  -  H*T   +    X    =  0  [5.13a]

              ^            z   y         x                          


    As has already been explained in section 5.4.2, 


      (L+2*G)*u  =  - (L+2*G)*T*alpha

              ^               x



      (L+G)*(j - k ) =  4*(L+G)*T*alpha        

             z   y              x


    Since X. the externally-applied force equals zero, the equation 

    reduces to:


      [ - (L+2*G) + 4*(L+G) - H ] * T  = 0



    Omission of the factor T , and introducing the definitions of L, G 


    and H, enable theis equation to be written as:


          3*L + 2*G - H = 0 


          3*P/[(1+P)*(1-2P)] + 1/(1+P) - 1/(1-2P)  = 0


          3*P + 1-2*P - 1 - P = 0


          (3 - 2 - 1)*P + 1 - 1 = 0


    which is evidently true.                                 


    (b) Linear temperature profile in a laterally-constrained block



    To be supplied.



    7.5 Non-uniform bending stress



    ...... to be continued ................



    8. Solutions in terms of 'stress-functions'



    8.1 Introduction



    Although the differential equations derived above, with the dependent

    variables u, v, w, i, j, k and T, provide a complete description of 

    the general three-dimensional stress-strain-distribution problem, it 

    is sometimes easier to solve other differential equations which can 

    be derived from them.


    Such equations are those in which the same laws are re-formulated in 

    terms of derived equations of which the dependent variables are 

    scalars, known as 'stress-functions.


    In general, three such functions are required; and their 

    differential equations are of fourth order. However there are two 

    practically important cases in which a single second-order 

    differential equation suffices.


    These cases involve the two-dimensional stress-strain distributions 

    caused by:


    (a) the torsion of straight bars, and 

    (b) the bending of beams.


    Both are of practical importance in engineering; and their study 

    also throws valuable light on how boundary conditions are to be 

    formulated in general three-dimensional problems when the rotation 

    equations have to be solved.


    They are discussed in sections 8.1 and 8.2 respectively.


    In both cases the treatment follows closely that which is to be 

    found in T&G, chapters 10 and 11, which itself is modelled on the 

    1855 treatise of St Venant.


    St Venant employed the so-called 'semi-inverse method', which 

    involved assuming analytical expressions for the strain components 

    and proving that these conformed with the fundamental differential 

    equations if the forces were appropriately applied.



    8.2 The torsion of straight bars of non-uniform cross-section



    a. Expressions for the strain components



    St Venant's assumptions were, when z measures the distance along the 

    length of the bar and x and y are the distances measured within its 



    (1) The u and v displacements within the cross-section are given by:

        u = - theta * z * y                    (8.2 - 1),

        v =   theta * z * x                    (8.2 - 2),



    (2) the w displacement is given by:

        w =   theta * psi(x,y)                 (8.2 - 3),


    wherein theta is the rotation angle per unit length, and psi is a 

    function to be determined, which describes the 'warping of  

    cross-sectional surfaces.


    These assumptions imply, by reason of the definitions of section 1



        ex = ey = ez = exy = 0                 (8.2 - 4)


        ezx = theta * ( psi - y )              (8.2 - 5)



        eyz = theta * ( psi + x )              (8.2 - 6)



    Expressed in terms of the components of rotation, these assumptions

    imply, again by reason of the definitions of section 1:


        i = theta * ( x - psi )                (8.2 - 7),



        j = theta * ( psi + y )                (8.2 - 8),



        k = - theta * 2 * z                    (8.2 - 9).



    b. Corresponding expressions for the stress components



    T&G show further that the components of stress can be 

    correspondingly expressed as:


        sx = sy = sz = sxy = 0                 (8.2 - 10)


        szx = G * theta * ( psi - y )          (8.2 - 11)



        syz = G * theta * ( psi + x )          (8.2 - 12)



    At the lateral boundary of the bar, because there are no external 

    forces, the szx and syz components of the shear stress are linked by 

    the relation:


        szx * l  + syz * m  =  0               (8.2 - 13)


     where i and m are the direction cosines of the outward-directed 

     normal to the boundary surface.   




    c. (Alternative) equations to be solved



    T&G show that the equations of section 2 imply that psi obeys the 

    differential equation:


        psi    +    psi   =   0                (8.2 - 14)

         xx          yy


    However, they also show that the boundary conditions are simpler for

    a similar equation, derived from (8.2-14), with the dependent variable 

    phi. This is:     


        phi    +    phi   =   F                (8.2 - 15)

         xx          yy


    wherein the 'stress-function' phi is related to the shear stresses 



        szx  =   phi                           (8.2 - 16)



        syz  = - phi                           (8.2 - 17)




         F   = - 2 * G * theta                 (8.2 - 18)


    The boundary condition for phi at the bar boundary then reduces 

    simply to:


         phi = 0                               (8.2 - 19) .


    Finally it is shown that the twisting moments at the two ends of the 

    bar are given by:


        moment = 2 * integral phi * d area     (8.2 - 20)    


    It will be found useful, below, to express the rotation components 

    in terms of the phi. This can be done as follows:


    * addition of 8.2-6 to 8.2-7 leads to:


        i =  2 * theta * x  -  eyz        

          =  2 * theta * x  -  syz / G       


    * subtraction of 8.2-5 from 8.2-8 leads to:


        j =  2 * theta * y  +  ezx

          =  2 * theta * y  +  szx / G


    * substitution for syz and szx from 8.2-16 and 8.2-17 then leads to:


        i = 2 * theta * x  -  phi / G             (8.2 - 21)



        j = 2 * theta * y  +  phi / G             (8.2 - 22)




    d. Solution by way of PHOENICS



    Equation (8.2-15), with boundary condition (8.2-19), is of course 

    easily solved by means of PHOENICS, no matter what shape the solid 

    posesses; and indeed the more general equation in which the material 

    properties can vary with x and y, namely:


          ( G phi )  +  ( G phi )  =  F        (8.2 - 23) 

               x    x        y    y   


    is equally easy.


    [ Note that such a case appears on p.160 of O.Zienkiewicz's 1967 

    book entitled The Finite-Element Method; it will be instructive to 

    compare the solution given there with that computed by way of



    All that is necessary is to:

    * declare a scalar variable called phi;

    * de-activate its convection and time-dependent terms by way of the 

      TERMS command;

    * impose the phi=0 condition on its boundary;

    * supply a uniform source per unit volume to cells within is 

      boundary; and

    * activate the linear-equation solver.    



    e. The bar of elliptic cross-section



    Analytical solutions are provided for several bar cross-sections 

    by T&G; these can also serve as tests of the accuracy of solutions 

    computed by PHOENICS.


    That for the bar of elliptical cross-section, of which the semi-axis 

    lengths are a and b, is expressible as:


               a**2 * b**2 * F   ( x**2   y**2     )

         phi = --------------- * ( ---- + ---- - 1 )    (8.2-24)

               2*(a**2 + b**2)   ( a**2   b**2     )



         F   = - 2 * moment * (a**2 + b**2) / ( pi * a**3 * b**3)  



         szx = - 2 * moment * y / ( pi * a * b**3 )     (8.2-26)


         syz =   2 * moment * x / ( pi * a**3 * b )     (8.2-27)


    Differentiation of 8.2-24 with respect first to x and then to y 

    leads to:


         phi = b**2 * F * x / (a**2 + b**2)    



         phi = a**2 * F * y / (a**2 + b**2)



    Substitution of these expressions into 8.2-21 and 8.2-22, with F 

    substituted from  8.2-18, leads to:


          i = 2 * theta * x * a**2 / ( a**2 + b**2 )    (8.2-28) 


          j = 2 * theta * y * b**2 / ( a**2 + b**2 )


    The third rotation component k, meanwhile, continues to obey the

    already-presented equation, valid for all cross-sections:                                                    


        k = - 2 * theta * z                    (8.2 - 9).


    It thus appears that the distributions of i, j and k are all 

    linear and have magnitudes which are proportional to the distance 

    from the origin of coordinates.


    e. Discussion



    It is interesting to consider how the stresses and strains could

    have been calculated WITHOUT recourse to the stress-function 

    technique, i.e. directly from the equations for the rotation 

    components and the displacements embodied in PHOENICS.


    (1) Differential equation 4.1 for k is easily solved; for the source 

    K is zero within the bar; the boundary conditions are k=0 at z=0 

    and k=-2*theta*z_end at the end of the bar; and there are no forces 

    at the lateral boundary which can act as sources of k. 


    [ Here it must be remembered the k is defined as TWICE the true 

    rotation, and that theta represents the TWIST PER UNIT LENGTH. ] 


    Therefore numerical solution of the k equation will lead to the 

    linear distribution of k expressed by 8.2-9.


    It is noteworthy that the gradient of k is finite only in the z 

    direction, and that this gradient does not appear as a source in any 

    of the differential equations for the displacements 5.13, 5.14 or 



    (2) Let us assume (and check later) that gradients of i and j are 

    finite only in the  x and y directions respectively and so also make 

    no appearance in the equations for the displacements.


    (3) Then the equations for u and v, namely  5.13 and 5.14, are also 

    easily solved numerically; for their values a z=0 and z=z_end are 

    known, from the imposed twist; and their are no forces in their 

    respective directions imposed anywhere else.


    Their numerical solutions will undoubtedly agree with equations 

    8.2-1 and 8.2-2 above.


    (4) There also appear to be no forces in the z-direction which could 

    affect the w equation. However, the boundary shear-stress conditions 

    represented by 8.2-11, 8.2.12 and 8.2-13 are relevant. 


    These may be usefully re-written in terms of ezx, eyz and gradients 

    of w as:

        ezx = w - theta*y ; eyz= w - theta*x ; ezx*l + eyl*m =0

              x                  y


    It thus follows that there exist easy-to-implement gradient boundary 

    conditions for the w equation which allows numerical solution of the 

    equation to proceed.


    (5) With u, v and w now computed, the rotations i and j can now be       

    evaluated; and it will undoubtedly be found that these quantities 

    have gradients in the x- and y-directions respectively, so 

    justifying the assumption (2).


    (6) It may be concluded that solution of the PHOENICS solid-stress 

    equation system can indeed be conducted straightforwardly for 

    problems in volving bars of arbitrary cross-section. 




    8.3 The bending of beams of non-uniform cross-section