**process anisotropic_failure#

Description#

This process is an extension of the weibull or batdorf models to the case of anisotropic brittle materials. A probability of failure \(P_f^i\) is evaluated on an arbitrary number of characteristic planes for the material (eg. particular crystallographic planes), and the total probability of failure \(P_f\) is then computed in the following way:

(187)#\[P_f = 1 - \prod_i \left( 1 - P_f^i \right)\]

A particular plane \(i\) is characterized by its normal vector \(\underline{n}_i\) and the weibull model parameters (strength function \(\sigma_u^i(\theta)\) and modulus \(m(\theta)\)) depend on an angle \(\theta\) that defines a particular direction in this plane:

../../_images/aniso.svg

The probability of failure \(P_f^i\) associated to failure on plane \(i\) is then calculated by means of one of the following equations:

  • Weibull mode:

    (188)#\[P_f^i = 1~-~exp\left(-\int_V \frac{dV}{V_0} ~\frac{1}{\theta_{max}}~\int_0^{\theta_{max}} \left(\frac{\sigma_{e}(\theta)}{\sigma_u^i(\theta)}\right)^{m^i(\theta)}~d\theta \right)\]

    In the above equation, the equivalent stress \(\sigma_{e}(\theta)\) is chosen as the normal stress in direction \(\theta\):

    (189)#\[\sigma_{e}(\theta) ~=~ \underline{d}_\theta ~\cdot~ \overline{\sigma} ~\cdot~ \underline{d}_\theta\]

    where \(\overline{\sigma}\) is the stress tensor at the current integration point.

    The modulus \(m^i(\theta)\) and the scaling function \(\sigma_u^i(\theta)\) are arbitrary functions of \(\theta\) that should be defined by the user. Coefficients of this function are automatically added as model coefficients by the post-processor, and need to be given in the material file. For example, one can express \(\sigma_0^i(\theta)\) in term of the strength \(s0\) (\(\theta=0\)) and \(s90\) (\(\theta=\frac{\pi}{2}\)) in the following way:

    (190)#\[\sigma_u^i(\theta) ~=~s0~cos^2(\theta)~+~s90~sin^2\left(\theta\right)\]

    with material coefficients \(s0\) and \(s90\) declared in the material file.

  • Batdorf mode: This mode is a slight modification of the previous one, where Batdorf-type statistics are preferred (integration over all possible crack sizes, characterized by their critical stress \(\sigma_c\)):

    (191)#\[P_f^i = 1~-~exp\left(-\int_V \frac{dV}{V_0} ~\frac{2}{\pi}~\int_0^{\frac{\pi}{2}}{\left\{ \int_0^{\sigma_I} P(\sigma_c,\theta)~g\left(\sigma_c,\sigma_e(\theta)\right) ~d\sigma_c~\right\}}~d\theta \right)\]

    where:

    • The equivalent stress \(\sigma_e(\theta)\) is taken as the normal stress in direction \(\theta\) (same definition as in the previous mode)

    • \(P(\sigma_c,\theta)\) is the crack density with critical stress \(\sigma_c\). As in the conventional batdorf model this quantity is expressed as a power law, but depends on \(\theta\) by means of the anisotropic scale parameter \(\sigma_u^i(\theta)\):

      (192)#\[P(\sigma_c,\theta)~=~\left( \frac{\sigma_c}{\sigma_u^i(\theta)} \right)^{m_i}\]
    • \(g\left(\sigma_c,\sigma_e(\theta)\right)\) is equal to 1 when the equivalent stress is greater than the critical stress \(\sigma_c\): \(g\left(\sigma_c,\sigma_e(\theta)\right)= 1\) if  \(\sigma_e(\theta)>=\sigma_c\), else \(g\left(\sigma_c,\sigma_e(\theta)\right)=0\)

  • Full Batdorf mode:

    This mode fully generalizes the Batdorf approach to the anisotropic case, and allows some flexibility in the choice of the equivalent stress \(\sigma_e\). Integration is taken over the whole space, using 2 angles \(\theta\) and \(\phi\) as represented in the following figure, where:

    • \(\overrightarrow{d}(\theta)\) is the direction at angle \(\theta\) in the crystallographic plane as before,

    • \(\phi\) is the angle between the \(z\) axis (normal to the crystallographic plane) and a particular \(\overrightarrow{n}(\theta,\phi)\) direction in plane \((\overrightarrow{d}(\theta),z)\)

    ../../_images/aniso_theta_phi.svg

    The probability of failure is then computed by the following equation:

    (193)#\[P_f = 1~-~exp\left[-\int_V \frac{1}{V_0}~dV~ \int_0^{\sigma_I} { d\sigma_c \left\{ \frac{1}{2 \pi} \int_0^\pi{ d\theta \int_0^\pi { P(\sigma_c,\theta,\phi)~g(\sigma_e(\theta,\phi),\sigma_c)~sin(\phi)~d\phi } } \right\} } \right]\]

    where:

    • for each gauss point, integration is performed on the half space only (which correspond to a solid angle of \(2 \pi\)) because of symmetry of the equivalent stress \(\sigma_e(\theta,\phi)\)

    • \(P(\sigma_c,\theta,\phi)\) is the density of cracks having a critical stress \(\sigma_c\):

      (194)#\[P(\sigma_c,\theta,\phi) ~=~ \frac{1}{\sigma_u(\theta,\phi)} \left( \frac{\sigma_c}{\sigma_u(\theta,\phi)} \right)^{m(\theta,\phi)}\]

      where coefficients \(\sigma_u\) and \(m\) depend on the angles \((\theta,\phi)\) defining the current direction \(\overrightarrow{n}(\theta,\phi)\).

    • \(g(\sigma_e(\theta,\phi),\sigma_c)\) is equal to 1 when the equivalent stress is higher than the crack critical stress:

      \(g(\sigma_e(\theta,\phi),\sigma_c) = 1\) if \(\sigma_e(\theta,\phi)>\sigma_c\), else \(g(\sigma_e(\theta,\phi),\sigma_c) = 0\)

    • The equivalent stress \(\sigma_e(\theta,\phi)\) is allowed to depend both on the normal stress in direction \(\overrightarrow{n}(\theta,\phi)\) and the norm of the shear stress vector acting on the facet normal to \(\overrightarrow{n}\):

      (195)#\[\begin{split}\begin{aligned} \sigma_n(\theta,\phi) & = & \overrightarrow{n}(\theta,\phi)~\cdot ~\overline{\sigma}~\cdot~\overrightarrow{n}(\theta,\phi)\\ \overrightarrow{\tau}(\theta,\phi) & = & \overline{\sigma}~\cdot~\overrightarrow{n}(\theta,\phi)~-~\sigma_n(\theta,\phi)~\overrightarrow{n}(\theta,\phi)\\ \tau(\theta,\phi) & = & \|{\overrightarrow{\tau}(\theta,\phi)}\| \end{aligned}\end{split}\]

      The implementation then allows to choose between several expression for the equivalent stress \(\sigma_e\):

      • normal stress (“ns”) : \(\sigma_e(\theta,\phi)~=~\sigma_n(\theta,\phi)\)

      • critical energy (“ce”) : \(\sigma_e(\theta,\phi)~=~\sqrt{\sigma_n(\theta,\phi)^2+\tau(\theta,\phi)^2}\)

      • maximal strength (“ms”) : \(\sigma_e(\theta,\phi)~=~\frac{1}{2} \left( \sigma_n(\theta,\phi)~+ ~\sqrt{\sigma_n(\theta,\phi)^2+\tau(\theta,\phi)^2} \right)\)

Syntax#

**process anisotropic_failure \(~\,\) *var name [ *mode weibull | batdorf | [ full_batdorf ( ns | ce | ms ) ] ] *plane \(~\,\) [ tmin tmin ] \(~\,\) [ tmax tmax ] \(~\,~\,\) normal <VECTOR> \(~\,~\,\) function <FUNCTION> \(~\,\) [ m_function <FUNCTION> ] \(~\,\) [ phi_function <FUNCTION> ] \(~\,\) [ phi_m_function <FUNCTION> ] [ *plane … ] [ *orientation | *local_frame <LOCAL_FRAME> ] [ *nb_constant_step steps ] [ *nb_sigc_step sc_steps ] [ *precision eps ]

name is the name of the variable used to store the stress in the results file (one would expect sig here)

The *mode command allows to select the model used to calculate the probability of failure. Possible options, as detailed before, are: weibull (default), batdorf or full_batdorf. In the full_batdorf mode an additional keyword is needed to define the choice of the equivalent stress expression. As described before, candidates are: ns (normal stress), ce (critical energy) or ms (maximum strength).

The *plane command is used to define the material (or crystallographic) planes characterizing anisotropy of the failure properties. Note that as many *planes definitions as needed may be included. The total probability of failure is then computed as:

(196)#\[P_f = 1 - \prod_i \left( 1 - P_f^i \right)\]

where \(P_f^i\) is the probability of failure computed on a particular plane (numbered \(i\) in the previous equation). Subcommands of the *plane keyword are the following ones:

  • the optional commands tmin and tmax allow to define the bounds (tmin and tmax values in degrees) for the \(\theta\) angle integration:

    (197)#\[P_f ~=~ 1 -exp \left\{ - \int_V \left(~\frac{1}{\theta_{max} - \theta_{min}}~ \int_{\theta_{min}}^{\theta_{max}} ... ~~d\theta \right) dV \right\}\]

    Default values are tmin=0 and tmax=90.

  • a VECTOR object is expected after the normal keyword to define the normal to the current plane.

    *plane
     normal (0.0 0.0 1.0)
    
  • a FUNCTION object is expected after the function keyword to define the \(\theta\) dependence of the strength weibull parameter \(\sigma_u(\theta)\).

    *plane
     function s0*(cos(theta)^2.0)+s90*(sin(theta)^2.0);
    

    As is conventional in Zebulon input files, the function definition must be ended by a “;” character, and all libc mathematical functions are allowed (sin, cos, sinh, log, exp, ..). Note that in this particular context a “theta” character string is expected in the definition for the name of the function variable. All other parameters in the function definition are assumed to be coefficients (eg. s0 and s90 in the previous example) and need to be given in the material file.

  • an optional m_function keyword can be used to declare a dependence of the \(m\) weibull modulus in term of the direction in the current plane. Conventions for the function definition are the same as in the previous case. Without this option, \(m\) is assumed to be a constant in all directions. The name of the coefficients expected in the material file are then the following one: m in the case of a single plane, or m1, m2 … when several planes are declared.

  • In the full_batdorf mode only, additional function definitions are needed to specify the values of the \(\sigma_u\) and \(m\) parameters in a particular \((\theta,\phi)\) direction. Corresponding keywords are phi_function (for \(\sigma_0\)) and phi_m_function (for \(m\)).

    *plane
     function s0*(cos(theta)^2.0)+s90*(sin(theta)^2.0);
     phi_function stheta*(sin(phi)^2) + sc*(cos(phi)^2);
    

    Functions in the previous example will define \(\sigma_u(\theta,\phi)\) in the following way:

    (198)#\[\sigma_u(\theta) = s0~cos^2(\theta)~+~s90~sin^2(\theta)\]
    (199)#\[\sigma_u(\theta,\phi) = \sigma_u(\theta)~sin^2(\phi)~+~sc~cos^2(\phi)\]

    where stheta is a predefined name denoting the result of the \(\sigma_0\) function in a \((\theta,\phi=0)\) direction, and s0, s90, sc are model coefficients that need to be declared in the material file. As before, those coefficients are arbitrarily named by the user, and automatically added to the model definition after parsing the function expressions.

The optional command *orientation or *local_frame command allows to transform the stress tensor \(\overline{\sigma}\) from the global frame to a local material frame prior to volume integration. Such a material frame is in general more convenient to define the anisotropic failure coefficients.

The default method used to perform integration over the \(\theta\) angle is a simple constant step method with a midpoint rule. The default number of steps is 20, and this value can be modified by means of the *nb_constant_step command.

Similarly, integration on the critical stress \(\sigma_c\) in the batdorf and full_batdorf modes is done by a constant step method and the number of steps involved can be specified by means of the *nb_sigc_step command (default value is 20).

Note that when of value of 0 is used as argument of the previous two commands, a second-order runge-kutta method will be used to perform the corresponding integrations. This method is more precise and will automatically adapt the size of the steps to verify a target integration error. However, CPU times may be too long for practical applications, in particular in the batdorf or full_batdorf modes, and this option should only be used in debug mode to check the implementation against analytical solutions. Experience has shown that a constant step method with a default value of 20 steps is in general sufficient to provide reasonable accuracy. However, if runge-kutta integration is defined (by setting *nb_constant_step 0 and/or *nb_sigc_step 0) the command *precision eps may to set the precision of the runge-kutta method (default value is eps=1.e-4).

For each plane \(i\) the probability of failure \(P_{ip}^i\) associated to each integration point is stored in the results file under the name Pname_planei, where name is the variable name (specified by the *var command) and i the plane index. The result of the volume integration of the probability of failure on each plane, and the total \(P_f\) for the complete set of planes is written in the .post file.

Examples#

Example 1

    **file integ % mandatory
    **process anisotropic_failure
     *mode weibull
     *var sig
     *plane
      tmin 0.0 tmax 90.0
      normal (0.0 0.0 1.0)
      function   sm*(cos(3.0*theta)^2.0)+sa*(sin(3.0*theta)^2.0);
      m_function mm*(cos(3.0*theta)^2.0)+ma*(sin(3.0*theta)^2.0);
     *plane
      tmin 0.0 tmax 90.0
      normal (0.0 1.0 0.0)
      function   sc*(cos(theta)^2)+sm*(sin(theta)^2);
      m_function mc*(cos(theta)^2.0)+mm*(sin(theta)^2.0);
     *plane
      tmin 0.0 tmax 90.0
      normal (1.0 0.0 0.0)
      function   sa*(cos(theta)^2)+sc*(sin(theta)^2);
      m_function ma*(cos(theta)^2.0)+mc*(sin(theta)^2.0);

    % material file
***post_processing_data
  **process anisotropic_failure
    V0  2.0
    mm  4.0
    ma  3.5
    mc  3.2
    sm  880.0
    sa  1200.0
    sc  2000.0
***return

Example 2

    **file integ % mandatory
    **process anisotropic_failure
     *mode full_batdorf ms
     *var sig
     *plane  % m,a
      tmin 0.0 tmax 90.0
      normal (0.0 0.0 1.0)
      function   sm*(cos(3.0*theta)^2.0)+sa*(sin(3.0*theta)^2.0);
      m_function mm*(cos(3.0*theta)^2.0)+ma*(sin(3.0*theta)^2.0);
      phi_function stheta*(sin(phi)^2) + sc*(cos(phi)^2);
      phi_m_function mtheta*(sin(phi)^2) + mc*(cos(phi)^2);

    % material file
***post_processing_data
**process anisotropic_failure
    V0  2.0
    mm  4.0
    ma  3.5
    mc  3.2
    sm  880.0
    sa  1200.0
    sc  2000.0
***return