******************************************************************** * * * The equations of displacements, stresses and strains in solids * * * ******************************************************************** by: D B Spalding Contents -------- 1. Nomenclature 2. Strain related to stress 3. Normal stress related to strain 4. Balance of forces 5. Differential equations for the displacements 6. The differential equations for the displacements; 2D 7. The differential equations for the displacements; 1D 8. Cylindrical geometry 9. The representation of the differential equations in PHOENICS 10. Concluding remarks -------------------------------------------------------------------- 1. Nomenclature --------------- * Cartesian coordinate directions: x, y, z * Radius, axial distance and angle in cylindrical-polar cordinates : r, a, z * Strains (extensions) in x, y, and z directions: ex, ey, ez * Linear thermal expansion: et * Normal stresses in x, y, z directions, divided by Young's modulus, presumed uniform: nx, ny, nz * Shear stresses in xy, xz, yz planes, divided by Young's modulus: sxy, sxz, syz * Displacements in x, y and z directions: u, v, w * Poisson's ratio: R * Various constant functions of Poisson's ratio: R1, R2, R3, R4 ... namely: R1 = (1 - R) / [(1 + R) (1 - 2 R)] R2 = R / [(1 + R) (1 - 2 R)] R3 = 1 / (1 - 2 R) R4 = 1 / [ 2 (1 + R)] * derived functions: R5 = R1 - R2 R2 / R1 = 1 / (1 + R) R6 = R2 - R2 R2 / R1 = R / [(1 + R) (1 - R)] R7 = R3 - R3 R2 / R1 = 1 / (1 - R) R8 = R1 - 2 R2 R2 / ( R1 + R2) ] R9 = R3 [ 1 - 2 R2 R2 / ( R1 + R2) ] R10 = 1 / [ 2 ( 1 + R ) ( 1 - 2 R) ] R11 = 1 / R4 = 2 (1 + R) R12 = 2 ( 1 + R ) / ( 1 - R ) R13 = 2 ( 1 + R ) (1 - R ) R14 = 2 ( 1 - R ) Note that, if R = 0, then R1 = R3 = R5 = R8 = R11 = 1.0, R2 = R6 = 0.0 , R4 = R10 = 0.5 and R12= 2.0 * First-order differentiation operators with respect to x, y, z, a: ddx., ddy., etc * Second-order differentiation operators with respect to x & x y & y z & z x & y x & z and y & z: d2dxx. d2dyy. d2dzz. d2dxy. d2dxz. d2dyz. * The divergence-of-gradient operator: div_grad.= d2dxx. + d2dyy. + d2dzz * The multiplication operator is implied by any space in an equation which is not adjacent to one of the operators: + - / = or a differential operator. * External forces per unit volume, including gravity, aceleration, etc, in the x, y and z directions: fx, fy , fz 2. Strain related to stress --------------------------- The following are expressions of Hooke's and Poisson's experimentally-based laws: 2.1 Normal strains related to normal stresses --------------------------------------------- ex = nx - R ( ny + nz ) + et (2.1-1) ey = ny - R ( nx + nz ) + et (2.1-2) ez = nz - R ( nx + ny ) + et (2.1-3) [ These correspond to the top part of matrix D on page 58 of Huebner & Thornton (H&T form now on), except that et (thermal strain) is not mentioned by them. In the absence of any stresses, they correspond to ex=ey=ez=et, which seems reasonable. ] 2.2 Shear stresses related to shear strains ------------------------------------------- sxy = R4 ( ddy.u + ddx.v ) = syx (2.2-1) sxz = R4 ( ddz.u + ddx.w ) = sxz (2.2 2) szy = R4 ( ddy.w + ddz.v ) = szy (2.2-3) [ The origin of these equations is not clear. They appear to correspond to the bottom part of the C matrix of H&T. But the coefficient s there are (1-2R)/2, whereas R4 = 1/{2(1+R)} ] 3. Normal stress related to strain normal strain ------------------------------------------------ 3.1 General 3D equations ------------------------ The following are derived from the above by algebraic re-arrangement: nx = R1 ex + R2 ( ey + ez ) - R3 et (3.1-1) ny = R1 ey + R2 ( ex + ez ) - R3 et (3.1-2) nz = R1 ez + R2 ( ex + ey ) - R3 et (3.1-3) The details of the derivation, for nx only, are as follows: Addition of R times equations 2.1-2 and 2.1-3 to (1 - R) times equation 2.1-1 leads to: (1 - R) ex + R ey + R ez = (1 - R) nx - (1 - R) R (ny + nz) + (1 - R) et + R ny - R R (nx + nz) + R et + R nz - R R (nx + ny) + R et (3.1-4) in which the multipliers of both ny and nz sum to zero. Hence, (1 - R) ex + R ey + R ez = (1 - R - 2 R R) nx + (1 + R) et (3.1-5) ie (1 - R) ex + R ey + R ez = (1 + R) (1 - 2 R) nx + (1 + R) et (3.1-6) ie nx = [ (1 - R)/ {(1 + R) (1 - 2 R)} ] ex + [ R / {(1 + R) (1 - 2 R)}] ( ey + ez) - [ 1 / (1 - 2 R)} ] et (3.1-7) ie nx = R1 ex + R2 ( ey + ez ) - R3 et (3.1-1) [ Note that the shear-stress equations (1.2-1.2.3) have not been used here.] 3.2 2D equations ---------------- If ez (for example) is zero everywhere, which implies that the process is two-dimensional. movement in the third dimension being prevented, then the equations for nx and ny become: nx = R1 ex + R2 ey - R3 et (3.2-1) ny = R1 ey + R2 ex - R3 et (3.2-2) nz = R2 ( ex + ey ) - R3 et (3.2-3) while the equation for nz remains as above. If however there is no constraint in the z direction, ez may be non-zero, while nz obeys: nz = 0 . (3.2-3) Then ez may be eliminated from the equations of section 3.1, with the result: nx = R5 ex + R6 ey - R7 et (3.2-4) ny = R5 ey + R6 ex - R7 et (3.2-5) The derivation, for nx only, is as follows: From 3.1-3, with nz = 0, leads to: ez = - (R2/R1) ( ex + ey ) + (R3/R1) et (3.2-6) Substitution of ez, from 3.2-6 into 3.1-1, yields: nx = R1 ex + R2 ( ey - (R2/R1) (ex + ey)) + (R2 R3 /R1) et - R3 et ie nx = (R1 - R2 R2/R1) ex + R2 (1 - R2/R1) ey - R3 (1 - R2/R1) et (3.2-7) ie, from the definitions of R5, R6 and R7: nx = R5 ex + R6 ey - R7 et (3.2-4) Expressed in terms of R, this is: nx = [ 1 / (1 + R) ] ex + { R / [(1 + R) (1 - R)] } ey - [ 1 / (1 - R) ] et (3.2-4) 3.3 1D equations ---------------- If ey and ez (for example) are zero everywhere, the equations for nx, ny and nz become: nx = R1 ex - R3 et (3.3-1) ny = R2 ex - R3 et (3.3-2) nz = R2 ex - R3 et (3.3-3) If however there are no constraints in the y and z direction, ey and ez may be non-zero, while ny and nz obey: ny = 0 and nz = 0 ; then ey and ez may be eliminated from the equations of section 3.2, with the result: nx = R8 ex - R9 et (3.3-4) The derivation is as follows. Equations 3.1-1. 3.1-2 and 3.1-3 become: nx = R1 ex + R2 ( ey + ez ) - R3 et (3.3-5) 0 = R1 ey + R2 ( ex + ez ) - R3 et (3.3-6) 0 = R1 ez + R2 ( ex + ey ) - R3 et (3.3-7) Since there is nothing to distinguish ey from ez, they may be taken as equal, with the results: nx = R1 ex + 2 R2 ey - R3 et (3.3-5) 0 = ( R1 + R2 ) ey + R2 ex - R3 et (3.3-6) Substitution for ey from 3.3-6 in 3.3-5 yields: nx = R1 ex + [ 2 R2 / (R1 + R2) ] [ - R2 ex + R3 et] - R3 et (3.3-7) ie nx = [ R1 - 2 R2 R2 / ( R1 + R2) ] ex - R3 [ 1 - 2 R2 R2 / ( R1 + R2) ] et (3.3-8) In terms of R, this is: 4. Balance of forces -------------------- The balances of the forces acting on unit volume in the three coordinate directions lead to: ddx. nx + ddy. sxy + ddz.sxz + fx = 0 (4.-1) ddy. ny + ddx. sxy + ddz.syz + fy = 0 (4.-2) ddz. nz + ddx. sxz + ddy.syz + fz = 0 (4.-3) 5. The differential equations for the displacements; 3D ------------------------------------------------------- 5.1 The x-direction displacement u ---------------------------------- Differentiation of the equation 4.-1 leads, after substitution from the equations of section 3.1, to: ddx. [ R1 ex + R2 ( ey + ez ) - R3 et ] + ddy. [ R4 ( ddy.u + ddx.v) ] + ddz. [ R4 ( ddz.u + ddx.w) ] + fx = 0.0 Since, by definition: ex = ddx. u ey = ddy. v ez = ddz. w the strains ex, ey, ez can be eliminated from the above differential equation, with the result: ddx. [ R1 ddx. u + R2 ( ddy. v + ddz. w ) - R3 et ] + ddy. [ R4 ( ddy. u + ddx. v) ] + ddz. [ R4 ( ddz. u + ddx. w) ] + fx = 0.0 If Poisson's ratio R is taken as constant, implying that all Rs are also constant, this equation can be written in terms of second-order operators, as follows: ddx. [ R1 ddx. u + R2 ( ddy. v + ddz. w ) - R3 et ] + R4 d2dyy. u + R4 ddx. ddy. v + R4 d2dzz. u + R4 ddx. ddz. w + fx = 0.0 with further rearrangement to: ddx. [ ( R1 - R4 ) ddx. u + ( R2 + R4 ) ( ddy. v + ddz. w ) ] - ddx. [ R3 et ] + R4 ( d2dxx. u + d2dyy. u + d2dzz. u ) + fx = 0.0 But, as a consequence of the above definitions, R1 - R4 = R2 + R4 = 1 / [2 (1 + R) (1 - 2 R)] = R10 Therefore, the equation can be reduced to: R10 ddx. dil - R3 ddx et + R4 div_grad. u + fx = 0 where dil = ddx. u + ddy. v + ddz. w , the dilatation. Re-arrangement so as to give prominence to the div_grad term, leads to: div_grad. u + ddx. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ] + fx / R4 = 0.0 (5.1-1) [ This should be compared with equation C.15 of H&T. In that equation, the multiplier of dil is 1/(1-2R). R10/R4 = 2(1+R)/{2(1+R)(1-2R)} i.e. 1/(1-2R) . So they equations agree. The multiplier of fx in H&T is 2(1+R)/E, which is again in agreement. However, H&T are silent about et .] 5.2 The y- and z-direction displacements v and w ------------------------------------------------ Replacing x successively by y and z, and u by v and w, leads to the corresponding equations: div_grad. v + ddy. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ] + fy / R4 = 0.0 div_grad. w + ddz. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ] + fz / R4 = 0.0 Note that, as a consequence of the definitions, R10 / R4 = 1 / ( 1.0 - 2 R ) = R3 and R3 / R4 = 2 ( 1 + R ) / ( 1 - 2 R ) The three equations can thus be written as: -------------------------------------------------------------------- |u| |x| |x| div_grad. |v| + R3 dd|y|. ( dil - R11 et ) + R11 f|y| = 0.0 |w| |z| |z| (5.2-1) -------------------------------------------------------------------- where it is to be noted that: R3 = 1 /( 1 - 2 R) R11 = 1 / R4 = 2 ( 1 + R ) 6. The differential equations for the displacements; 2D ------------------------------------------------------- 6.1 When displacements in the third direction are zero ------------------------------------------------------ Let the zero-displacement direction be z. Then the equation of section 5.2 still holds, for directions x and y, with the proviso that ddz. terms are absent from both the div_grad and dil definitions. 6.2 When stressess in the third direction are zero -------------------------------------------------- When nz, say, is zero, the analysis must start from nx = R5 ex + R6 ey - R7 et rather than from nx = R1 ex + R2 ey - R3 et . Thereafter, apart from the substitutions of Rn+4 for Rn, and from the absence of z-direction derivatives, the analysis proceeds as in section 6. The resulting equations are: -------------------------------------------------------------------- |u| |x| |x| | div_grad. |v| + R3 dd|y|. ( dil - R11 et ) + R11 f|y| = 0.0 (6.2-1) -------------------------------------------------------------------- whereby it should be noted that R11 remains in place, but R13 has replaced R12. 7. The differential equations for the displacements; 1D ------------------------------------------------------- 7.1 When displacements in the second and third directions are zero ------------------------------------------------------------------ Let the zero-displacement directions be y and z. Then the equation of section 5.2 still holds, for directions x and y, with the proviso that ddz. terms are absent from both the div_grad and dil definitions. 7.2 When stressess in the second and third directions are zero -------------------------------------------------------------- When ny and nz, say, are zero, the analysis must start from: nx = R8 ex + R9 et . Thereafter, the analysis proceeds as in section 6. The resulting equation is: ------------------------------------------------------------------ | | div_grad. u + R11 ddx. ( dil + R14 et ) + fx / R4 = 0.0 | (7.2-1) ------------------------------------------------------------------ 8. Cylindrical geometry ----------------------- 8.1 Nomenclature ---------------- When the grid is cylindrical polar, and x, y and z are replaced by a (ie angle), r (ie radial distance) and z (now axial distance), * the strains become: ea, er, ez * the normal stresses become: na, nr, nz * the displacements become: u, v, w 8.2 Strain related to stress ---------------------------- The following relations replace those of section 2 above: ea = na - R ( nr + nz) + et er = nr - R ( na + nz) + et ez = nz - R ( na + nr) + et sar = R4 ( ddr.u - u / r + dda.v / r) = sra saz = R4 ( ddz.u + dda.w / r) = sza szr = R4 ( ddr.w + ddz.v ) = srz 8.3 Normal stress related to strain, general 3D ----------------------------------------------- Because the e-to-n relations are the same as before (ie in section 2), so are the n-to-e relations (ie as in section 3). 8.4 Balance of forces --------------------- In the angular direction, the balance of forces becomes: dda. na / r + ddr. sar + ddz. saz + fa + ( sra + sar ) / r = 0.0 . In the radial direction, the balance of forces becomes: ddr. ( nr r ) / r + dda. sar / r + ddz. szr - na / r + fr = 0.0 , which may also be expressed as: ddr. nr + dda. sar / r + ddr. szr + ( nr - na ) / r + fr = 0.0 . In the axial direction, the balance of forces becomes: ddz. nz + ddr. ( srz r) / r + dda. saz / r + fz = 0.0 , which may also be written as: ddz. nz + ddr. srz + dda. saz / r + srz / r + fz = 0.0 . Thus, as compared with the corresponding equations appearing section 3, the following extra terms appear: - in the equation for the angular (formerly x) direction: none; - in the equation for the radial (formerly y) direction: + ( nr - na ) / r ; - in the equation for the axial (formerly, and still, z) direction: + srz / r . 8.5 The differential equations for the displacements ---------------------------------------------------- Since the differences in the differential equations from their cartesian-coordinate counterparts all derive from the above- mentioned extra terms, it can be deduced that the differential equations in cylindrical-polar coordinates can be expressed as: div_grad. u + R11 dda. ( dil + R12 et ) / r + fa / R4 = 0.0 div_grad. v + R11 ddr. ( dil + R12 et ) + [fr + ( nr - na ) / r ] / R4 = 0.0 div_grad. w + R11 ddz. ( dil + R12 et ) + (fz + srz / r ) / R4 = 0.0 9. The representation of the differential equations in PHOENICS -------------------------------------------------------------- 9.1 Comparison of the equations for displacements and velocities ---------------------------------------------------------------- The momentum equations which are solved by PHOENICS, when the viscosity is set to unity and the convection terms are neglected, can be represented as: ---------------------------------------------------------- |u| |x| |x| div_grad. |v| - dd|y|. ( p ) + mom_so|y| = 0.0 |w| |z| |z| (9.1-1) ---------------------------------------------------------- where p stands for pressure, and mom_so stands for momentum source per unit volume. Comparison of these equations for those for the displacements shows that they can be made identical if the pressure is related to the dilatation and to the thermal expansion by: - p = R3 ( dil - R11 et ) ie - p = [1 / (1 - 2 R)] [ dil - 2 ( 1 + R ) et ] (9.1-2) and if: |x| = |x| mom_so|y| = R11 f|y| |z| = |z| ie: |x| = |x| mom_so|y| = 2 ( 1 + R ) f|y| (9.1-3) |z| = |z| 9.2 The dilatation in PHOENICS ------------------------------ In PHOENICS, the dilatation, ie ddx. u + ddy. v + ddz. w , appears as a volumetric source in the mass-conservation equation. Therefore the required relationship can be contrived by providing, for the part of the domain which is occupied by solids, a volumetric source term calculated as: vol_sor = dil = - (1 / R3) p + R11 et ie: vol_sor = dil = - (1 - 2 R) p + 2 (1 + R ) et (9.2-1) 9.3 Normal-stress boundary conditions ------------------------------------- The expression for nx in equation 3.1-1 can be re-written in terms of the dilatation as follows: nx = ( R1 - R2 ) ex + R2 ( ex + ey + ez ) - R3 et ie nx = ( R1 - R2 ) ex + R2 dil - R3 et (9.3-1) Replacement of dil by p, from equation 9.1-2, then leads to: nx = ( R1 - R2 ) ex + R2 [- (1 - 2 R) p + 2 (1 + R ) et ] - R3 et (9.3-2) This can be rearranged with ex on the left-hand side and replaced by its equivalent, ddx.u, as: ddx. u = [ 1 / (R1 - R2) ] nx + [ R2 (1 - 2 R) / (R1 - R2 )] p + [ R3 - 2 R2 (1 + R) / ( R1 - R2 ) ] et (9.3-3) which becomes, after simplification: ddx. u = ( 1 + R) nx + R p + ( 1 + R ) et (9.3-4) Since ddx.u ie the gradient of the displacement u, this equation can be used as the normal-stress boundary condition of the differential equation for u. It can be represented by way of a source term in the computational cell within the solid and adjacent to the fluid, as illustrated below. . | . | . /| . | . | . /| . | . | . /| -----*------>------*------>-----*---->=====> normal stress . | . | . /| . | solid | solid /| fluid . | . | . /| Corresponding equations can be derived for the v and w boundary conditions, namely: ddy. v = ( 1 + R) ny + R p + ( 1 + R ) et (9.3-5) ddz. w = ( 1 + R) nz + R p + ( 1 + R ) et (9.3-6) A consequence is that, when a direct stress S is to be exerted at a boundary, the RATCH and COVAL statements in the RHOENICS Q1 file take the form: PATCH(name,north/east/high,,,,,,) COVAL(name,v1/u1/w1,fixflu,S*(1+R)) A further consequence is that, within EARTH, sources per unit area are applied to cells adjacent to solid-fluid boundaries which are equal to: R p + ( 1 + R ) et . 10. Concluding remarks ---------------------- Further topics to be dealt with in the continuation of the current document include:- * allowance for non-uniform material properties; * the relationship of the dilatation to the PHOENICS pressure variable; * the solution algorithm; * a series of exact solutions of the equations, which can be used as tests; * a series of practically-interesting problems, which can be used so as to interest engineers in the use of the technique. It should also be remarked that all the above derivations require to be checked carefully and independently.