Softening Law Implementation
Introduction
Softening Laws are used to implement damage mechanics in softening materials. This documentation explains how to write code for user-defined softening laws.
Softening Law Definition
Implementation of softening laws starts by defining a model that gives strength as a function damage parameter δ. The law is 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 dimensionless [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 factor to connect total energy released to particle volume and crack orientation [1]:
[math]\displaystyle{ s = {A_c \over V_p \sigma_{init}} }[/math]
Note that although softening laws are traditionally monotonically decreasing functions, that is not a thermodynamic requirement. The only limitation on softening laws is that its derivative must satisfy:
[math]\displaystyle{ \frac{\partial f(\delta,s)}{\partial \delta} \le \frac{f(\delta,s)}{\delta} }[/math]
If you plot [math]\displaystyle{ f(\delta,s) }[/math] as a function of [math]\displaystyle{ \delta }[/math], this condition means no line from [math]\displaystyle{ \bigl(\delta,f(\delta,s)\bigr) }[/math] back to the origin can cross over the [math]\displaystyle{ f(\delta,s) }[/math] plot. Thermodynamically, it means that dissipated energy is nonnegative and that stiffness is monotonically decreasing. The function [math]\displaystyle{ f(\delta,s) }[/math], however, may increase to a peak before decreasing to zero at failure.
Writing Softening Law Code
The following sections list the tasks needed to write and implement a user-defined softening law in NairnMPM.
Class Source Code
The first step is to create source files for the new softening law class. The parent class for all softening laws are in SofteningLaw.cpp and SofteningLaw.hpp files. Normally, you will create new files that subclass this class (which can be done by starting with new files or by duplicating files from an existing softening law). You can pick any unique name for the softening law source files. If beneficial, a new law cal alternative subclass any other existing softening law.
Editing Core Code
Almost all coding will be done in the new softening law 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.
- 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.
- In void MaterialBase::SetSofteningLaw(char *lawName,int mode), instantiate your 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:
- 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. - double MySofteningLaw::GetFFxn(double delta,double s) const
Return [math]\displaystyle{ f(\delta,s) }[/math]. - 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, which is available to all subclasses of SofteningLaw). - double MySofteningLaw::GetFpFxn(double delta,double gScaling) const
Return derivative of the softening law or [math]\displaystyle{ \partial f(\delta,s)/\partial\delta }[/math]. - 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]
Note: that recent code changes have eliminated this method. New cohesive laws need not implement it. - double MySofteningLaw::GetGoverGc(double delta,double gScaling) const
Return ratio of energy released to the total toughness. It is defined by
[math]\displaystyle{ \frac{G}{G_c} = \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.
Note: this method is only used for cubic failure surfaces (and then only in 3D for isotropic softening materials). Cohesive laws that do not override this method are limited to other modes (but cubic failure surface is no longer recommended). - double MySofteningLaw::GetEtaStability(void) const
Return dimensionless stability factor, [math]\displaystyle{ \eta }[/math] [1]. It is given by
[math]\displaystyle{ \frac{1}{\eta} = sG_c\max\left(-\frac{\partial f(\delta,s)}{\partial\delta}\right) }[/math]
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. - double MySofteningLaw::GetPhiFxn(double delta,double s) const
Return the energy dissipation function [1] defined by
[math]\displaystyle{ \varphi(\delta,s) = f(\delta,s) - \delta\frac{\partial f(\delta,s)}{\partial\delta} }[/math]
Note that the above requirement for a valid softening law is equivalent to [math]\displaystyle{ \varphi(\delta,s)\ge0 }[/math] for all [math]\displaystyle{ \delta }[/math]. - double MySofteningLaw::GetRdFxn(double delta,double s,double ei) const
Return function that gives ratio of damage variable [math]\displaystyle{ D }[/math] evolution to damage variable [math]\displaystyle{ \delta }[/math] evolution. It is found from:
[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 Optional 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, use numerical solutions to key problems. If any specific law can find closed-form solutions instead, implementing one or more of the following optional methods can improve calculation speed for new laws. You can refer to Linear Softening or Exponential Softening coding to see examples of overriding these methods.
- 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 [1]:
[math]\displaystyle{ D = \frac{\delta}{\delta + \varepsilon_i f(\delta,s)} }[/math]
The purpose of this method is to invert this expression to find [math]\displaystyle{ \delta }[/math] in terms [math]\displaystyle{ D }[/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. - double MySofteningLaw::GetDDelta(double de,double ei,double delta,double D,double s) const
A core calculation for finite amount of damage evolution is to solve:
[math]\displaystyle{ d\varepsilon = d\delta + \varepsilon_i \bigl(f(\delta+d\delta,s) - f(\delta,s)\bigr) }[/math]
for [math]\displaystyle{ d\delta }[/math] where [math]\displaystyle{ d\varepsilon }[/math]=de is an input strain increment. The generic code solves this numerically using Newton's method with bracketing. The initial brackets are:
[math]\displaystyle{ Dd\varepsilon \le d\delta \le d\varepsilon + \varepsilon_if(\delta,s) }[/math]
The lower limit is found because [math]\displaystyle{ \varphi(\delta,s)\ge0 }[/math]. 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_n - {{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 that implements all requred methods, some softening laws may have more efficient solutions. Any new softening law can override this method to return a closed-form solution or provide a more efficient numerical solution. - double MySofteningLaw::GetRelFFxn(double delta,double x,double gScaling) const
This method returns [math]\displaystyle{ f(\delta+d\delta,s)/f(\delta,s) }[/math] that is needed in the Newton's method calculations described for GetDDelta(). If that method is overridden, this method provides no function. If it is not overridden, however, this method might provide some benefit (although to date, it has only been useful for Exponential Softening where finding the ratio is more efficient than evaluating the numerator and denominator separately). - double MySofteningLaw::GetRelFpFxn(double delta,double x,double gScaling) const
This method returns [math]\displaystyle{ [\partial f(\delta+d\delta,s)/\partial \delta]/f(\delta,s) }[/math] that is needed in the Newton's method calculations described for GetDDelta(). If that method is overridden, this method provides no function. If it is not overridden, however, this method might provide some benefit (although to date, it has only been useful for Exponential Softening where finding the ratio is more efficient than evaluating numerator and denominator separately). - double MySofteningLaw::GetDDeltaElastic(double delta,double sigmaAlpha,double s,double D,double E1md) const
This methods solves for elastic change in damage variable [math]\displaystyle{ \delta }[/math] when another parameter (such as pressure) affects the damage evolution (it implements some unpublished damage mechanics methods). The numerical task is to solve
[math]\displaystyle{ {\rm E1md*}(\delta +d\delta_e) = {\rm D*sigmaAlpha*}f(\delta+d\delta_e,s) }[/math]
for [math]\displaystyle{ d\delta_e }[/math]. The generic code solves this equation numerically. Any new softening laws can override this method to return a closed-form solution or a more-efficient numerical solution.
Including Additional Softening Law Properties
The parent class SofteningLaw.cpp mains two softening law properties for toughness, Gc, and minimum value for failure, min, used by some laws. New custom laws are likely to need additional parameters. Because softening laws are dimensionless, new parameters should be dimensionless as well. To allow users to input new parameters, the new class will need the following overrides:
MySofteningLaw::MySofteningLaw() : SofteningLaw()
{ myVariable = 0.;
}
char *MySofteningLaw::InputSofteningProperty(char *xName,int &input,double &gScaling)
{ if(strcmp(xName,"MyParameter")==0)
{ input=DOUBLE_NUM;
return (char *)&myVariable;
}
return SofteningLaw::InputSofteningProperty(xName,input,gScaling);
}
void MySofteningLaw::PrintSofteningProperties(double sigmainit)
{ SofteningLaw::PrintSofteningProperties(sigmainit);
MaterialBase::PrintProperty("MyParameter",myVariable,"");
cout << endl;
}
The first method adds a constructor (that was not needed before) and it initializes your law's parameter (which here is in a class variable myVariable). The second method overrides the method to decode input parameters. If the current one matches the name for your parameter, return a pointer of appropriate type to the class variable. If not, pass the method onto your law's super class. The third method will print information about the value used for your parameter in the simulation output file. This override calls the parent class to print its parameters and then prints this law's new variable. This scheme can be expanded to add any number of parameters with any names.
To support new parameters in input XML files, you have to add multiple definitions for that parameter to the DTD file (because of the way softening law properties are set in material definitions). The best approach is to replicate the definitions for a previous softening law parameter as follows:
- Open the NairnMPM.dtd file and search for all occurrences of "-Gc". They will show up in various forms such as "I-Gc", "II-Gc", etc. (if any).
- For each occurrence of "-Gc", add a corresponding entry for your parameter (e.g., for "I-Gc", add "I-MyParameter", where you use the name used in InputSofteningProperty() above).
Compiling with make Command
To enable compiling code with the new softening law in Linux, Unix, or Mac without using XCode, the new softening law should be added to the makefile. The best approach is to search the makefile for an existing softening law (such as LinearSoftening) and use all those entries as templates for adding your new softening law.