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