Softening Law Implementation

Introduction

Softening Laws are used to implement damage mechanics in softening materials.

Softening Law Definition

Implementation of softening laws starts by define a model that gives strength as a function damage parameter δ. The law if defined in the form

      [math]\displaystyle{ F(\delta,s) = \sigma_{init}f(\delta,s) }[/math]

Because the initiation stress is provided by the choosen initiation law, the softening law depends only on [math]\displaystyle{ f(\delta,s) }[/math] that is equal to 1 when [math]\displaystyle{ \delta=0 }[/math] and [math]\displaystyle{ s }[/math] is a scaling function to connecting total energy released to particle volume and crack orientation:

      [math]\displaystyle{ s = {A_c \over V_p \sigma_{init}} }[/math]

Writing Softening Law Code

The first steps are to create source files for the new softening law class, to give the law a unique name and number, and to allow the main NairnMPM code to instantiate this softening law when it is needed.

Class Source Code

The parent class for all softening laws are in NewMaterial.cpp and NewMaterial.hpp files. Normally, you will create new files that subclass this class (which can be done by starting with new files or duplicating files from an existing softening law). You can pick any unique name for the softening law source files.

Editing Core Code

Almost all coding will be done in the new softening laws class files, but for that law to be recognized as an option in NairnMPM, one source file in the core code has to be edited first. These should be the only changes needed outside the new softening laws's class files.

1. Include the new softening law's header file at the top of NairnMPM/src/Materials/MaterialBaseMPM.cpp:

    #include "Materials/MySofteningLaw.hpp"
   

   or the appropriate relative path to the new law's header file.

2. In void MaterialBase::SetSofteningLaw(char *lawName,int mode), instantiate you new law based on user input of a name or a number. Pick any name and number that does not conflict with currently implemented softening laws. The new code will look like:

    else if(strcmp(lawName,"MyFavoriteLaw")==0 || strcmp(lawName,"9")==0)
    {   sLaw = new MySofteningLaw();
    }
   

Implementing Required Methods

Once source code is created and made available to core code, the next steps are to implement all the required methods for softening laws. The required methods are:

1. const char *MySofteningLaw::GetSofteningLawName(void)
   Return a string with name for the law (to be displayed in output files). It should be unique name among all implemented softening laws.
2. double MySofteningLaw::GetFFxn(double delta,double s) const
   Return [math]\displaystyle{ f(\delta,s) }[/math].
3. double MySofteningLaw::GetDeltaMax(double s) const
   Return [math]\displaystyle{ \delta_{max} }[/math] or the value of [math]\displaystyle{ \delta }[/math] for which [math]\displaystyle{ f(\delta_{max},s)=0 }[/math]. Some softening laws (like Exponential Softening) may never reach zero, but are advised to implement a critical value where [math]\displaystyle{ f(\delta_{max},s)={\rm min} }[/math] or where the softening law decreases to some small value (specified in min).
4. double MySofteningLaw::GetFpFxn(double delta,double gScaling) const
   Return derivative of the softening law or [math]\displaystyle{ \partial f(\delta,s)/\partial\delta }[/math].
5. double MySofteningLaw::GetGToDelta(double delta,double gScaling) const
   Return normalized energy released for material point that has evolved to [math]\displaystyle{ \delta }[/math]. It is found from:
         [math]\displaystyle{ \Omega(\delta,s) = \int_0^\delta f(u,s) du - \frac{1}{2}\delta f(\delta,s) }[/math]
6. double MySofteningLaw::GetGoverGc(double delta,double gScaling) const
   Return ratio of energy released to the total toughness. It is defined by
         [math]\displaystyle{ \frac{\Omega(\delta,s)}{\Omega(\delta_{max},s)} }[/math]
   Although this result could be found from
        GetGToDelta(delta,s)/GetGToDelta(GetDeltaMax(s),s)
   that approach is less efficient (probably for any softening law). As a result, this implementation requires implementation of a separate method to do the calculation.
7. double MySofteningLaw::GetMaxSlope(double gScaling) const
   Return the maximum value for the negative slope or
         [math]\displaystyle{ \max\left(-\frac{\partial f(\delta,s)}{\partial\delta}\right) }[/math]
   This calculation is used to verify cell sizes are sufficently small for stable calculations.
8. double MySofteningLaw::GetEtaStability(void) const
   Return dimensionless stability factor, [math]\displaystyle{ \eta }[/math] . It is given by s*G_c*GetMaxSlope(s), but the scaling factor s is not provided as a parameter. None of the current laws need s to find [math]\displaystyle{ \eta }[/math]. If that ever changes for a future softening law, this method will change to include s as an argument.
9. double MySofteningLaw::GetPhiFxn(double delta,double s) const
   Return the energy dissipation function defined by
         [math]\displaystyle{ \varphi(\delta,s) = f(\delta,s) - \delta\frac{\partial f(\delta,s)}{\partial\delta} }[/math]
   Note that the only requirement for valid softening laws is that [math]\displaystyle{ \varphi(\delta,s)\ge0 }[/math] for all [math]\displaystyle{ \delta }[/math].
10. double MySofteningLaw::GetRdFxn(double delta,double s,double ei) const
    Return function that relates evolution of damage variable [math]\displaystyle{ D }[/math] to damage variable [math]\displaystyle{ \delta }[/math] defined by:
          [math]\displaystyle{ \frac{dD}{d\delta} = \frac{\varepsilon_i\varphi(\delta,s)}{\bigl(\delta + \varepsilon_i f(\delta,s)\bigr)^2} }[/math]
    where [math]\displaystyle{ \varepsilon_i }[/math] is initiation strain provided in a parameter.

Implementing Option Methods

The required methods in the previous section are sufficient to implement a new softening law. The parent class, SofteningLaw.cpp implements calculations that work for any softening law that provides these methods. Some of those calculations, however, need numerical solutions to varios problems. If any specific law can find closed-form solutions instead, implementing one or more of the following options methods can improve calculation speed when uses you custom softening law.

1. double MySofteningLaw::GetDeltaFromDamage(double d,double s,double ei,double deltaGuess)
   The damage variable [math]\displaystyle{ D }[/math] is related to the damage variable [math]\displaystyle{ \delta }[/math] by :       [math]\displaystyle{ D = \frac{\delta}{\delta + \varepsilon_i f(\delta,s)} }[/math]
   The purpose of this function is to invert this expression to find [math]\displaystyle{ D }[/math] in terms of [math]\displaystyle{ \delta }[/math] and s. The general code inverts it numerically using Newton's method. If it can be inverted analytical for your softening law, override this method and implement the analytical solution.

Compiling with make Command

Softening Law Implementation

It is easy to add new softening laws, although not clear if that change would be much different then using the softening laws already available. As explained in Ref. , the addition of a softening only requires code to return:

      [math]\displaystyle{ f(\delta,s), \quad f'(\delta,s), \quad A(\delta,s), \quad{\rm and}\quad\max\bigl(-f'(\delta,s)\bigr) }[/math]

One internal calculation is to evolve damage. That calculation can be done from the above softening law properties by solving the discrete differential equation

      [math]\displaystyle{ d\varepsilon = d\delta + \varepsilon_0 \bigl(f(\delta+d\delta,s) - f(\delta,s)\bigr) }[/math]

where [math]\displaystyle{ d\delta }[/math] is a maximum cracking strain increment associated with strain increment [math]\displaystyle{ d\varepsilon }[/math] during damage evolution. The task is to solve for [math]\displaystyle{ d\delta }[/math], which represents the evolution of damage in the current time step. The generic softening law code solves this equation using Newton's method with bracketing. The initial brackets are:

      [math]\displaystyle{ d\varepsilon \le d\delta \le d\varepsilon + \varepsilon_0f(\delta,s) }[/math]

The lower limit is found because the second term is less than or equal to zero for monotonic softening. The upper limit is if failure occurs such that [math]\displaystyle{ f(\delta+d\delta,s)=0 }[/math]. If [math]\displaystyle{ d\delta_n }[/math] is the Newton's method sequence approaching the solution, the next increment by Newton's method is derived to be:

      [math]\displaystyle{ d\delta_{n+1} = d\delta - {{d\delta-d\varepsilon\over \varepsilon_0f(\delta,s)} + {f(\delta+d\delta,s)\over f(\delta,s)}-1\over {1\over \varepsilon_0f(\delta,s)} + {f'(\delta+d\delta,s)\over f(\delta,s)}} }[/math]

When the next increment is outside the current brackets, a custom bracketing calculation is used to make the solution always stable. Although generic code can solve this differential equation for any softening law with provided properties listed above, some softening laws may have more efficient solutions. For example Linear Softening can exactly solve the differential equation without needing Newton's method and Exponential Softening can express the solution in terms of branch zero of the Lambert W function for real arguments (although evaluation of that function needs a numerical method). Any new softening laws can override generic functions for Newton's method if they have a more efficient approach.

References