The influence of interpolation errors on finite-element calculations involving stress-curvature proportionalities

Finite Elements in Analysis and Design 40 (2004) 1873 – 1884 www.elsevier.com/locate/!nel

The in$uence of interpolation errors on !nite-element calculations involving stress-curvature proportionalities A.E. Segalla;∗ , M.J. Sipicsb a

Engineering Science and Mechanics, The Pennsylvania State University, University Park, PA, 16802, USA b AGERE Systems, Allentown, PA, USA Received 3 April 2003; accepted 30 November 2003

Abstract The presence of non-linear axial gradients of pressure/temperature in a !nite-element model can invoke an often overlooked proportionality between the resulting curvature and bending stresses. Because these stresses can be signi!cant, the use of polynomials and cubic-splines to interpolate any gradients to a !nite-element mesh must be carefully weighed against their tendency to undulate through the data. As shown for a test case involving an interpolated pressure-distribution with arti!cially induced errors, the resulting polynomial oscillation can indeed induce signi!cant variations of both sign and magnitude in the !nite-element calculations. In contrast, a constrained B-spline with smoothing provided more reasonable stress predictions. ? 2004 Elsevier B.V. All rights reserved. Keywords: Polynomials; B-splines; Interpolation; Errors; Curvature; Bending stress; Axial gradients

1. Introduction There are many areas in applied mechanics where various order derivatives of experimental data are required either directly or indirectly. For instance, three-dimensional photoelastic studies of cracks [1,2] require the direct analysis of the !rst derivative of measured fringe-order data in order to determine the resulting stress state. Second-order derivatives, denoted herein as curvature, show up frequently in studies of beams and plates [3,4], thermoelasticity [5–8], and when pressure gradients exist in thick-walled tubes [9]. In these instances, the curvature of the deformation and underlying pressure and temperature gradients induce additional bending stresses. Concavity and the corresponding positive second-derivative cause compressive stresses on the exposed surface while ∗

Corresponding author. Tel.: +1-814-865-7829; fax: +1-814-865-9974. E-mail address: [email protected] (A.E. Segall).

0168-874X/$ - see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2003.11.006

1874

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

convexity switches the sign. In either case, the direct proportionality between the curvature and bending stress imposes a potential, albeit often overlooked complication and potential source for errors if not matched by the model. Although not a problem for a uniform !eld, this can become a problem when a varying temperature or pressure !eld must be described analytically by some form of curve-!t as is often the case for !nite-element models where the nodes and/or elements outnumber the points of measurement. Perhaps the most commonly used method of interpolation of a data set or function f(z) involves least-squares polynomials [10] f(z) =

N


an z n ;

(1)

n=1

where z is the axial coordinate and an are the polynomial coeEcients up to order N determined through a minimization of the sum of the squares of the residuals. Least-squares methods of interpolation oHer the inherent advantage of smoothing the data to some degree because the curve is not forced to pass through each data point. However, least-squares polynomials often require a relatively high number of terms to adequately describe complex data resulting in a curve that undulates around the original data points. While the curve itself may be adequate for interpolation, the changes in slope and curvature about the data points may have an increasingly deleterious eHect if an underlying proportionality to the higher-order derivatives exists. This is an inherent and often un-noticed drawback to traditional polynomials in that any errors in the data are signi!cantly magni!ed in their derivatives [11]. To counter these drawbacks, piece-wise polynomials such as the popular cubic-splines [12] are often used because of their inherent ability to model complex data using a limited number of terms. The so-called natural cubic-spline is a piece-wise polynomial with end conditions f(z) = an + bn (z − zn ) + cn (z − zn )2 + dn (z − zn )3 ;

zn 6 z 6 zn+1 ;

(2)

where the coeEcients an , bn , cn , and dn are determined for each of the data intervals with 0 6 n 6 m and m as the order of the polynomials (m = 4 for a cubic spline). For strictly interpolation purposes, cubic splines oHer many advantages in that the curve must pass through each data point regardless of the complexity of the data. However, the forcing of the curve through each data point has the potential to magnify any inherent errors in the data when derivatives are required. One obvious solution to all of the problems discussed so far is a combination of least-squares smoothing and the versatility of a piece-wise polynomial [13,14]. More recent advances using constrained B-spline formulations and a least-squares smoothing algorithm [15,16] add the ability to control the shape of the interpolating curve f(z) =

m


An Bn (z);

(3)

n=1

where each function, Bn is a spline (basis set) that is zero over as much of the interval as possible   0 z 6 zi − 3 (4) Bn (z) = A1 (z − zi−3 )3 + A2 (z − zi−2 )3 + A3 (z − zi−1 )3 + A4 (z − zi )3

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

1875

and An are interval constants. For a constrained B-spline representation of data, various combinations of the function f(z), concavity, convexity, and/or the sign of the slope can be imposed throughout the curve. Hence, the constrained B-spline and interpolations polynomial can be “forced” to reasonably approximate the data without many of the drawbacks previously discussed. Since a !nite-element implementation of a gradient will usually require some form of interpolation to assign values to individual nodes and/or elements, the type of interpolating curve chosen for the analysis can have a profound in$uence on the resulting stress calculations. In order to better understand the implications of these interpolation issues, a study of commonly used curve-!tting methods was undertaken. This paper explores the potential for both constrained and smoothed formulations and evaluates the advantages relative to traditional cubic splines and polynomials for !nite-element calculations where stress-curvature proportionalities exist. 2. Stress-curvature considerations Since the intent of the current study is to explore the potential implications of interpolation and !nite-element modeling, a practical example involving an axisymmetric pressure gradient was used. For the generalized case shown in Fig. 1, an axisymmetric pressure vessel or thick-walled pipe is subjected to axial pressure-gradients on the internal and external surfaces. As derived by an earlier analytical study [9], the ensuing stress-state can be determined using a truncated series expansion based in-part on Lame’s well-known analytical solution     1 A B 2  r = D(z) 2 − − S(z) +  D (z) B; r 2 − A2 + (A2 − 2 B1 ) 2 r 2 2 r  

  r  S  (z) a 2 b2 + r 2 − a2 + b 2 + 2 + ··· ; (5a) + A1 ln a 16(1 + ) r     1 −A B a2  2 2 −  = D(z) − (z) 3B r + (A − A ) − (A − a B ) S(z) +  D 1 1 2 2 1 r2 2 2 r2 

 r  S  (z) a 2 b2 2 2 2 3r − (a + b ) − 2 + + ··· ; (5b) + A1 ln 16(1 + ) a r

Fig. 1. Thick-walled tube with generalized pressure distributions acting on the internal and external surfaces.

1876

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

  r  S  (z) + z = D (z) 4B1 r 2 + (A1 − 2A2 ) + 2A1 ln a 8(1 + ) ×[2r 2 − (a2 + b2 )] + · · · ;

(5c)

where r is the radial stress,  the circumferential or hoop stress, z the axial stress, a is the inner radius, b is the outer radius, r is the radial coordinate,  is Poisson’s ration, and the primes represent diHerentiation with respect to the axial coordinate z. The diHerentiable quantities D(z) and S(z) represent the axial pressure distribution based on the inner, p(z) and outer, q(z) pressures, respectively D(z) = q(z) − p(z);

(6a)

S(z) = q(z) + p(z):

(6b)

For the current analysis, the outer pressure distribution will be assumed to be zero so that q(z) = 0. In addition, the recurring terms seen in the expressions are de!ned as follows: a 2 b2 ; b2 − a 2 A ; A1 = 2(1 − )

A0 =

A2 = (b2 + a2 )B1 ;   b b2 + a 2 A A1 ; + 2 ln B= 2 2 b −a a a B : B1 = 16(1 + )

(7a) (7b) (7c) (7d) (7e)

An examination of Eq. (5) reveals that the !rst terms represent Lame’s well-known solution for a thick-walled vessel under a known pressure diHerence, whereas the second term reveals the underlying proportionality to the curvature of the pressure pro!les. As long as the axial pressure pro!le is constant and/or linear such that D (z) = S  (z) = 0, the solution will be independent of the axial coordinate. However, if the pro!les are in fact non-linear, then curvature-induced bending-stresses will also be present in the solution. Moreover, if the interpolating curves used to model the pressure distributions p(z) and q(z) are not capable of accurately capturing the essence of the pressure distribution along with its curvature, potentially signi!cant polynomial errors can result. 3. Finite-element modeling For the comparison of interpolating methods, an axisymmetric !nite-element model of a thickwalled cylinder was !rst developed. Using the ANSYS (ANSYS Inc., Canonsburg, PA) generalpurpose code, a relatively !ne mesh was developed using quadrilateral isoparametric-elements with an inner radius of =8, an outer radius of =4, and a half-length of  units; by using the

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

1877

Fig. 2. Normalized axial stress for a sinusoidal pressure distribution on the internal surface of a thick-walled tube.

axisymmetric formulation and half-length symmetry, a uniform mesh of 800 elements (10 through the thickness and 80 along the axis) was found to be suEcient for modeling the stresses and maintaining an error norm well below 5% for all calculations. Materials properties were assumed to be that of steel with an elastic modulus of 208 GPa (30 Msi) with a Poisson’s ration of 0.3. As a precaution, the mesh was !rst tested under a uniform pressure and compared to Lame’s solution. Under this constant pressure scenario, the FEA results were within 1% of the closed-form solution. A second test was also performed using a sinusoidal pressure distribution p(z) = P sin(z + =2) acting on the entire -long internal surface. A sinusoidal function was chosen for the this portion of the study because the absolute values of the induced displacement (zeroth derivative) and curvature (second derivative) assume the same form of the original sinusoid and therefore, minimizes the need for the higher-order terms alluded to in Eq. (6). As shown by Figs. 2 and 3, excellent agreement was seen between theory and FEA predictions for the axial and hoop stress distributions. The discrepancy near the edge results from the fact that the analytical solution assumes an in!nite length while the FEA re$ects the !nite length and corresponding edge eHects. In order to investigate the in$uence of the various interpolation schemes for a more diEcult interpolation free from edge-eHects, a symmetric sinusoid acting only over a portion of the total length was used. Accordingly, the applied pressure distribution on the internal surface of the -long FEA model was assumed to vary as sin(z + =2) for a length of =2 while the outer pressure was held at zero such that  |D(z)| = |S(z)| = P sin(z + =2); 0 6 z 6 ; (8a) 2  6 z 6 ; (8b) |D(z)| = |S(z)| = 0; 2 q(z) = 0;

0 6 z 6 :

(8c)

1878

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

Fig. 3. Normalized circumferential (hoop) stress for a sinusoidal pressure distribution on the internal surface of a thick-walled tube.

Fig. 4. Sinusoidal pressure distribution acting on the internal surface of a thick-walled tube with added random errors.

In order to simulate the errors inherent in empirical data, a random number algorithm was used to generate absolute errors of ±0:05 in addition to relative errors of up to ±10% into the data generated by the sinusoid. Fig. 4 shows the analytical form of sin(z + =2) for the interval 0 6 z 6 =2 and the equally spaced data points complete with randomly induced errors.

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

1879

4. Results Once the error-laden data was generated, three diHerent curves were then used to !t the sin(z+=2) data as shown by Fig. 5. The !rst curve was a standard polynomial consisting of 6 integral order terms Dp (z) =

5


an z n

(9)

n=0

determined using standard least-squares methods. An alternative !t was also achieved by using a natural cubic-spline polynomial Ds (z) = an + bn (z − zn ) + cn (z − zn )2 + dn (z − zn )3

(10)

with the end derivatives set to zero and −1 for the data intervals zn 6 z 6 zn+1 and 0 6 n 6 4, respectively. As discussed earlier, this type of curve is commonly used and can be readily found in many commercially available software packages and numerical libraries. Finally, a constrained B-spline representation with least-squares smoothing [16] was also used to interpolate the data B

D (z) =

m


An Bn (z):

(11)

n=1

Fig. 5. Sinusoidal pressure distribution with added random errors, approximations based on various interpolation algorithms, and the resulting curve instability or “undulate.”

1880

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

Fig. 6. First derivatives of the imposed sinusoidal pressure distribution, approximations based on various interpolation algorithms, and the magni!ed curve instability or “undulate.”

For the B-spline curve shown in Fig. 5, convexity was maintained throughout the curve by imposing the following constraint: d2 DB (z)  ¡ 0; 0 6 z 6 (12a) dz 2 2 with the additional symmetry condition dDB (z) = 0; z = 0 (12b) dz required to force the curve maxima to occur at the center line (z = 0). As shown in Fig. 5, all curves were found to be reasonably capable of describing the function sin(z + =2) and the underlying data used to de!ne the curve. However, the plot of their !rst and second derivatives shown in Figs. 6 and 7 demonstrates how the curves diHer in terms of the predicted slope and curvature. The relatively large $uctuations in the slope and curvature of the polynomial, Dp (z) and cubic spline Ds (z) !t are caused by the inherent undulation as the curves are forced to be near (polynomial) or through (cubic spline) the data points. In contrast, the B-spline must only maintain the proper convexity as the curve is allowed to pass near, but not necessarily through the data points. In terms of the resulting stress components, the comparative FEA studies revealed a number of interesting !ndings pertaining to the potential interpolation errors. For all of the curve !ts used in this study, the hoop stresses show reasonable agreement with theory as indicated in Fig. 8. The relative unimportance of the interpolating polynomial error is because the curvature only manifests itself through a Poisson’s eHect and is therefore of a lesser degree. On the other hand, the type of interpolating curve did have a signi!cant in$uence on the axial stresses predicted by the model. As shown in Fig. 9, the inherent undulation around the data points for the least-squares polynomial has produced spurious axial stress distributions on the internal surface of the vessel. The magnitude

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884 Constrained Natural Polynomial Sin(z+π/2)

2

2

dz

10

d f (z)

1881

0

-10 0

0.25

0.5

z/L

Fig. 7. Second derivatives of the imposed sinusoidal pressure distribution, approximations based on various interpolation algorithms, and the magni!ed curve instability or “undulate.”

Fig. 8. Normalized circumferential (hoop) stress for a partial sinusoidal pressure distribution on the internal surface resulting from various interpolation algorithms.

of the errors in both sign and magnitude become evident when the stress state is compared to the theoretical predictions. While the use of relatively low-order polynomial has produced a reasonably good curve-!t, the resulting !nite-element implementation is clearly lacking when curvature eHects become dominant.

1882

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

Fig. 9. Normalized axial stress for a partial sinusoidal pressure distribution on the internal surface resulting from various interpolation algorithms.

Similar results were seen with the natural cubic spline polynomial as also shown in Fig. 9. For the natural cubic-spline, the inherent undulation around the data points has again produced spurious axial stress distributions on the internal surface of the vessel. As with the least-squares polynomial just examined, the magnitude of the errors in both sign and magnitude become evident when the stress state is also compared to the theoretical solution. Because the spline curve is forced to pass through each individual data point, the largest errors are concentrated there. Somewhat more reasonable approximations are made for the interval between the data points. However, the overall stress distribution must be considered erroneous and would ultimately produce poor design information. As anticipated based on the smoothed and controlled curvature, the best results were seen with the constrained B-spline as shown in Fig. 9. Because of the smoothing of the least-squares algorithms, the inherent polynomial error around the data points that can result in spurious stress distributions has been greatly reduced. The spurious oscillations of stress magnitude and sign have been mostly eliminated with the stress distribution assuming a reasonable form. While the resulting stresses still show a small tendency to oscillate around the early data points, the overall distribution is relatively smooth. However, there are two additional artifacts of the FEA and theoretical solution that warrant further discussion. Firstly, the theoretical axial stresses do not always show good agreement with the FEA predictions except in the interval up to =4. This is primarily due to the pressure and curvature in$ection caused by the abrupt end of the sinusoidal pressure at =2. It is this type of behavior that shows the limitations of the closed-form solution and highlights the advantages of FEA provided the curvature is adequately interpolated. Otherwise, the errors due to interpolation can easily oHset any advantages of the FEA method if care is not taken. Another important and related point is that the

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

1883

pressure values were averaged to individual element faces along the inner radius of the tube model. Hence, the implementation of the pressure to the element faces and not individual nodes may have inadvertently smoothed the pressure distributions to some degree. Similar smoothing would not be expected for thermal solutions because the temperatures are imposed directly to the nodes de!ning each surface element. Hence, the interpolation errors examined in this study could be signi!cantly worse for thermal stress models. 5. Conclusions An investigation was undertaken on the implications of using various interpolation schemes for !nite-element modeling of curvature-induced stresses. When compared to theoretical predictions for a thick-walled vessel under a non-linear axial pressure gradient, a constrained B-spline representation with least-squares smoothing was found to produce the most accurate results. In contrast, both a least-squares polynomial and a natural cubic spline were found to produce relatively large variations of both the magnitude and sign of the predicted stresses. These variations can be attributed to curve oscillation as the curves were forced near (polynomial) or through (spline) the data points. While these curve oscillations may appear to be relatively small, they can cause signi!cant errors when the bending stresses are proportional to the curvature of the modeled pro!le. Since interpolation is often required for a !nite-element implementation of pressure and temperature !elds, signi!cant errors can result. Hence, the type of curve and the degree of control employed by the numerical algorithm can have a profound in$uence on the accuracy of the !nite-element predictions.

References [1] A.E. Segall, J.C. Conway Jr., Scattered-light analysis of a surface-$awed plate subjected to cylindrical bending, Exp. Mech. 27 (1987) 64–67. [2] A.E. Segall, J.C. Conway Jr., Scattered-light determination of mode I stress-intensity factors by a fringe gradient independent technique, Exp. Mech. 30 (1990) 173–176. [3] S. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, NY, 1951. [4] C.W. Lee, Thermal stresses in a thick plate, Int. J. Solids Struct. 6 (1970) 605–615. [5] A.E. Segall, J.R. Hellmann, Analysis of gas-!red ceramic radiant tubes during transient heating: part II—thermoelastic stress analysis, ASTM J. Testing Eval. 20 (1992) 25–32. [6] A.E. Segall, J.R. Hellmann, R.E. Tressler, Thermal shock and fatigue behavior of ceramic tubes, Proceedings of the 10th Biennial ASME Conference on Reliability, Stress Analysis, and Failure Prevention, 1993, pp. 81–91. [7] C.W. Lee, Thermoelastic stresses in thick-walled cylinders under axial temperature gradient, J. Appl. Mech. 88 (1966) 467–469. [8] K.W. Yang, C.W. Lee, Thermal stresses in thick-walled circular cylinders under axisymmetric temperature distribution, ASME paper 71-PVP-16, 1971. [9] C.W. Lee, A theory for thick-walled cylinders under axisymmetric loading, Proceedings of the Fourth National Congress on Applied Mechanics, ASME, New York, NY, 1962, pp. 667–675. [10] T.J. Riolin, An Introduction to the Approximation of Functions, Blaisdell, Waltham, MA, 1969. [11] S. Whitaker, R.L. Pigford, Numerical diHerentiation of experimental data, Ind. Eng. Chem. 52 (1960) 185–187. [12] J.L. Walsh, J.H. Ahlberg, E.H. Nilson, Best approximate properties of the spline !t, J. Math. Mech. 11 (1962) 225–234.

1884

A.E. Segall, M.J. Sipics / Finite Elements in Analysis and Design 40 (2004) 1873 – 1884

[13] D.G. Berghaus, J.P. Cannon, Obtaining derivations from experimental data using smoothed-spline functions, Exp. Mech. 13 (1973) 38–42. [14] R.E. Rowlands, T. Libor, I.M. Daniel, P.G. Rose, Higher order numerical diHerentiations of experimental information, Exp. Mech. 14 (1973) 105–112. [15] D.E. Amos, Computation with Splines and B-Splines, Sandia Report No. SAND78-1968, 1978. [16] R.J. Hanson, Constrained least-squares curve !tting to discrete data using B-splines: A Users Guide, Sandia Report No. SAND78-1291, 1978.