`***compute_G_by_gth`#

Description#

This command is used to calculate the energy release rate G along a 3D crack front using the the energetic G-\(\theta\) method. Various options may be activated to:

smooth the crack front by defining the number of points of the spline curve used to build the variational problem solved to obtain the energy release rate (G) and crack virtual extension (\(\theta\)) at each control point of the underlying spline,
output results in terms of conventional stress-intensity factors (SIF) \(K_I\), \(K_{II}\), \(K_{III}\) instead of the default \(G\) values,
define propagation laws based on the previous SIF calculations, that will generate crack advance values that may be used to drive crack propagation by means of the ***auto_remesh command,
define various crack bifurcation criteria for out-of-plane crack propagations.

Note that this command may be applied either to conforming cracks FE models (ie. the crack is explicitly introduced in the mesh, and the crack front is defined by a liset) or xfem ones (in this case the discontinuity is defined by means of level-sets in a ***xfem_crack_mode block of commands, see ). Note also, that the input for this command may be built automatically by the interactive Zcracks and Zxfem commands.

Syntax#

The syntax is as follows:

***compute_G_by_gth component-name [ **crack_front liset-name ] | [ **xfem ] \(~\,\) **nbnodes nbnodes [ **theta ra [ rb ] ] | [ **elem_radius nbe ] [ **elset elset_name ] [ **exclude_BE ] [ **lip lipname ] [ **compute_K ] [ **compute_Ki ] [ **fatigue preload period ] [ **behavior btype ] [ **Delta_N deltan ] [ **h hmax ] [ **optim_dir prec scan ] [ **vectorial ]

The different option-keys are

liset-name is the liset (line set of nodes) defining the crack front for conforming mesh modelizations. This liset is usually automatically created in the FE mesh by the Zcracks interactive script, and maintained during propagation analysis.
alternatively the command **xfem should be used for xfem models.
as defined on the next figure, the parameter nbnodes (**nbnodes command) defines the number of points on the spline curve built internally to smooth out the crack front. Note that a negative values can be given for the nbnodes parameter, in which case the number of spline control points is calculated as a fraction of the total number of crack front nodes (number of nodes in liset liset-name in the conforming mode).
the **theta command defines the size of the integration domain used to build the G-\(\theta\) variational problem. As shown on the figure, this domain has the shape of a ring, centered around the crack front, where ra is the radius of the circular section of the ring. The optional parameter rb may also be entered, and correspond to the size of an inner volume, close to the crack front: in this case, the mechanical fields on the first rows of integration points/elements near the front are not considered precise enough, and are eliminated from the integration domain. Note that in practice, the outer radius ra of the integration volume, should be chosen to stay within the volume (ie. it should not reach the outer free surface) to avoid erroneous \(G-\theta\) solutions.
the command **elem_radius is an alternative way to specify the size of the integration domain. In this case the outer radius ra is calculated to correspond to nbe times the average size of crack front elements.
the optional command **elset can be used to speed up the G-\(\theta\) volume definition, ie. only candidate elements in elset elset_name will be tested for addition in the volume integration.
in general, values near free surfaces (at both ends of opened crack fronts) are less precise than inside the volume along the crack front. The optional command **exclude_BE may be used to reset those values, and calculate SIFs at the front ends by extrapolation of values obtained on inner points.
command *lip is not used in the xfem mode. lipname is a name of an nset in the FE mesh, containing nodes on the crack lips. In the conforming mode this nset is used to calculate the tangential plane at each crack front point, an the orientation of the front normals at this point.
when the optional **compute_K command is specified, the default \(G\) energy release rate output, is replaced by the mode I \(K_I\) stress intensity factor. This \(K_I\) value is obtained from \(G\) using the following equation:

(51)#\[K_I~=~\sqrt{\frac{E~G}{1-\nu^2}}\]

Note

In order to evaluate the previous equation, the \(G-\theta\) module needs to retrieve the elasticity coefficients of elements close to the current crack front node, and that the influence of external parameters (eg. temperature) on those parameters is accounted for. Note also, that when this option is used, the propagation law behavior automatically takes \(K\) values as input, instead of the default \(G\) energy release rate, such that coefficients of those laws should be modified accordingly.
the optional **compute_Ki command activates the calculation of the \(K_I\), \(K_{II}\), \(K_{II}\) stress intensity factors by means of an interaction integral computation based on the Westergaard analytical solution. This choice has an impact on the branching criterion for out of plane propagation, that automatically switches to a modal mixity angle \(\alpha\) calculated from those stress intensity factors instead of the default \(Gmax\) based criterion. In this case, \(\alpha\) is computed using the following equation:

(52)#\[\alpha ~=~2~ arctan \left( \frac{-K_{II}/K_I + \sqrt{ \left(K_{II}/K_I\right)^2 + 8}}{4} \right)\]
command **fatigue preload period is mandatory for fatigue crack propagation, and defines the cycle period (parameter period) and a preloading initial time value at the beginning of the calculation during which crack propagation is ignored (parameter preload). Depending on the above definitions, crack propagation will occur at times \(t~=~\)preload\(~+~n_{cyc}~\)period, where \(n_{cyc}\) is the cycle number.
command **behavior btype defines the model btype selected to compute crack advance, that will drive crack propagation and automatic remeshing procedures (by means of the ***auto_remesh command). Currently the only model available is the Paris law for fatigue crack propagation:

(53)#\[\theta ~=~ \frac{da}{dN} ~=~ C \left( \Delta G \right)^m\]

where \(\Delta G\) is the \(G\) variation during the cycle (or \(K\) if the **compute_K keyword is activated) computed at each point along the crack front.
command **Delta_N dn defines a multiplication factor dn applied to the crack advance \(\theta\) calculated with the propagation law specified by the **behavior command. This option may be important to insure that a significant crack advance will occur at each simulated cycle in the FE analysis, and that automatic remeshing is indeed pertinent.
command **h hmax is another way to control the value of the multiplication factor dn applied to the \(\theta\) values. When a positive value is given, hmax is the target crack advance when remeshing occurs, and the multiplication dn is calculated accordingly:

(54)#\[dn~~~such~that~~~hmax~=~\displaystyle \max_{i} \left\{ dn~\frac{\theta_i}{\theta_{max}} \right\}\]

where \(\theta_i\) is the crack advance at point \(i\) on the crack front, and \(\theta_{max}\) the max value of \(\theta_i\) on all crack front points. When a negative values is specified, hmax defines a target crack advance of hmax times the mean size of elements along the crack front.
the command **optim_dir prec scan triggers out-of-plane propagation. Without this keyword, the crack advance direction always lie in a plane corresponding to the initial crack definition: in the xfem mode this plane is characterized by the first levelset definition, and in the conforming mode the **lip nset allows this tangential plane calculation at each crack front point.

When a \(G_{max}\) bifurcation criterion is used (which happens when no keyword **compute_Ki has been specified, see above), 2 methods are available to compute the angle \(\alpha_{max}\) (at each crack front point) for which \(G=G_{max}\).

First one (default case) is a dichotomy search where parameter prec is the precision, and scan the angle variation used during the scanning of angle values needed to initiate the dichotomy algorithm. Angle values should be entered in degrees, and a typical definition is:

**optim_dir 1. 15.

for a precision of 1 degree in the calculation of the modal mixity angle, and a scan of values every 15 degree to calculate the initial dichotomy interval.

The second method is activated with the **vectorial keyword, and is based on the calculation of \(G\) in 2 orthogonal directions. The first direction, along the tangential crack direction, gives a so-called \(G_I\) value. Then a second \(G_{II}\) value is calculated in a direction orthogonal to this tangential plane. The branching angle \(\alpha\) for the new propagation direction is then defined as: \(\alpha~=~arctan\left(\frac{G_{II}}{G_I}\right)\).

Examples#

Example 1

% Conforming mode example:
***compute_G_by_gth FRONT0
 % G-theta problem definition
 **crack_front FRONT0
 **nbnodes -8
 **elem_radius 3
 **elset NEW
 **lip lip
 % SIF output options
 **compute_K
 **compute_Ki
 **excludeBE
 % crack propagation law definition
 **fatigue 0.000000e+00 2.000000e+00
 **behavior paris
  C 1.000000e-06
  m 3.000000e+00
 **Delta_N 100
 **h -2.0
 % out-of-plane propagation
 **optim_dir 1.000000e+00 1.500000e+01

Example 2

% XFEM mode example (note, no out-of plane propagation)
***compute_G_by_gth gtheta_a
   **xfem
   **theta 0.3 0.0
   **nbnodes 3
   **behavior paris
     C 335.298
     m 3.74
   **Delta_N 2000

***compute_G_by_gth

On this page

***compute_G_by_gth#

Description#

Syntax#

Examples#

`***compute_G_by_gth`#