Closed-Form Solutions

Closed-Form Solutions February 21, 2001 Introduction During most of its historical development, the science of Mechanics of Materials relied principally on closed-form (not computational) mathematical theorists. Much of their work represents mathematical intuition and skill of a very high order, challenging even for advanced researchers of today. This theory is taught primarily in graduate subjects, but is outlined here both to provide some background that will be useful in the Module on Fracture and as a preliminary introduction to these more advanced subjects. Governing equations We have earlier shown (see Module 9) how the spatial gradients of the six Cauchy stresses are related by three equilibrium equations that can be written in pseudovector form as LT σ = 0 (1) These are augmented by six constitutive equations which can be written for linear elastic materials as (see Module 11) σ = D (2) and six kinematic or strain-displacement equations (Module 8)  = Lu (3) These fifteen equations must be satisfied by the fifteen independent functions (three displacements u, six strains , and six stresses σ). These functions must also satisfy boundary conditions on displacement u = û on Γu (4) where Γu is the portion of the boundary on which the displacements u = û are prescribed. The remainder of the boundary must then have prescribed tractions T = T̂ , on which the stresses must satisfy Cauchy’s relation: σn̂ = T̂ on 1 ΓT (5) In the familiar cantilevered beam shown in Fig. 1, the region of the beam at the wall constitutes Γu , having specified (zero) displacement and slope. All other points on the beam boundary make up ΓT , with a load of P at the loading point A and a specified load of zero elsewhere. Figure 1: Cantilevered beam. With structures such as the beam that have simple geometries, solutions can be obtained by the direct method we have used in earlier modules: an expression for the displacements is written, from which the strains and stresses can be obtained, and the stresses then balanced against the externally applied loads. (Problem 2 provides another example of this process.) In situations not having this geometrical simplicity, the analyst must carry out a mathematical solution, seeking functions of stress, strain and displacement that satisfy both the governing equations and the boundary conditions. Currently, practical problems are likely to be solved by computational approximation, but it is almost always preferable to obtain a closed-form solution if at all possible. The mathematical result will show the functional importance of the various parameters, such as loading conditions or material properties, in a way a numerical solution cannot, and is therefore more useful in guiding design decisions. For this reason, the designer should always begin an analysis of load-bearing structures by searching for closed-form solutions of the given, or similar, problem. Several compendia of such solutions are available, the book by Roark 1 being a useful example. However, there is always a danger in performing this sort of “handbook engineering” blindly, and this section is intended partly to illustrate the mathematical concepts that underlie many of these published solutions. It is probably true that most of the problems that can be solved mathematically have already been completed; these are the classical problems of applied mechanics, and they often require a rather high level of mathematical sophistication. The classic text by Timoshenko and Goodier2 is an excellent source for further reading in this area. The Airy stress function Expanding the kinematic or strain-displacement equations (Eqn. 3) in two dimensions gives the familiar forms: x = y = γxy = 1 2 ∂u ∂x ∂v ∂y ∂v ∂u + ∂x ∂y W.C. Young, Roark’s Formulas for Stress and Strain, McGraw-Hill, New York, 1989. S. Timoshenko and J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1951. 2 (6) Since three strains (x , y , γxy ) are written in terms of only two displacements (u, v), they cannot be specified arbitrarily; a relation must exist between the three strains. If x is differentiated twice by dx, y twice by dy, and γxy by dx and then dy we have directly ∂ 2 x ∂ 2 y ∂ 2 γxy + = ∂y 2 ∂x2 ∂x ∂y (7) In order for the displacements to be so differentiable, they must be continuous functions, which means physically that the body must deform in a compatible manner, i.e. without developing cracks or overlaps. For this reason Eqn. 7 is called the compatibility equation for strains, since the continuity of the body is guaranteed if the strains satisfy it. The compatibility equation can be written in terms of the stresses rather than the strains by recalling the constitutive equations for elastic plane stress: 1 (σx − νσy ) E 1 y = (σy − νσx ) E 1 2(1 + ν) = τxy = τxy G E x = γxy (8) Substituting these in Eqn. 7 gives ∂2 ∂2 ∂ 2 τxy (σ − νσ ) + (σ − νσ ) = 2(1 + ν) x y y x ∂y 2 ∂x2 ∂x ∂y (9) Stresses satisfying this relation guarantee compatibility of strain. The stresses must also satisfy the equilibrium equations, which in two dimensions can be written ∂σx ∂τxy + =0 ∂x ∂y ∂τxy ∂σy + =0 ∂x ∂y (10) As a means of simplifying the search for functions whose derivatives obey these rules, G.B. Airy (1801–1892) defined a stress function φ from which the stresses could be obtained by differentiation: ∂2φ ∂y 2 ∂2φ = ∂x2 ∂2φ = − ∂x ∂y σx = σy τxy (11) Direct substitution will show that stresses obtained from this procedure will automatically satisfy the equilibrium equations. This maneuver is essentially limited to two-dimensional problems, but with that proviso it provides a great simplification in searching for valid functions for the stresses. 3 Now substituting these into Eqn. 9, we have ∂4φ ∂4φ ∂4φ + 2 + ≡ ∇2 (∇2 φ) ≡ ∇4 φ = 0 ∂y 4 ∂x4 ∂x2 ∂y 2 (12) Any function φ(x, y) that satisfies this relation will satisfy the governing relations for equilibrium, geometric compatibility, and linear elasticity. Of course, many functions could be written that satisfy the compatibility equation; for instance setting φ = 0 would always work. But to make the solution correct for a particular stress analysis, the boundary conditions on stress and displacement must be satisfied as well. This is usually a much more difficult undertaking, and no general solution that works for all cases exists. It can be shown, however, that a solution satisfying both the compatibility equation and the boundary conditions is unique; i.e. that it is the only correct solution. Stresses around a circular hole Figure 2: Circular hole in a uniaxially stressed plate To illustrate the use of the Airy function approach, we will outline the important work of Kirsch3 , who obtained a solution for the influence on the stresses of a hole placed in the material. This is vitally important in analyzing such problems as rivet holes used in joining, and the effect of a manufacturing void in initiating failure. Consider a thin sheet as illustrated in Fig. 2, infinite in lateral dimensions but containing a circular hole of radius a, and subjected to a uniaxial stress σ. Using circular r, θ coordinates centered on the hole, the compatibility equation for φ is  4 ∇ φ= ∂2 1 ∂ 1 ∂2 + + 2 2 2 ∂r r ∂r r ∂θ  ∂ 2 φ 1 ∂φ 1 ∂2φ + + 2 2 2 ∂r r ∂r r ∂θ  =0 (13) In these circular coordinates, the stresses are obtained from φ as 1 ∂φ 1 ∂2φ + 2 2 r ∂r r ∂θ ∂2φ = ∂r 2   ∂ 1 ∂φ = − ∂r r ∂θ σr = σθ τrθ 3 G. Kirsch, VDI, vol. 42, 1898; described in Timoshenko & Goodier, op. cit.. 4 (14) We now seek a function φ(r, θ) that satisfies Eqn. 13 and also the boundary conditions of the problem. On the periphery of the hole the radial and shearing stresses must vanish, since no external tractions exist there: σr = τrθ = 0, r=a (15) Far from the hole, the stresses must become the far-field value σ; the Mohr procedure gives the radial and tangential stress components in circular coordinates as σr = σθ = τrθ = σ 2 (1 + cos 2θ) σ 2 (1 − cos 2θ) σ 2 sin 2θ ⎫ ⎪ ⎬ ⎪ ⎭ r→∞ (16) Since the normal stresses vary circumferentially as cos 2θ (removing temporarily the σ/2 factor) and the shear stresses vary as sin 2θ, an acceptable stress function could be of the form φ = f (r) cos 2θ (17) When this is substituted into Eqn. 13, an ordinary differential equation in f (r) is obtained:  d2 1 d 4 + − 2 2 dr r dr r  d2 f 1 df 4f + − 2 2 r dr r dr  =0 This has the general solution 1 +D (18) r2 The stress function obtained from Eqns. 17 and 18 is now used to write expressions for the stresses according to Eqn. 14, and the constants determined using the boundary conditions in Eqns. 15 and 16; this gives f (r) = Ar 2 + Br 4 + C σ a4 σ a2 σ A = − , B = 0, C = − , D= 4 4 2 Substituting these values into the expressions for stress and replacing the σ/2 that was temporarily removed, the final expressions for the stresses are σr = σ 2 σθ = σ 2 τrθ   σ = − 2 a2 1− 2 r a2 1+ 2 r    σ + 2 σ − 2  3a4 4a2 1+ 4 − 2 r r  3a4 1+ 4 r 3a4 2a2 1− 4 + 2 r r  cos 2θ  cos 2θ (19)  sin 2θ As seen in the plot of Fig. 3, the stress reaches a maximum value of (σθ )max = 3σ at the periphery of the hole (r = a), at a diametral position transverse to the loading direction (θ = π/2). The stress concentration factor, or SCF, for this problem is therefore 3. The x-direction stress falls to zero at the position θ = π/2, r = a, as it must to satisfy the stress-free boundary condition at the periphery of the hole. 5 Figure 3: Stresses near circular hole. (a) Contours of σy (far-field stress applied in y-direction). (b) Variation of σy and σx along θ = π/2 line. Note that in the case of a circular hole the SCF does not depend on the size of the hole: any hole, no matter how small, increases the local stresses near the hole by a factor of three. This is a very serious consideration in the design of structures that must be drilled and riveted in assembly. This is the case in construction of most jetliner fuselages, the skin of which must withstand substantial stresses as the differential cabin pressure is cycled by approximately 10 psig during each flight. The high-stress region near the rivet holes has a dangerous propensity to incubate fatigue cracks, and several catastrophic aircraft failures have been traced to exactly this cause. Note also that the stress concentration effect is confined to the region quite close to the hole, with the stresses falling to their far-field values within three or so hole diameters. This is a manifestation of St. Venant’s principle4 , which is a common-sense statement that the influence of a perturbation in the stress field is largely confined to the region of the disturbance. This principle is extremely useful in engineering approximations, but of course the stress concentration near the disturbance itself must be kept in mind. When at the beginning of this section we took the size of the plate to be “infinite in lateral extent,” we really meant that the stress conditions at the plate edges were far enough away from the hole that they did not influence the stress state near the hole. With the Kirsch solution now in hand, we can be more realistic about this: the plate must be three or so times larger than the hole, or the Kirsch solution will be unreliable. Complex functions In many problems of practical interest, it is convenient to use stress functions as complex functions of two variables. We will see that these have the ability to satisfy the governing equations automatically, leaving only adjustments needed to match the boundary conditions. For this reason, complex-variable methods play an important role in theoretical stress analysis, and even in this introductory treatment we wish to illustrate the power of the method. To outline a few necessary relations, consider z to be a complex number in Cartesian coordinates x and y or polar coordinates r and θ as z = x + iy = reiθ 4 (20) The French scientist Barré de Saint-Venant (1797–1886) is one of the great pioneers in mechanics of materials. 6 where i = the form √ −1. An analytic function f (z) is one whose derivatives depend on z only, and takes f (z) = α + iβ (21) where α and β are real functions of x and y. It is easily shown that α and β satisfy the Cauchy-Riemann equations: ∂α ∂β = ∂x ∂y ∂α ∂β =− ∂y ∂x (22) If the first of these is differentiated with respect to x and the second with respect to y, and the results added, we obtain ∂2α ∂2α + 2 ≡ ∇2 α = 0 ∂y ∂x2 (23) This is Laplace’s equation, and any function that satisfies this equation is termed a harmonic function. Equivalently, α could have been eliminated in favor of β to give ∇2 β = 0, so both the real and imaginary parts of any complex function provide solutions to Laplace’s equation. Now consider a function of the form xψ, where ψ is harmonic; it can be shown by direct differentiation that ∇4 (xψ) = 0 (24) i.e. any function of the form xψ, where ψ is harmonic, satisfies Eqn. 12, and many thus be used as a stress function. Similarly, it can be shown that yψ and (x2 +y 2 )ψ = r 2 ψ are also suitable, as is ψ itself. In general, a suitable stress function can be obtained from any two analytic functions ψ and χ according to φ = Re [(x − iy)ψ(z) + χ(z)] (25) where “Re” indicates the real part of the complex expression. The stresses corresponding to this function φ are obtained as σx + σy = 4 Re ψ  (z) σy − σx + 2 iτxy = 2 [zψ  (z) + χ (z)] (26) where the primes indicate differentiation with respect to z and the overbar indicates the conjugate function obtained by replacing i with −i; hence z = x − iy. Stresses around an elliptical hole In a development very important to the theory of fracture, Inglis5 used complex potential functions to extend Kirsch’s work to treat the stress field around a plate containing an elliptical rather than circular hole. This permits crack-like geometries to be treated by making the minor axis of the ellipse small. It is convenient to work in elliptical α, β coordinates, as shown in Fig. 4, defined as x = c cosh α cos β, y = c sinh α sin β 5 (27) C.E. Inglis, “Stresses in a Plate Due to the Presence of Cracks and Sharp Corners,” Transactions of the Institution of Naval Architects, Vol. 55, London, 1913, pp. 219–230. 7 Figure 4: Elliptical coordinates. where c is a constant. If β is eliminated this is seen in turn to be equivalent to x2 y2 + = c2 2 cosh α sinh2 α On the boundary of the ellipse α = α0 , so we can write c cosh α0 = a, (28) c sinh α0 = b (29) where a and b are constants. On the boundary, then x2 y 2 + 2 =1 (30) b a2 which is recognized as the Cartesian equation of an ellipse, with a and b being the major and minor radii . The elliptical coordinates can be written in terms of complex variables as z = c cosh ζ, ζ = α + iβ (31) As the boundary of the ellipse is traversed, α remains constant at α0 while β varies from 0 to 2π. Hence the stresses must be periodic in β with period 2π, while becoming equal to the far-field uniaxial stress σy = σ, σx = τxy = 0 far from the ellipse; Eqn. 26 then gives 4 Re ψ  (z) = σ + χ (z)] = σ 2[zψ  (z) ζ→∞ (32) These boundary conditions can be satisfied by potential functions in the forms 4ψ(z) = Ac cosh ζ + Bc sinh ζ 4χ(z) = Cc2 ζ + Dc2 cosh 2ζ + Ec2 sinh 2ζ where A, B, C, D, E are constants to be determined from the boundary conditions. When this is done the complex potentials are given as 4ψ(z) = σc[(1 + e2α0 ) sinh ζ − e2α0 cosh ζ] 8   π 1 4χ(z) = −σc (cosh 2α0 − cosh π)ζ + e2α0 − cosh 2 ζ − α0 − i 2 2 The stresses σx , σy , and τxy can be obtained by using these in Eqns. 26. However, the amount of labor in carrying out these substitutions isn’t to be sneezed at, and before computers were generally available the Inglis solution was of somewhat limited use in probing the nature of the stress field near crack tips. 2 Figure 5: Stress field in the vicinity of an elliptical hole, with uniaxial stress applied in ydirection. (a) Contours of σy , (b) Contours of σx . Figure 5 shows stress contours computed by Cook and Gordon6 from the Inglis equations. A strong stress concentration of the stress σy is noted at the periphery of the hole, as would be expected. The horizontal stress σx goes to zero at this same position, as it must to satisfy the boundary conditions there. Note however that σx exhibits a mild stress concentration (one fifth of that for σy , it turns out) a little distance away from the hole. If the material has planes of weakness along the y direction, for instance as between the fibrils in wood or many other biological structures, the stress σx could cause a split to open up in the y direction just ahead of the main crack. This would act to blunt and arrest the crack, and thus impart a measure of toughness to the material. This effect is sometimes called the Cook-Gordon toughening mechanism. The mathematics of the Inglis solution are simpler at the surface of the elliptical hole, since here the normal component σα must vanish. The tangential stress component can then be computed directly:  2α0 (σβ )α=α0 = σe sinh 2α0 (1 + e−2α0 ) −1 cosh 2α0 − cos 2β The greatest stress occurs at the end of the major axis (cos 2β = 1):   a (σβ )β=0,π = σy = σ 1 + 2 (33) b This can also be written in terms of the radius of curvature ρ at the tip of the major axis as   σy = σ 1 + 2 6 a ρ  J.E. Gordon, The Science of Structures and Materials, Scientific American Library, New York, 1988. 9 (34) This result is immediately useful: it is clear that large cracks are worse than small ones (the local stress increases with crack size a), and it is also obvious that sharp voids (decreasing ρ) are worse than rounded ones. Note also that the stress σy increases without limit as the crack becomes sharper (ρ → 0), so the concept of a stress concentration factor becomes difficult to use for very sharp cracks. When the major and minor axes of the ellipse are the same (b = a), the result becomes identical to that of the circular hole outlined earlier. Stresses near a sharp crack Figure 6: Sharp crack in an infinite sheet. The Inglis solution is difficult to apply, especially as the crack becomes sharp. A more tractable and now more widely used approach was developed by Westergaard7 , which treats a sharp crack of length 2a in a thin but infinitely wide sheet (see Fig. 6). The stresses that act perpendicularly to the crack free surfaces (the crack “flanks”) must be zero, while at distances far from the crack they must approach the far-field imposed stresses. Consider a harmonic function φ(z), with first and second derivatives φ (z) and φ (z), and first and second integrals φ(z) and φ(z). Westergaard constructed a stress function as Φ = Re φ(z) + y Im φ(z) (35) It can be shown directly that the stresses derived from this function satisfy the equilibrium, compatibility, and constitutive relations. The function φ(z) needed here is a harmonic function such that the stresses approach the far-field value of σ at infinity, but are zero at the crack flanks except at the crack tip where the stress becomes unbounded:  σy = σ, x → ±∞0, ∞, x = ±∞ −a < x < +a, y = 0 These conditions are satisfied by complex functions of the form σ φ(z) =  1 − a2 /z 2 7 (36) Westergaard, H.M., “Bearing Pressures and Cracks,” Transactions, Am. Soc. Mech. Engrs., Journal of Applied Mechanics, Vol. 5, p. 49, 1939. 10 This gives the needed singularity for z = ±a, and the other boundary conditions can be verified directly as well. The stresses are now found by suitable differentiations of the stress function; for instance ∂2Φ = Re φ(z) + y Im φ (z) ∂x2 In terms of the distance r from the crack tip, this becomes σy =  σy = σ  a θ θ 3θ 1 + sin sin · cos 2 2r 2 2  + ··· (37) where these are the initial terms of a series approximation. Near the crack tip, when r  a, we can write  (σy )y=0 = σ a K ≡√ 2r 2πr (38) √ √ where K = σ πa is the stress intensity factor, with units of Nm−3/2 or psi in. (The factor π seems redundant here since it appears to the same power in both the numerator and denominator, but it is usually included as written here for agreement with the older literature.) We will see in the Module on Fracture that the stress intensity factor is a commonly used measure of the driving force for crack propagation, and thus underlies much of modern fracture mechanics. The √ dependency of the stress on distance from the crack is singular, with a 1/ r dependency. The K factor scales the intensity of the overall stress distribution, with the stress always becoming unbounded as the crack tip is approached. Problems 1. Expand the governing equations (Eqns. 1—3) in two Cartesian dimensions. Identify the unknown functions. How many equations and unknowns are there? 2. Consider a thick-walled pressure vessel of inner radius ri and outer radius ro , subjected to an internal pressure pi and an external pressure po . Assume a trial solution for the radial displacement of the form u(r) = Ar + B/r; this relation can be shown to satisfy the governing equations for equilibrium, strain-displacement, and stress-strain governing equations. (a) Evaluate the constants A and B using the boundary conditions σr = −pi @ r = ri , σr = −po @ r = ro (b) Then show that   pi (ro /r)2 − 1 + po [(ro /ri )2 − (ro /r)2 ] σr (r ) = − (ro /ri )2 − 1 3. Justify the boundary conditions given in Eqns. 14 for stress in circular coordinates (σr , σθ , τxy appropriate to a uniaxially loaded plate containing a circular hole. 4. Show that the Airy function φ(x, y) defined by Eqns. 11 satisfies the equilibrium equations. 11 Prob. 2 5. Show that stress functions in the form of quadratic or cubic polynomials (φ = a2 x2 + b2 xy + c2 y 2 and φ = a3 x3 + b3 x2 y + c3 xy 2 + d3 y 3 ) automatically satisfy the governing relation ∇4 φ = 0. 6. Write the stresses σx , σy , τxy corresponding to the quadratic and cubic stress functions of the previous problem. 7. Choose the constants in the quadratic stress function of the previous two problems so as to represent (a) simple tension, (b) biaxial tension, and (c) pure shear of a rectangular plate. Prob. 7 8. Choose the constants in the cubic stress function of the previous problems so as to represent pure bending induced by couples applied to vertical sides of a rectangular plate. Prob. 8 9. Consider a cantilevered beam of rectangular cross section and width b = 1, loaded at the free end (x = 0) with a force P . At the free end, the boundary conditions on stress can be written σx = σy = 0, and  h/2 −h/2 τxy dy = P 12 The horizontal edges are not loaded, so we also have that τxy = 0 at y = ±h/2. (a) Show that these conditions are satisfied by a stress function of the form φ = b2 xy + d4 xy 3 (b) Evaluate the constants to show that the stresses can be written P xy σx = , I σy = 0, τxy P = 2I  2 h 2  −y 2 in agreement with the elementary theory of beam bending (Module 13). Prob. 9 13