3D Crack Analysis Using MSC.
Marc
Chris Timbrell, Gerry Cook, Ramesh Chandwani
Zentech International Limited, 103 Mytchett Road, Camberley,
Surrey, GU16 6ES, U.K.
http://www.zentech.co.uk
Introduction
This paper describes software which is interfaced to MSC.Marc for analysis of
arbitrary 3D cracks. The software, Zencrack, can be used in three ways:
1. For generating 3D finite element meshes containing multiple crack fronts
from a finite element model of an uncracked component.
2. For determining the distribution of the maximum energy release rates and
stress intensity factors along crack fronts.
3. For automatically calculating fatigue or sustained load crack growth in a
general 3D body under arbitrary loading.
Items 2 and 3 use the results from the LORENZI option of a MSC.Marc
analysis. The raw data from this option is processed to calculate energy release
rates along the crack front.
There are many complex issues to be addressed in 3D crack growth prediction.
This paper presents the key concepts of the software implementation and the
interface to MSC.Marc. A high level overview of the procedure is given in
Figure 1.
Mesh generation for 3D cracks
A schematic diagram of the steps taken to generate a cracked mesh is given in
Figure 2. The user must generate a finite element mesh for an uncracked
component. This must consist of 20 noded brick elements in the region(s)
where the crack(s) are required. An ancillary file is used to define the crack
related data.
A crack-block approach is used to introduce the crack front. The term crack-
block refers to a collection of 20 noded brick elements stored as a unit cube.
The arrangement of these crack-blocks is such that they contain either a quarter
circular or through crack front on one face. Part of this face is allowed to open
up under loading forming the opening crack face within the crack-block.
1
Examples of unit cube crack-blocks are shown in item 2 of Figure 2 with their
opening crack faces shaded for clarity.
The elements forming the crack front in each crack-block are modelled using
collapsed 20 noded brick elements, with an option to control the midside node
positions extending radially from the crack front to allow generation of quarter
point nodes if required. The crack front itself is seen as a line of nodes on the
crack-block surface. The internal mesh of the crack-block coarsens away from
the crack front such that some of the crack-block faces can be matched with a
standard 20 noded brick element. The other highly populated faces can be left
as free surfaces, symmetry surfaces, or connected to other compatible crack-
blocks.
The meshing procedure is one of replacement of one or more 20 noded brick
elements in the user supplied uncracked mesh by crack-blocks, as shown in
Figure 2. The elements to be replaced can be of arbitrary shape, but the faces
which become the initial crack plane should be as flat as possible to minimise
element distortion. During the mapping process to introduce the crack-blocks
the user can control the size of the generated crack front section for each crack-
block. Crack-blocks can be connected together to form distinct crack fronts of
the required size in the cracked mesh. In the first example of Figure 3 two
crack-blocks have been merged together to form a single semi-elliptic crack
front. This is done by replacement of elements 29 and 30 in the uncracked
mesh. The crack orientation is defined by key node positions 177, 162, 219.
The initial crack size is defined by ratios along the key edges defined by these
nodes.
In the simplest case a cracked mesh may contain only a single crack-block (and
therefore a single crack front). If there is a single crack-block then only one side
of the crack (i.e. one crack face) is modelled and symmetry constraints should
be applied. These can be applied in the uncracked mesh and automatically
updated by Zencrack. Figure 4 shows an example of update of symmetry
conditions from the uncracked mesh to the cracked mesh for a semi-elliptic
crack. The program can also update pressure loads in the crack region as
shown in the example of Figure 5. (Note that in the updated mesh each element
is listed in a separate entry on the DIST LOADS option since the face numbers
vary from one element to the next – this results in multiple boundary condition
sets when the cracked mesh is read by Mentat).
2
If both sides of the crack (i.e. both crack faces) are to be modelled, then pairs of
crack-blocks are used with a face-to-face match of the crack-blocks. This
approach also allows modelling of fully embedded cracks. Examples of
embedded cracks in a cube block are shown in Figure 6 (with one side of the
model removed to allow the cracks to be seen).
All merging of crack-blocks with one another and with the uncracked mesh is
carried out automatically to create a new finite element mesh containing the
required crack fronts. If requested, this new cracked mesh can be submitted for
f.e. analysis or a crack growth analysis can be carried out with this mesh as the
starter crack.
If the crack grows too large to be contained within the crack-blocks a restart
analysis can be carried out. This takes the history of crack growth calculated in
the analysis and continues growth with a newly modelled mesh containing
larger crack-blocks.
Crack growth prediction
Two parameters are generally used in L.E.F.M. (linear elastic fracture
mechanics) to describe conditions at a crack front:
KI KII KIII the stress intensity factors for mode I, II and III loading
respectively.
G the energy released at a crack front due to an increase in crack
surface area.
It is known that for mode I loading, K and G can be related by the material
properties (E and ν) and the state of stress :
G K2
=
1 − (α ν ) 2
E
where :
α=0 for plane stress
α=1 for plane strain
3
Crack growth is generally evaluated using one of the parameters K or G. For
fatigue crack growth the simplest equation relating loading and material data to
crack growth and fatigue cycles is the Paris equation of form:
da
=C∆Km
dn
where:
a = crack growth
n = number of fatigue cycles
C, m = material constants
∆K = stress intensity factor range
Generalising these expressions allows the crack growth rate to be written as a
function of ∆ K or ∆ G:
da da
=f(∆K) or =f(∆G)
dn dn
Zencrack always uses the second of these forms. The energy release rate G
can be accurately calculated using the finite element method. Since the energy
release rate is a fundamental quantity describing the formation of new crack
surfaces under any state of stress, the difficulties of assessing the state of
stress do not arise.
A Paris equation in terms ∆ K for mode I loading can be re-formulated in terms
of ∆ G. It is argued that when other modes of loading exist, the energy release
rate can always be used to describe crack growth from the mode I test data.
The direction of crack growth in a general case is assumed to be in the direction
of maximum energy release rate. The energy release rate approach has an
added advantage in that it automatically takes care of effects such as
anisotropy, plasticity, creep strains and other material non-linearities.
In order to take full account of such effects, the LORENZI extended J-integral
virtual extension technique is used in MSC.Marc.
Using the energy release rate distribution, the direction and magnitude of crack
growth is evaluated using an integration scheme and the crack growth law. This
integration is controlled by tolerances on the allowable “da” and “dn” increments
4
between one finite element analysis and the next. Figure 7 shows a typical
curve of growth against load cycles for a point on a crack front. It is clear from
this figure that early on in the analysis the increment should be controlled by
“dn” and that the control should switch to “da” as the crack extends.
An automatic crack advancement scheme in Zencrack updates the crack front
and generates a new finite element model containing the updated crack front.
Interface to MSC.Marc
Update of the uncracked model to incorporate the crack front(s) requires special
treatment of a number of MSC.Marc keywords:
SETNAME updated to allow for new crack related node and elements
sets
COORDINATES new nodes added
CONNECTIVITY new elements added, replaced elements removed
FIXED DISP update supports in crack region (e.g. Figure 4)
DIST LOADS update pressure loads in crack region (e.g. Figure 5)
A key feature for the cracked mesh is the LORENZI option. To allow analysis of
an arbitrary 3D crack the treatment of data for this option is non-trivial. Several
virtual crack extension directions are defined in order allow calculation of the
maximum energy release rate and its direction. This is shown schematically in
Figure 8. In addition, the raw values generated by the LORENZI option are local
nodal energy release rate values. These must be integrated along the crack
front paying attention to the element shape functions and the applied virtual
crack extension. The scheme developed by Zentech for this processing is
described in Ref. 1.
Zencrack recovers the results from the LORENZI option via the user subroutine
IMPD:
• During mesh generation Zencrack writes a file specifying the number of
nodes on the crack front and the number of LORENZI evaluations per node
• During the Marc analysis this file is read during the first call to IMPD in order
to establish the quantity of LORENZI output to be extracted from each
analysis increment
5
• LORENZI results data is extracted directly from the relevant common block
and written to a new file
• On completion of the Marc analysis all the LORENZI values are read from
the generated data file by Zencrack
Example - Compact Tension Specimen
An ASTM standard compact tension specimen has been modelled using
dimensions shown in Figure 9. A straight crack has been generated in the
model as can be seen in the mesh plots Figure 10 and Figure 11. The mesh is a
quarter symmetry model (half length and half thickness). Load is applied as
displacement control to the rigid surface load pin. A deformable contact body is
defined around part of the load hole. The load magnitude is recovered by
summation of reaction forces. The analysis provides an energy release rate
distribution which can be converted to stress intensity factor if an assumption is
made regarding stress conditions along the crack. The plane strain and plane
stress bounds are shown in Figure 12 compared with a reference solution (Ref.
2).
Example – Growth of an inclined crack in a plate
This analysis demonstrates crack growth prediction in a plate under cyclic axial
load. The plate contains an initial defect inclined at 43 degrees to the load
direction as shown in Figure 13. The full specification for this problem is given in
Ref. 3.
A displaced plot of the initial defect mesh is shown in Figure 14. Contour plots
showing the crack at two points during the growth prediction are shown in
Figure 15. As may be expected the crack very quickly turns to take up a mode I
configuration.
Comparison with experimental results is shown in Figure 16 and shows
excellent agreement. A typical plot of growth magnitude against number of load
cycles is shown in Figure 17.
6
References
Ref. 1 “Maximum Energy Release Rate Distribution From A Generalised
3D Virtual Crack Extension Method”, P.W.Claydon, Zentech
International Ltd., Engineering Fracture Mechanics, 42 (6), 1992
Ref. 2 Ewalds, H.L. & Wanhill, R.J.H., Fracture Mechanics, Edward
Arnold (Publishers) Limited, 3rd imprint 1986 (ISBN 0-7131-3515-
8)
Ref. 3 “Fatigue crack propagation in titanium under general in-plane
loading : I – Experiments”, Pustejovsky, M.A., Eng. Frac. Mech.,
Volume 11, pg 9-15, 1979
7
USER INPUT
an existing finite element mesh of
an uncracked component
(MSC.Marc .dat input file) ZENCRACK MSC.Marc
creates finite element mesh finite element analysis
USER INPUT of cracked component
supplementary data
e.g. crack location & size, crack ZENCRACK
growth law evaluates crack growth
ZENCRACK
updates crack front
geometry and f.e. model
ZENCRACK
stop no yes
next f.e. analysis ?
Figure 1 – Overview of crack growth scheme
2.
1.
Summary :
Five stage solution
for 3D crack modelling
and crack growth 3.
prediction 4.
1. Uncracked mesh
2. Selected crack-blocks
3. Mapped crack-blocks
4. Cracked mesh
5. Cracked mesh after crack growth 5.
Figure 2 – Example of meshing using crack-blocks
8
Element 29 30
Replaced by crack-block sq158x6.sup sq158x6.sup
Crack orientation via nodes 177 219 177 162
Initial crack size ratios 0.4 0.2 0.2 0.4
30
177
162
29
219
Element 26 25
Replaced by crack-block sq103x4.sup st88x5.sup
Crack orientation via nodes 157 105 157 162
Initial crack size ratios 0.4 0.3 0.3 0.3
26
105
162
157
25
Figure 3 – Definition of an initial crack
9
Figure 4 – Update of FIXED DISP data for a cracked mesh
10
Figure 5 – Example of update of DIST LOADS data for a cracked mesh
11
Figure 6 – Fully embedded cracks (one side shown)
12
Crack length, a
crack growth
“control box” “da”
“dn” B
Number of cycles, N
Figure 7 – Controls on crack growth increment
Local energy release rate, G
vs
Inclination to crack plane, θ (in the normal plane)
G
Gmax
G4 G5
G3 G6
G2
θ
θmax
G1 G7
Figure 8 – Distribution of energy release rate around a crack front
13
a=48.7, W=100
Figure 9 – ASTM standard compact tension specimen
14
Figure 10 – Compact tension specimen
15
Figure 11 – Compact tension specimen
16
9500
9000
Stress intensity factor
8500
Plane stress assumption
Plane strain assumption
Reference
8000
7500
7000
0 20 40 60 80 100
% from crack mid point
Figure 12 – Stress intensity factor for compact tension specimen
17
σ
h/W = 2.0 , a/W = 0.176
Stress ratio = 0.1
2h a Crack at 43° to load direction
2 distinct crack fronts
4 crack-blocks
Figure 13 – Plate with inclined crack
Figure 14 – Initial defect for 43 degree inclined crack model
18
Figure 15 – Crack model at two stages during crack growth prediction
19
Figure 16 – Comparison with experimental results
Surface node - growth vs cycles
10
9
8
Cummulative growth, a
7
6
5
4
3
2
1
0
0 5000 10000 15000 20000
Number of cycles, N
Figure 17 – Typical growth curve
20