A COMPARISON OF FINITE ELEMENT CODES

AND
RECOMMENDED INVESTIGATION METHODOLOGY

Michael A. Porter Dennis H. Martens

Dynamic Analysis Black & Veatch Pritchard, Inc.
Leawood, Kansas Overland Park, Kansas

Charles S. Hsieh
Black & Veatch Pritchard, Inc.
Overland Park, Kansas

ABSTRACT
Significantly different results attained from the use of three Finite
Element codes used in the analysis of a large complex model are
discussed. Building on previous work by the authors regarding the
comparison of stress results from several commercial FE codes used
on a simple model, this paper recommends steps for an investigation Shell (Brick)
methodology to aid in ascertaining results which are most
representative, useful and correct.
Tubes (Plate)
DESCRIPTION OF PROBLEM
The current commercial finite element (FE) codes represent an
effective analysis tool for the investigation of critical aspects of
pressure containment equipment. The need to formulate an adequate
model, apply the correct loadings and confirm the results represents a Knuckle (Brick)
major portion of the engineer’s analysis task. Once a valid model is
constructed, the variability of results produced by the different FE Tubesheet (Brick)
codes remains a concern. To evaluate the magnitude of these
variations, we will examine a relatively large and complex model. This Figure 1 - Pictorial View of Heat Exchanger Model
model will be “solved” using three commercial FE codes. The results,
in terms of the indicated Stress Intensity, will be compared to FEA MODEL
investigate the variability of the results. The model used for this project was that of a shell and tube heat
For this comparison, the same model was used with each FE exchanger. Figure 1 illustrates the model pictorially. The various
code, except for some minor variations required to adapt to the components of the exchanger model are indicated along with the type
specific code. The results derived from the three codes were compared of element used in the model. Not indicated are the beam elements
to closed-form solutions for pressure and thermal displacements. The that were used to attach the tubes to the tubesheet. This form of
variability of the indicated displacements and stresses were then connection resulted in some high indicated stresses at the attachment
compared to assess the ability of the commercial codes to give points. This modeling technique was employed as a means of limiting
reasonably accurate and consistent results. This methodology of the model size. Since prior work had shown that the stress in the tubes
establishing a valid model of the problem, verifying the model by was not a primary concern, the stress at these junctures was not
closed-form solutions, and comparing and interpreting the results evaluated.
forms a basic protocol that can be extended to other FE analysis work. Figure 2 illustrates a cross-section through the knuckle portion of
the model. As may be seen, the material thickness in the knuckle
portion of the vessel was less than the shell thickness and greater than
the tubesheet thickness. Although a formal mesh convergence study

1997 ASME
 contribute to the axial movement of the tubesheet. Thus, we would
Shell expect that the actual extension will be less than computed above.
When the FE results of the two models were compared with the
values computed above, the displacements in the code B model were
found to be almost an order of magnitude greater. Further
examination of this model indicated that the Young’s Modulus (E)
value that had been entered in code B (by the author, unfortunately)
was low by an order of magnitude. This value was corrected and the
Knuckle displacements and stresses were re-computed. The results from codes
(radius) A and B using the correct E value were in close agreement, as will be
discussed in the next section.
The use of the computed displacements facilitated identification
of the incorrect exponent on the material E value. Reliance on a check
Tubesheet of the hoop stress alone, as is often done when checking vessel
analyses, would not have identified the problem. As may be seen in
Appendix Equation 1, the hoop stress is independent of the material E
value. The hoop stress in the model with the low E value was
substantially the same as in the model with the correct E value! The
maximum stress in the knuckle, however, was low in the low E model
Figure 2 - Cross Section Through Knuckle Region of Model by nearly 30%. Thus, even though the hoop stress may “look” right in
a model, this does not necessarily mean that the model is correct.
was not conducted, several iterations of this geometry were employed More checks must be conducted. Checking the displacements is a key
to assure that the elements in this area were reasonably behaved. place to start.
The completed model has approximately 11,390 nodes that define In order to avoid further “operator entry error” problems between
some 8,500 elements. Of this total, approximately 4,500 elements the models, the program FEMAP (1997) was used to translate one of
were the brick elements used to model the shell, knuckle and tubesheet the models into the native format for each of the other two codes.
portions of the model. Symmetrical constraints were applied to the cut Following the solutions phase, the model results were examined using
surfaces indicated in the model, so that only 1/8 of the exchanger (as the FEMAP post-processing facilities.
illustrated) was required for the analysis. The assembled model had
approximately 41,000 degrees of freedom. The solution times ranged
from about 3 to 6 hours on Pentium 90-100 computers, depending on
the code used.

COMPUTATION OF DISPLACEMENTS AND STRESSES

The original model was constructed using one of the popular PC
based FE packages (code A). As a means of checking the results, a
geometrically similar model (using most of the geometry of the
original model) was constructed using another FE package (code B).
The results (stress) obtained from the two models differed by
approximately 30% in the critical knuckle region! In order to try to
resolve the differences, several checks of the results were performed
using closed-form “hand” approximations.
In the appendix, equations used to approximate the hoop stress,
radial expansion and axial expansion are presented. All these
equations may be found in Roarks’s Formulas for Stress and Strain.
Note that such a reference, in the authors’ opinion, is a mandatory tool Figure 3 - Indicated Stress Intensity
for anyone conducting FE analyses.
(Displaced Shape - Pressure Only)
For this vessel, the computed hoop stress using the equations in
the appendix is approximately 16,200 psi. The computed radial
displacement is approximately 0.035” and the computed axial INTERPRETATION OF RESULTS
extension for the tubes and tubesheet is approximately 0.026”. Figure 3 illustrates the deflected shape of the shell portion of the
However, the above analysis ignores the restraint that the shell will vessel along with the indicated stress intensity as computed by one of
the packages. As may be seen, there is a reverse bending that takes

1997 ASME
place in the knuckle region. The highest stresses are indicated in the of the knuckle, as a function of the distance in inches from the point of
region where this bending is at a maximum. The indicated hoop stress highest indicated stress, was plotted.
in the vessel away from the tubesheet was approximately 16,000 psi
independent of the FE code used to compute the deflections and
stresses. In addition, the region where the highest stress was indicated
was approximately the same for all three FE codes examined. Figure 4
illustrates the deflected shape along with the magnitude of the
displacements. All three FE codes indicated a radial displacement of
approximately 0.036” and an axial displacement of approximately
A
0.024”, indicating good agreement with the hand computations. Thus,
on initial examination, the same model run on different processors
seems to produce consistent displacement results.

Figure 5 - Location of Stress Reporting Line

Figure 5 illustrates a section of the model through the knuckle

near the region of highest indicated stress. The line A-B passes
through the point of the highest indicated stress in the model. In
Figure 6, we illustrated the stresses along this line (at nodal points)
indicated by the three FE codes. The horizontal axis in Figure 6 is the
distance away from the point of highest stress. As may be seen, the
stress intensities indicated by two of the codes appear to be in close
agreement, while those indicated by the third code (code C) are
Figure 4 - Magnitude of Indicated Displacements considerably lower.
(Displaced Shape - Pressure Only)
In s id e S tre s s In te n s ity

However, if we look at the magnitude of the maximum indicated 30000

Sy = 29.1 ksi @ 600 F
o

Stress Intensity (in the knuckle region), we find that the agreement is
less than satisfactory. As illustrated in Table 1, the results from one of 25000
1.5*Sm

the codes differs from the others by nearly 30%. Code A

Sm

Code C' Code B

Stress Intensity

20000
Table 1
Maximum Indicated Values
15000
Code C

Code Stress Radial Axial

Intensity Displacement Displacement 10000
A 25,963 0.036 0.026
B 25,760 0.036 0.026
C 19,137 0.035 0.024 5000
-6 -4 -2 0 2 4 6
Distance From Highest Stress Point - Inches

In an earlier paper, Porter and Martens (1996) illustrated that, Figure 6 Stress Intensity Contours Through Point of
while the maximum indicated Stress Intensity on a model comprised of Highest Stress
all plate elements might differ significantly depending on the FE code
used, the agreement away from this point of maximum value was quite Code C, as indicated by the lower curve, uses an approximation
good. As a check to see if such agreement could be found for the recommended by Cook (1981) to avoid the stiffening effect on linear
brick elements being used in this case, the stress on the inside surface isoparametric elements caused by parasitic shear in bending. When
this approximation in code C is replaced with an incompatible mode

1997 ASME
formulation (see curve Code C’ on Figure 6), all three codes report Table 2
nearly the same value. Note that testing the different element Results of “Bar” Model
formulations in code C was facilitated by code C having the shortest
run time. We would expect similar results with the other codes. In the Code Maximum Indicated Deflection
following section, assuming that we had only the original (lower Stress
curve, code C) values for code C, we will discuss the results. A 21,855 0.496
B 21,803 0.496
C 21,704 0.496
ATTEMPTED RESOLUTION OF DIFFERENCES
The nearly 30% difference between the highest stress intensity As may be seen, the differences between the results indicated by
indicated by two of the codes and that indicated by the third code is, the three codes are quite small and in good agreement with the closed
without question, too great a difference to ignore. The conservative, form solution. Unlike the heat exchanger model, in the Bar model
“safe” assumption would be that there is something wrong with code C there is no large difference between the stress reported by code C and
or with the element formulation, and that the values reported by codes the other two codes. Thus, we have not identified a difference in the
A and B are correct. The possibility exists, however, that code C is, in codes that would explain the differences reported in the heat
fact, reporting the correct value and that the use of codes A and B exchanger models (except that the models only agree when the same
would result in an overly conservative design. In order to attempt to element formulation is used). If an axial load is added to the Bar
resolve this matter, an additional verification model was constructed. model to more closely approximate the loading of the elements in the
knuckle region of the heat exchanger model, the indicated stresses get
closer rather than farther apart, even though different element
formulations are being used.
The indicated differences in the reported stresses based on
element formulation is an unresolved issue. As users, we need better
guidance on which element formulations to use and when to use them.
It obviously makes a difference. But, as may be seen from the bar
problem, the conditions under which various formulations may affect
the answer is not always clear.

CONCLUSIONS
Modeling technique is very important in the development and use
of FE models. Without a carefully considered, constructed and tested
model, the results are questionable at best. As illustrated by the
unintentional use of a low E value in the original code B heat
exchanger model, the stresses can “seem” to be “OK,” but in fact they
can be low by a significant degree. Careful model construction and
model verification are essential.
Figure 7 - Bar Model Even with carefully constructed and tested models, however, the
answer is still dependent on the code and element formulation used for
Figure 7 illustrates the “Bar” model that was constructed to test the solution. The vessel design using codes A and B would be more
the brick elements in the three codes. This model has dimensions of conservative than that using code C (and a different element
40” long by 6” wide and 1.5” thick. The elements are meshed so that formulation). At this point, there is no conclusive proof that the
the aspect ratios of the elements are approximately the same as those results from any of the codes (or formulation) are actually correct. All
used in the heat exchanger model. The bar is fixed at one end and three codes seem to give the same answer when examining a simple
loaded with 1,200 LB on the other end, as illustrated. Using the problem (e.g. Bar) for which there is a handy closed-form solution.
familiar Mc/I equation for the maximum stress at the point of The loading and element shapes seem to be approximately the same in
attachment, we get a computed stress of 21,331 psi. The maximum both the simple model and the more complex model, yet the stress
deflection at the end, computed from d = Fl3/3EI is 0.505”. Table 2 intensity results differ dramatically for the more complex model.
illustrates the stresses and deflections computed using the three FE At this point, the phrase “caveat emptor” (Let the Buyer Beware)
codes. would seem to apply to the use of commercial FE codes. Testing
(using a simple model against known results) seems to indicate that all
three of the codes examined give the same results. However, when
used in an actual vessel analysis, the results differ dramatically. If
code C (with the first element formulation) is reporting stresses that

1997 ASME
are 30% low, then a design at 1.5 Sm (ASME, 1996) would likely would have verifiable solutions and be made available to other
result in stresses that exceed yield. engineers to test the codes that they are using. PVP should consider
How many of us really know the performance level of the codes taking the lead in the development of such a suite.
that we and/or our vendors are using? And with the codes that we use,
do most of us really know which element formulation to use for each REFERENCES
particular application?
ASME Boiler and Pressure Vessel Code Section VIII Divisions 1 and
2, 1996, The American Association of Mechanical Engineers, New
RECOMMENDATIONS York, NY
The authors recommend that the engineer/user carefully select an
FE code that gives consistent and accurate results for the type of Cook, Robert D., 1981, Concepts and Applications of Finite Element
problem being investigated. Additionally, the engineer using the code Analysis, John Whiley & Sons, Inc., New York, NY, pp 190-195
needs to be aware of the element options available and the
consequences of their use. The variability of FE code results can be FEMAP, 1997, Enterprise Software Products, Inc., 415 Eagleview
significant and, in many cases very difficult to confirm. The engineer Boulevard, Exton, PA. 19341
must be cautious, think through and verify the results of an FE
analysis before using them as a basis for a design. Such checking of Porter, M. A. and Martens, D. H., 1996, “A Comparison of the Stress
the FE solution is essential to the achievement of a practical and safe Results from Several Commercial Finite Codes with ASME Section
result. VIII, Division 2 Requirements,” PVP Vol. 336, ASME, pp 361-371
The authors recommend that the following methodology be used
for FE analysis: “Roark’s Formulas for Stress and Strain,” Sixth Edition, 1989,
McGraw-Hill Book Company, N.Y. pp 518-519
1. Select the FE code that you will be using based on its proven
applicability
2. Determine the design parameters that must be modeled to assure a
valid analysis. Consider:
• pressure loadings that must be applied to model
elements to simulate actual conditions
• thermal gradients and profiles that can be expected
during various operating condition
3. Select the type of elements that will react correctly to the above
design parameters and give reasonable results on the first attempt.
Convergence of critical stressed areas is required.
4. Confirm the model displacements for pressure and thermal
parameters.
• use closed-form solutions for pressure induced strain
• use closed-form solutions for thermal induced
expansion
5. Question the results for all critical displacement and stress areas
of the model. Confirm that the results appear to be intuitively
correct.
• movements are in the correct direction and magnitude
• stresses on individual elements of the model are
distributed adequately
Finally, as users of FE codes, we need to have better information
about the validity of the codes we select. All three of the codes
examined in this paper are advertised in the pages of Mechanical
Engineering magazine as being suitable for the type of analysis
conducted in this study. Although we have not tested other codes
using this problem, we suspect that the scatter of answers might
increase with the number of codes tested.
The authors recommend that the engineer using FE for design
analysis develop a suite of real life application problems for
verification of the FE codes being used. Ideally, these problems

1997 ASME
 APPENDIX - Check Calculations

Hoop Stress
The hoop stress in the shell, away from the ends, should be
approximately equal to:

Sh = Pr/t Equation 1

where: P = Pressure (790 psi)

r = radius of vessel (59”)
t = vessel wall thickness (2.875”)

Radial Expansion
The expansion of the vessel in the radial direction may be related to
the hoop stress by:

dr = Sh * r/E Equation 2

where: E = Young’s Modulus for the vessel material

(29 x 106 psi)

Axial Extension
If we assume that the axial extension of the vessel proportional to the
pressure load acting on the tubesheet and the axial stiffness of the
tubes, we may compute this extension using:

F = P*(Ats -Atb) Equation 3

Sa = F/Atb Equation 4

e = Sa/E Equation 5

dl = e * l Equation 6

where: F = Pressure force on tubesheet

Ats = Area of tubesheet ( 10,936 sqin)
Atb = Area of tube bundle (1,268 sqin)
e = Strain in tubes (in/in)
l = Length of tunes (126”)
dl = Axial extension of tubes and tubesheet

1997 ASME