******************************************************************** * * * The equations for displacements, stresses and strains in solids * * * ******************************************************************** by: D B Spalding (latest update 04.04.06) Contents -------- 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 x likwise du/dy is written as u y likewise du/dz is written as u z * 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] x sy = L*e + 2*G*v - H*T [2.5] y sz = L*e + 2*G*w - H*T [2.6] z 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] s m = y [3.1-5] s n = z [3.1-6] s 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] x 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] y 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] z 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 differential. (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] ^ and 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 elasticity. 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, thus: k = k1 * x/x1 where the length of the beam is x1 and the rotation at x=x1 equals k1. 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 components: 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 i.e. 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 weaker. 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 and: (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 ! ! ! ! ! ! --> ! ! ! ! ^ ! !/////////! ! ! ! ! ! !/////////! ! ! !--dyn----!----k----^----kn---^----k----!---^-----!------------------ ! ! ! NW/////////NE ! ! ! ! ! ! W !/// P ///! E ! ! ! ! v ! --> ! --> ! --> ! dyp ! ! ^ ! !\\\\\\\\\! ! ! ! ! ! ! !\\\\\\\\\! ! ! ! !--dys----!----k----^----ks---^----k----!---v-----!------------------ ! ! ! SW\\\\\\\\\SE ! ! ! ! ! !\\\ S \\\! ! ! ! v ! ! --> ! ! ! ! ! !<--dxp-->! ! ! ! ! <--dxw--><-- dxe--> ! ! Consider the terms of equation 5.2-1 which concern displacements, namely: (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 xy-plane. 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 thus: sxy = G*exy = G*( u + v ) [2.7] y x kn = 2*dy/dy - sxy whence: 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 wherein: * 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 ry u stands for [d(u)/dx]/r rx ru stands for [d(ru)/dy]/r ry v stands for [d(v)/dx]/r rx 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, namely: * the dimensions which are NOT the 'one', or one of the 'two', need NOT BE CONSIDERED AT ALL. 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 displacements. 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 movement. 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) i.e. ex = ( - L*e + H*T) / ( 2*G + bx ) [5.16] ey = ( -by*ey - L*e + H*T ) / (2*G) i.e. ez = ( - L*e + H*T) / ( 2*G + by ) [5.17] ez = ( -bz*ez - L*e + H*T ) / (2*G) i.e. 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 . x 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] x 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] xx (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] x 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] zz (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: z 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 xz j = 4*a/l [5.27] z A similar argument shows that k has the same value; and the upshot y 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 . x Similar conclusions can be drawn about the equations for the y and x directions. 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, or * 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 ) wherefore 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 definition. ** 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] whence 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 expansion. 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, is: 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 is: 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 substitutions: 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 z-directions; (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 ) 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] x 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 and (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 x Omission of the factor T , and introducing the definitions of L, G x and H, enable theis equation to be written as: 3*L + 2*G - H = 0 i.e. 3*P/[(1+P)*(1-2P)] + 1/(1+P) - 1/(1-2P) = 0 i.e. 3*P + 1-2*P - 1 - P = 0 i.e. (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 cross-section: (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), and (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 above: ex = ey = ez = exy = 0 (8.2 - 4) and ezx = theta * ( psi - y ) (8.2 - 5) x eyz = theta * ( psi + x ) (8.2 - 6) y 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), y j = theta * ( psi + y ) (8.2 - 8), x 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) and szx = G * theta * ( psi - y ) (8.2 - 11) x syz = G * theta * ( psi + x ) (8.2 - 12) y 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 by: szx = phi (8.2 - 16) y syz = - phi (8.2 - 17) x and 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) x j = 2 * theta * y + phi / G (8.2 - 22) y 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 PHOENICS ]. 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) (8.2-25) 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) x and phi = a**2 * F * y / (a**2 + b**2) y 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) and j = 2 * theta * y * b**2 / ( a**2 + b**2 ) (8.2-29) 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 5.15. (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 -----------------------------------------------------