Compute Quasi-birttle Linear Elastic Phase Field Fracture Stress

Description

This material implements the unified phase-field model for mechanics of damage and quasi-brittle failure from Jian-Ying Wu Wu (2017).

Crack Surface Energy

The regularized functional $var element = document.getElementById("moose-equation-d3f8871c-44e3-40fd-a19e-d5b898d91a35");katex.render("A_d", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is given as

$var element = document.getElementById("moose-equation-6fd7fc00-0f83-49f6-aba4-ecd05b533945");katex.render("A_d(d):= \\int_{\\mathcal{B}}\\gamma(d,\\nabla d) dV", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The crack surface density function $var element = document.getElementById("moose-equation-ae77d60f-64a2-4b7c-aa03-dae56b775282");katex.render("\\gamma(d,\\nabla d)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is expressed in terms of the crack phase-field $var element = document.getElementById("moose-equation-058e766d-2c9f-4506-bac6-45b3a42093e8");katex.render("d", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and its spatial gradient $var element = document.getElementById("moose-equation-c965f33c-029d-4a85-bf7a-6ffa63e23017");katex.render("\\nabla d", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ as $var element = document.getElementById("moose-equation-e563a4c1-329e-4d3a-be3b-f0878dfbecca");katex.render("\\gamma(d,\\nabla d) = \\frac{1}{c_0}(\\frac{1}{l}\\alpha(d) + l|\\nabla d|^2)~~~~~~~\\text{with}~~~~~~~c_0 = 4\\int_0^1\\sqrt{\\alpha(\\beta)}d\\beta", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where the geometric function $var element = document.getElementById("moose-equation-7cf0b9d9-2e29-40a4-a777-a6245d7f1f6e");katex.render("\\alpha(d)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ characterizes homogeneous evolution of the crack phase-field. $var element = document.getElementById("moose-equation-98dd60fe-8130-4853-9f4b-c3eef8b01c6c");katex.render("l", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is an internal length scale regularizing the sharp crack. $var element = document.getElementById("moose-equation-df183490-36d5-4844-869f-8ec12f15c207");katex.render("c_0", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is a scaling parameter such that the regularized functional $var element = document.getElementById("moose-equation-6f55f63d-6880-4370-8b82-c638db62ac74");katex.render("A_d(d)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ represents the crack surface.

The crack geometric function $var element = document.getElementById("moose-equation-f12c4c82-37d9-41aa-beb3-8cf477cd3c83");katex.render("\\alpha(d)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ generally satisfies the following properties, $var element = document.getElementById("moose-equation-48d678c5-b62c-4d1f-a721-e1cfa58bbb03");katex.render("\\alpha(0) = 0~~~~~~~~,~~~~~~~~\\alpha(1) = 1", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

In the classical phase-field models the modeling crack geometric function have been widely adopted for brittle fracture_energy $var element = document.getElementById("moose-equation-2b53f9f3-b623-4679-bf02-066aeba206ec");katex.render("\\alpha(d) = \\begin{cases} d~~~~\\text{Linear function} \\\\ d^2~~~\\text{Quadratic function}\\end{cases}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Elastic Energy

The elastic energy is defined as $var element = document.getElementById("moose-equation-89ca4000-8221-4a8a-8141-21948d9a7a61");katex.render("\\Psi(\\varepsilon,d) = \\omega(d)\\Psi_0(\\varepsilon)", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The monotonically decreasing energetic function $var element = document.getElementById("moose-equation-66c81181-3a62-4b63-a6f7-34db1a708d8d");katex.render("w(d)\\in[0,1]", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ describes degradation of the initial strain energy $var element = document.getElementById("moose-equation-ed77d7b3-93fb-4e39-9a69-f0d8aa34cfb9");katex.render("\\Psi_0(\\varepsilon)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ as the crack phase-field evolves, satisfying the following properties Miehe et al. (2015) $var element = document.getElementById("moose-equation-428060f1-3ffd-47e4-9f73-64e40c0011ff");katex.render("\\omega'(d) < 0~~~~~\\text{and}~~~~~\\omega(0)=1,~~~~~\\omega(1) = 0,~~~~~\\omega'(1)=0", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The variation of the elastic energy gives constitutive relations $var element = document.getElementById("moose-equation-27168b6d-1e25-4411-bc8f-7c0ae401e107");katex.render("\\boldsymbol{\\sigma} = \\omega(d)\\frac{\\partial{\\Psi_0}}{\\partial{\\varepsilon}},~~~~~~~~Y=-\\frac{\\partial\\Psi}{\\partial d} = -\\omega'(d)\\mathcal{Y}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where the thermodynamic force $var element = document.getElementById("moose-equation-6861ee0b-f08f-43fb-bef3-c001c2e12cdc");katex.render("Y", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ drives evolution of the crack phase-field with the reference energy $var element = document.getElementById("moose-equation-f4255076-c505-4e98-80be-e8ccf2b15fbf");katex.render("\\mathcal{Y}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ related to the strain field $var element = document.getElementById("moose-equation-1a698681-06f5-4dcb-9d4d-794d4853fa58");katex.render("\\varepsilon", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ .

Energetic Degradation Function

A genetic expression for degradation function $var element = document.getElementById("moose-equation-dbc94f4b-8cbe-4a24-953a-d0e4966e4f0a");katex.render("\\omega(d)", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is given as $var element = document.getElementById("moose-equation-8f041a08-d343-4a65-bbf6-8b45733cc1ba");katex.render("\\omega(d) = \\frac{1}{1+\\phi(d)} = \\frac{(1-d)^p}{(1-d)^p + Q(d)},~~~~~~\\phi(d) = \\frac{Q(d)}{(1-d)^p}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ for $var element = document.getElementById("moose-equation-5f66421a-e414-4fc2-8be4-014a81a5f8f0");katex.render("p>0", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and continuous function Q(d) > 0. Jian-Ying Wu considers following polynomials $var element = document.getElementById("moose-equation-af5b4449-86f1-44b8-9467-d292c1e33c01");katex.render("Q(d) = a_1d + a_1a_2d^2 + a_1a_2a_3d^3 + ...", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where the coefficients $var element = document.getElementById("moose-equation-b161217c-5569-4335-ba4f-6b8fa47d1bb0");katex.render("a_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ are calibrated from standard material properties.

The energetic function recovers some particular examples used in the literature, such as $var element = document.getElementById("moose-equation-aa10a95b-e6e3-43b6-8771-9e0703977364");katex.render("\\omega(d) = (1-d)^2", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ when $var element = document.getElementById("moose-equation-8c2a9e0d-75be-4dba-a6e3-1ed05266ae04");katex.render("p=2", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-11a7c3d9-5878-4c4a-9cb6-43ef0c214327");katex.render("Q(d) = 1 - (1-d)^p", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ .

Decomposition Approaches

The elastic energy is usually decomposed additively to distinguish between tensile and compressive contributions. Three decomposition approaches are implemented.

Strain Spectral Decomposition

The total strain energy density is defined as $var element = document.getElementById("moose-equation-aab1a087-02b8-41fb-b296-8876b6e39bc9");katex.render("\\Psi = \\omega(d) \\Psi^{+} +\\Psi^{-},", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where $var element = document.getElementById("moose-equation-61081736-3cd2-4a06-9945-fad943da96fa");katex.render("\\psi^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the strain energy due to tensile stress, $var element = document.getElementById("moose-equation-a480ec58-804b-4694-b76c-85fa9c7b6002");katex.render("\\psi^{-}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the strain energy due to compressive stress. $var element = document.getElementById("moose-equation-7a114451-9d11-44c5-a1b2-0a865c3f25b3");katex.render("\\Psi^{\\pm} = \\lambda \\left< \\varepsilon_1 + \\varepsilon_2 + \\varepsilon_3 \\right>^2_{\\pm}/2 + \\mu \\left(\\left< \\varepsilon_1 \\right>_{\\pm}^2 + \\left< \\varepsilon_1 \\right>_{\\pm}^2 + \\left< \\varepsilon_1 \\right>_{\\pm}^2 \\right),", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where $var element = document.getElementById("moose-equation-8a9daefe-2171-41f9-aed8-f7752ca41155");katex.render("\\epsilon_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the $var element = document.getElementById("moose-equation-5ebe839a-55fc-4699-8c49-a1374aaa27bf");katex.render("i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ th eigenvalue of the strain tensor and $var element = document.getElementById("moose-equation-9f7d24a1-4056-4c8a-83b9-c671cb618282");katex.render("\\left< \\right>_{\\pm}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is an operator that provides the positive or negative part.

To be thermodynamically consistent, the stress is related to the deformation energy density according to $var element = document.getElementById("moose-equation-569b1a72-c8af-4c1c-a51a-18506c9a62b3");katex.render("\\boldsymbol{\\sigma} = \\frac{\\partial \\Psi}{\\partial \\boldsymbol{\\varepsilon}}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ Thus, $var element = document.getElementById("moose-equation-031793bb-09c1-4d2a-8401-a70e344c7569");katex.render("\\boldsymbol{\\sigma}^{\\pm} = \\frac{\\partial \\Psi^{\\pm}}{\\partial \\boldsymbol{\\varepsilon}} = \\sum_{a=1}^3 \\left( \\lambda \\left< \\varepsilon_1 + \\varepsilon_2 + \\varepsilon_3 \\right>_{\\pm} + 2 \\mu \\left< \\varepsilon_a \\right>_{\\pm} \\right) \\boldsymbol{n}_a \\otimes \\boldsymbol{n}_a,", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where $var element = document.getElementById("moose-equation-78de6b62-0238-47a6-a6a8-889daa543a0a");katex.render("\\boldsymbol{n}_a", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the $var element = document.getElementById("moose-equation-897c5d5c-5077-47c8-a48f-9df39c304e18");katex.render("a", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ th eigenvector. The stress becomes $var element = document.getElementById("moose-equation-fe82cba0-3570-40f3-8f24-04abb05fedf6");katex.render("\\boldsymbol{\\sigma} = \\left[(1-c)^2(1-k) + k \\right] \\boldsymbol{\\sigma}^{+} - \\boldsymbol{\\sigma}^{-}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Strain Volumetric and Deviatoric Decomposition

The approach is based on the orthogonal decomposition of the linearized strain tensor in its spherical and deviatoric components: $var element = document.getElementById("moose-equation-1e9ef6f3-c856-4c18-a7b3-67a975baf1b8");katex.render("\\varepsilon = \\varepsilon_S + \\varepsilon_D,~~~~\\varepsilon_D = \\frac{1}{n}tr(\\varepsilon)I,~~~~\\varepsilon_D = \\varepsilon - \\frac{1}{n}tr(\\varepsilon)I", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ werhe $var element = document.getElementById("moose-equation-6cc2fd10-b72e-406d-884e-6114796d6ad3");katex.render("I", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ deontes the n-dimensional identity tensor.

$var element = document.getElementById("moose-equation-ca18e463-8eda-4953-a5b1-3afbc5b4e7ca");katex.render("\\psi^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-cd569167-3d6b-41fe-819f-a77ea97d785e");katex.render("\\psi^{-}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is defined as $var element = document.getElementById("moose-equation-fb9cf649-8184-4736-afd4-e3d622c695e6");katex.render("\\psi^{+} = \\frac{1}{2}\\kappa\\left<tr(\\varepsilon)^2\\right>_{+} + \\mu\\varepsilon_{D}\\cdot\\varepsilon_{D}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ $var element = document.getElementById("moose-equation-6bd6f4c8-d084-4201-8a4e-b7505480cdd6");katex.render("\\psi^{-} = \\frac{1}{2}\\kappa\\left<tr(\\varepsilon)^2\\right>_{-}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The stress is defined as $var element = document.getElementById("moose-equation-1ced15ac-1017-4c02-afee-2fcc4bdcc5e6");katex.render("\\boldsymbol{\\sigma}^- = \\kappa \\left<tr(\\varepsilon) \\right>", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-a0e97b8b-0c64-4021-b17c-66f38bb3a991");katex.render("\\boldsymbol{\\sigma}^+ = \\boldsymbol{\\mathcal{C}}\\varepsilon -\\boldsymbol{\\sigma}^-", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Stress Spectral Decomposition

$var element = document.getElementById("moose-equation-ef2556c6-ea7c-438d-8378-2b825d86c3ef");katex.render("\\psi^{+}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-9323b897-58bf-4dfb-a271-b823298a4bcd");katex.render("\\psi^{-}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is defined as $var element = document.getElementById("moose-equation-08982e05-56b6-47f9-babf-8aa2ef96c624");katex.render("\\psi^{+} = \\frac{1}{2} \\boldsymbol{\\sigma}^{+} : \\boldsymbol{\\varepsilon}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ $var element = document.getElementById("moose-equation-83f62080-3005-408f-9fb5-5fe5d1df784d");katex.render("\\psi^{-} = \\frac{1}{2} \\boldsymbol{\\sigma}^{-} : \\boldsymbol{\\varepsilon}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The compressive and tensile parts of the stress are computed from postive and negative projection tensors (computed from the spectral decomposition) according to $var element = document.getElementById("moose-equation-9c11b839-9bb4-4de4-96d3-746c9388f851");katex.render("\t\\boldsymbol{\\sigma}^+ = \\mathbf{P}^+ \\boldsymbol{\\sigma}_0", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ $var element = document.getElementById("moose-equation-33faba89-4c66-411f-9068-8db6ed92edd6");katex.render("\t\\boldsymbol{\\sigma}^- = \\mathbf{P}^- \\boldsymbol{\\sigma}_0,", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

To be thermodynamically consistent, the stress is related to the deformation energy density according to $var element = document.getElementById("moose-equation-789d56b6-180c-4109-9094-ba4cada402df");katex.render(" \\boldsymbol{\\sigma} = \\frac{\\partial \\psi_e}{\\partial \\boldsymbol{\\varepsilon}} = \\omega{d}\\frac{\\partial \\psi^+}{\\partial \\boldsymbol{\\varepsilon}} + \\frac{\\partial \\psi^-}{\\partial \\boldsymbol{\\varepsilon}}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ Since $var element = document.getElementById("moose-equation-94680e17-fdff-4b44-9fe1-bd9a4c60517d");katex.render("\t\\frac{\\partial \\psi^+}{\\partial \\boldsymbol{\\varepsilon}} = \\boldsymbol{\\sigma}^+", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ $var element = document.getElementById("moose-equation-26fc5899-664a-4768-9341-3979b0acd1f2");katex.render("\t\\frac{\\partial \\psi^-}{\\partial \\boldsymbol{\\varepsilon}} = \\boldsymbol{\\sigma}^-,", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ then, $var element = document.getElementById("moose-equation-453feeb1-7a4f-4d00-be17-134587c9edf0");katex.render("\t\\boldsymbol{\\sigma} = \\omega(d)\\boldsymbol{\\sigma}^+ + \\boldsymbol{\\sigma}^-", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The Jacobian matrix for the stress is $var element = document.getElementById("moose-equation-3153da3d-b83e-437d-b58a-16992cfd0ed1");katex.render(" \\boldsymbol{\\mathcal{J}} = \\frac{\\partial \\boldsymbol{\\sigma}}{\\partial \\boldsymbol{\\varepsilon}} = \\left((\\omega(d) \\boldsymbol{\\mathcal{P}}^+ + \\boldsymbol{\\mathcal{P}}^- \\right) \\boldsymbol{\\mathcal{C}},", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where $var element = document.getElementById("moose-equation-e67916d4-4cde-4e7d-85a4-e37cb6bf020b");katex.render("\\boldsymbol{\\mathcal{C}}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the elasticity tensor.

Note that stress spectral decomposition approach can be used for anisotropic elasticity tensor.

Evolution Equation (Allen-Cahn)

To avoid crack healing, a history variable $var element = document.getElementById("moose-equation-2b3e4366-96dc-4f39-8dda-0e1fa61351bc");katex.render("H", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is defined that is the maximum energy density over the time interval $var element = document.getElementById("moose-equation-944dbd99-d1ed-43c3-ae82-5c3d0744530d");katex.render("t=[0,t_0]", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ , where $var element = document.getElementById("moose-equation-75d08c1c-43a4-4733-a5c7-0933bbf5f17f");katex.render("t_0", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the current time step, i.e. $var element = document.getElementById("moose-equation-b9261926-5ca8-4bf1-bdb0-4968ef6286bc");katex.render("H = \\max_t (\\Psi^{+})", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Now, the total free energy is redefined as: $var element = document.getElementById("moose-equation-81025b9d-b0ee-4b29-9603-deedb9350903");katex.render("\\begin{aligned} F =& \\omega(d) H +\\Psi^{-} + \\frac{1}{c_0}(\\frac{g_c}{l}\\alpha(d) + l|\\nabla d|^2) \\\\ =& f_{elastic} + f_{fracture} + \\frac{g_c l}{c_0} {|{\\nabla d}|}^2 \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ with $var element = document.getElementById("moose-equation-58c1276e-0e4b-4f81-bce3-58c1483ff6f1");katex.render("f_{loc} = f_{elastic} + f_{fracture}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-81fc3f3a-a5ab-487b-a849-0d971a97b3e5");katex.render("f_{fracture} = \\frac{1}{c_0}\\frac{g_c}{l}\\alpha(d),~~~~~f_{elastic} = \\omega(d)H +\\Psi^{-}.", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

Its derivatives are $var element = document.getElementById("moose-equation-39ea5895-56a6-4f3e-900f-bcb0bed93b25");katex.render("\\begin{aligned} \\frac{\\partial f_{elastic}}{\\partial d} =& \\omega'(d)H\\\\ \\frac{\\partial^2 f_{elastic}}{\\partial d^2} =& \\omega''(d)H \\\\ \\frac{\\partial f_{fracture}}{\\partial d} =& \\frac{1}{c_0}\\frac{g_c}{l} \\alpha'(d) \\\\ \\frac{\\partial^2 f_{fracture}}{\\partial d^2} =& \\frac{1}{c_0}\\frac{g_c}{l} \\alpha''(d). \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

To further avoid crack phase-field going to negative, $var element = document.getElementById("moose-equation-439221f5-eb44-4c84-93b9-c508f603aeb7");katex.render("H", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ should overcome a barrier energy. The barrier energy $var element = document.getElementById("moose-equation-1d4d6f2a-39a8-470a-a251-8da7123780b0");katex.render("f_{barrier}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is determined by $var element = document.getElementById("moose-equation-5e577322-bbe2-4516-972a-8f7b03c41fee");katex.render("f_{barrier} = -\\frac{f_'{fracture}}{\\omega'(d)}~/text{at}~d=0", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and the $var element = document.getElementById("moose-equation-7a0c3fca-35ca-486b-84e9-420d1a43c73c");katex.render("H", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is modified as $var element = document.getElementById("moose-equation-4e3c162e-606a-4e21-b91b-caa22112268d");katex.render("H = \\max_t (\\Psi^{+}, f_{barrier})", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$

The evolution equation for the damage parameter follows the Allen-Cahn equation $var element = document.getElementById("moose-equation-b7c3bb5e-6777-4ff7-a141-d7bab38751ad");katex.render("\\dot{d} = -L \\frac{\\delta F}{\\delta d} = -L \\left( \\frac{\\partial f_{loc}}{\\partial d} - \\nabla \\cdot \\kappa \\nabla d \\right),", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ where $var element = document.getElementById("moose-equation-05f4ec44-44aa-40f4-a57a-04a997c3934d");katex.render("L = (g_c \\tilde\\eta)^{-1}", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ and $var element = document.getElementById("moose-equation-986c8c7f-e341-4d78-93f5-6725303adeb1");katex.render("\\kappa = 2g_cl/c_0", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ .The $var element = document.getElementById("moose-equation-70fb7fa6-a228-4fc7-b20f-6279986121ca");katex.render("\\tilde\\eta = \\eta/g_c", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is scaled by the $var element = document.getElementById("moose-equation-880fd2c4-7613-4568-ba7b-30da6280778e");katex.render("g_c", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ which is consistent with the definition given by Miehe at.al Miehe et al. (2015).

This equation follows the standard Allen-Cahn and thus can be implemented in MOOSE using the standard Allen-Cahn kernels, TimeDerivative, AllenCahn, and ACInterface. There is now an action that automatically generates these kernels: NonconservedAction. See the PhaseField module documentation for more information.

Example Input File Syntax

[./damage_stress]
  type = ComputeLinearElasticPFFractureStress
  c = c
  E_name = 'elastic_energy'
  D_name = 'degradation'
  F_name = 'local_fracture_energy'
  decomposition_type = strain_spectral
[../]

/opt/civet/build_0/moose/modules/combined/test/tests/phase_field_fracture/crack2d_iso.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
  [./TensorMechanics]
    [./Master]
      [./mech]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'stress_yy'
        save_in = 'resid_x resid_y'
      [../]
    [../]
  [../]
[]

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.04 1e-4'
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = strain_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = 2
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = 2
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'

  nl_rel_tol = 1e-8
  l_max_its = 10
  nl_max_its = 10

  dt = 1e-4
  dtmin = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

Input Parameters

cName of damage variable
C++ Type:std::vector
Options:
Description:Name of damage variable

Required Parameters

D_namedegradationName of material property for energetic degradation function.
Default:degradation
C++ Type:MaterialPropertyName
Options:
Description:Name of material property for energetic degradation function.
E_nameelastic_energyName of material property for elastic energy
Default:elastic_energy
C++ Type:MaterialPropertyName
Options:
Description:Name of material property for elastic energy
F_namelocal_fracture_energyName of material property for local fracture energy function.
Default:local_fracture_energy
C++ Type:MaterialPropertyName
Options:
Description:Name of material property for local fracture energy function.
barrier_energyName of material property for fracture energy barrier.
C++ Type:MaterialPropertyName
Options:
Description:Name of material property for fracture energy barrier.
base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Options:
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
blockThe list of block ids (SubdomainID) that this object will be applied
C++ Type:std::vector
Options:
Description:The list of block ids (SubdomainID) that this object will be applied
boundaryThe list of boundary IDs from the mesh where this boundary condition applies
C++ Type:std::vector
Options:
Description:The list of boundary IDs from the mesh where this boundary condition applies
computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Options:
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeSubdomainProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE ELEMENT SUBDOMAIN
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeSubdomainProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
decomposition_typenoneDecomposition approaches. Choices are: strain_spectral strain_vol_dev stress_spectral none
Default:none
C++ Type:MooseEnum
Options:strain_spectral strain_vol_dev stress_spectral none
Description:Decomposition approaches. Choices are: strain_spectral strain_vol_dev stress_spectral none
use_current_history_variableFalseUse the current value of the history variable.
Default:False
C++ Type:bool
Options:
Description:Use the current value of the history variable.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Options:
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Options:
Description:The seed for the master random number generator

Advanced Parameters

output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector
Options:
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
outputsnone Vector of output names were you would like to restrict the output of variables(s) associated with this object
Default:none
C++ Type:std::vector
Options:
Description:Vector of output names were you would like to restrict the output of variables(s) associated with this object

Outputs Parameters

Input Files

modules/combined/test/tests/phase_field_fracture/crack2d_iso_wo_time.i
modules/combined/test/tests/phase_field_fracture/crack2d_aniso_hist_false.i
modules/combined/test/tests/phase_field_fracture/void2d_iso.i
modules/combined/test/tests/phase_field_fracture/crack2d_vol_dev.i
modules/combined/test/tests/phase_field_fracture/crack2d_linear_fracture_energy.i
modules/combined/test/tests/phase_field_fracture/crack2d_iso.i
modules/combined/test/tests/phase_field_fracture/crack2d_aniso.i

modules/combined/test/tests/phase_field_fracture/crack2d_iso_wo_time.i

#This input does not add time derivative kernel for phase field equation
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./mech]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'stress_yy'
        save_in = 'resid_x resid_y'
      [../]
    [../]
  [../]
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
  [../]
[]

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
  [./ACBulk]
    type = AllenCahn
    variable = c
    f_name = F
  [../]

  [./ACInterface]
    type = ACInterface
    variable = c
    kappa_name = kappa_op
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.04 1e-4'
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
  [../]
  [./elastic]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = strain_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = 2
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = 2
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'

  nl_rel_tol = 1e-8
  l_max_its = 10
  nl_max_its = 10

  dt = 1e-4
  dtmin = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/crack2d_aniso_hist_false.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 40
    ny = 20
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./All]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'strain_yy stress_yy'
        planar_formulation = PLANE_STRAIN
      [../]
    [../]
  [../]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.05 1e-6'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '127.0 70.8 70.8 127.0 70.8 127.0 73.55 73.55 73.55'
    fill_method = symmetric9
    euler_angle_1 = 30
    euler_angle_2 = 0
    euler_angle_3 = 0
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = stress_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '1.0e-6'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./av_stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./av_strain_yy]
    type = SideAverageValue
    variable = disp_y
    boundary = top
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_factor_mat_solving_package'
  petsc_options_value = 'lu superlu_dist'

  nl_rel_tol = 1e-8
  l_tol = 1e-4
  l_max_its = 100
  nl_max_its = 10

  dt = 2e-6
  num_steps = 5
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/void2d_iso.i

[Mesh]
  type = FileMesh
  file = void2d_mesh.xda
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./All]
        add_variables = true
        strain = SMALL
        additional_generate_output = stress_yy
      [../]
    [../]
  [../]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        mobility = L
        kappa = kappa_op
      [../]
    [../]
  [../]
[]

[Functions]
  [./tfunc]
    type = ParsedFunction
    value = t
  [../]
  [./void_prop_func]
    type = ParsedFunction
    value = 'rad:=0.2;m:=50;r:=sqrt(x^2+y^2);1-exp(-(r/rad)^m)+1e-8'
  [../]
  [./gb_prop_func]
    type = ParsedFunction
    value = 'rad:=0.2;thk:=0.05;m:=50;sgnx:=1-exp(-(x/rad)^m);v:=sgnx*exp(-(y/thk)^m);0.005*(1-v)+0.001*v'
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = tfunc
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = left
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'l visco'
    prop_values = '0.01 0.1'
  [../]
  [./pfgc]
    type = GenericFunctionMaterial
    prop_names = 'gc_prop'
    prop_values = 'gb_prop_func'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
    elasticity_tensor_prefactor = void_prop_func
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'fracture_energy'
    decomposition_type = strain_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./fracture_energy]
    type = DerivativeParsedMaterial
    f_name = fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  solve_type = PJFNK

  petsc_options_iname = '-pc_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      lu           1'

  nl_rel_tol = 1e-9
  nl_max_its = 10
  l_tol = 1e-4
  l_max_its = 40

  dt = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/crack2d_vol_dev.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
  [./TensorMechanics]
    [./Master]
      [./mech]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'stress_yy'
        save_in = 'resid_x resid_y'
      [../]
    [../]
  [../]
[]

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.04 1e-4'
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = strain_vol_dev
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = 2
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = 2
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'

  nl_rel_tol = 1e-8
  l_max_its = 10
  nl_max_its = 10

  dt = 1e-4
  dtmin = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/crack2d_linear_fracture_energy.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
  [./TensorMechanics]
    [./Master]
      [./mech]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'stress_yy'
        save_in = 'resid_x resid_y'
      [../]
    [../]
  [../]
[]

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.04 1e-4'
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l * 3 / 4'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
  [../]
  [./elastic]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'fracture_energy'
    barrier_energy = 'barrier'
    decomposition_type = strain_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./fracture_energy]
    type = DerivativeParsedMaterial
    f_name = fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = '3 * gc_prop / (8 * l) * c'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
  [./barrier_energy]
    type = ParsedMaterial
    f_name = barrier
    material_property_names = 'gc_prop l'
    function = '3 * gc_prop / 16 / l'
  [../]
[]

[Postprocessors]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = 2
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = 2
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'

  nl_rel_tol = 1e-8
  l_max_its = 10
  nl_max_its = 20

  dt = 1e-4
  dtmin = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/crack2d_iso.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 20
    ny = 10
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
  [./TensorMechanics]
    [./Master]
      [./mech]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'stress_yy'
        save_in = 'resid_x resid_y'
      [../]
    [../]
  [../]
[]

[AuxVariables]
  [./resid_x]
  [../]
  [./resid_y]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = top
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.04 1e-4'
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '120.0 80.0'
    fill_method = symmetric_isotropic
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = strain_spectral
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '0.0'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./resid_x]
    type = NodalSum
    variable = resid_x
    boundary = 2
  [../]
  [./resid_y]
    type = NodalSum
    variable = resid_y
    boundary = 2
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm      31                  preonly       lu           1'

  nl_rel_tol = 1e-8
  l_max_its = 10
  nl_max_its = 10

  dt = 1e-4
  dtmin = 1e-4
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

modules/combined/test/tests/phase_field_fracture/crack2d_aniso.i

#This input uses PhaseField-Nonconserved Action to add phase field fracture bulk rate kernels
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 40
    ny = 20
    ymax = 0.5
  []
  [./noncrack]
    type = BoundingBoxNodeSetGenerator
    new_boundary = noncrack
    bottom_left = '0.5 0 0'
    top_right = '1 0 0'
    input = gen
  [../]
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [./TensorMechanics]
    [./Master]
      [./All]
        add_variables = true
        strain = SMALL
        additional_generate_output = 'strain_yy stress_yy'
        planar_formulation = PLANE_STRAIN
      [../]
    [../]
  [../]
  [./PhaseField]
    [./Nonconserved]
      [./c]
        free_energy = F
        kappa = kappa_op
        mobility = L
      [../]
    [../]
  [../]
[]

[Kernels]
  [./solid_x]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_x
    component = 0
    c = c
  [../]
  [./solid_y]
    type = PhaseFieldFractureMechanicsOffDiag
    variable = disp_y
    component = 1
    c = c
  [../]
  [./off_disp]
    type = AllenCahnElasticEnergyOffDiag
    variable = c
    displacements = 'disp_x disp_y'
    mob_name = L
  [../]
[]

[BCs]
  [./ydisp]
    type = FunctionDirichletBC
    variable = disp_y
    boundary = top
    function = 't'
  [../]
  [./yfix]
    type = DirichletBC
    variable = disp_y
    boundary = noncrack
    value = 0
  [../]
  [./xfix]
    type = DirichletBC
    variable = disp_x
    boundary = right
    value = 0
  [../]
[]

[Materials]
  [./pfbulkmat]
    type = GenericConstantMaterial
    prop_names = 'gc_prop l visco'
    prop_values = '1e-3 0.05 1e-6'
  [../]
  [./elasticity_tensor]
    type = ComputeElasticityTensor
    C_ijkl = '127.0 70.8 70.8 127.0 70.8 127.0 73.55 73.55 73.55'
    fill_method = symmetric9
    euler_angle_1 = 30
    euler_angle_2 = 0
    euler_angle_3 = 0
  [../]
  [./define_mobility]
    type = ParsedMaterial
    material_property_names = 'gc_prop visco'
    f_name = L
    function = '1.0/(gc_prop * visco)'
  [../]
  [./define_kappa]
    type = ParsedMaterial
    material_property_names = 'gc_prop l'
    f_name = kappa_op
    function = 'gc_prop * l'
  [../]
  [./damage_stress]
    type = ComputeLinearElasticPFFractureStress
    c = c
    E_name = 'elastic_energy'
    D_name = 'degradation'
    F_name = 'local_fracture_energy'
    decomposition_type = stress_spectral
    use_current_history_variable = true
  [../]
  [./degradation]
    type = DerivativeParsedMaterial
    f_name = degradation
    args = 'c'
    function = '(1.0-c)^2*(1.0 - eta) + eta'
    constant_names       = 'eta'
    constant_expressions = '1.0e-6'
    derivative_order = 2
  [../]
  [./local_fracture_energy]
    type = DerivativeParsedMaterial
    f_name = local_fracture_energy
    args = 'c'
    material_property_names = 'gc_prop l'
    function = 'c^2 * gc_prop / 2 / l'
    derivative_order = 2
  [../]
  [./fracture_driving_energy]
    type = DerivativeSumMaterial
    args = c
    sum_materials = 'elastic_energy local_fracture_energy'
    derivative_order = 2
    f_name = F
  [../]
[]

[Postprocessors]
  [./av_stress_yy]
    type = ElementAverageValue
    variable = stress_yy
  [../]
  [./av_strain_yy]
    type = SideAverageValue
    variable = disp_y
    boundary = top
  [../]
[]

[Preconditioning]
  [./smp]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient

  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_factor_mat_solving_package'
  petsc_options_value = 'lu superlu_dist'

  nl_rel_tol = 1e-8
  l_tol = 1e-4
  l_max_its = 100
  nl_max_its = 10

  dt = 5e-5
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

References

Christian Miehe, Lisa-Marie Schänzel, and Heike Ulmer. Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering, 294:449 – 485, 2015. URL: http://www.sciencedirect.com/science/article/pii/S0045782514004423, doi:https://doi.org/10.1016/j.cma.2014.11.016.

@article{Miehe2015,
    author = {Miehe, Christian and Sch{\"a}nzel, Lisa-Marie and Ulmer, Heike},
    title = "Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids",
    journal = "Computer Methods in Applied Mechanics and Engineering",
    volume = "294",
    pages = "449 - 485",
    year = "2015",
    issn = "0045-7825",
    doi = "https://doi.org/10.1016/j.cma.2014.11.016",
    url = "http://www.sciencedirect.com/science/article/pii/S0045782514004423"
}

Jian-Ying Wu. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. Journal of the Mechanics and Physics of Solids, 103:72 – 99, 2017. URL: http://www.sciencedirect.com/science/article/pii/S0022509616308341, doi:https://doi.org/10.1016/j.jmps.2017.03.015.

@article{JYWu2017,
    author = "Wu, Jian-Ying",
    title = "A unified phase-field theory for the mechanics of damage and quasi-brittle failure",
    journal = "Journal of the Mechanics and Physics of Solids",
    volume = "103",
    pages = "72 - 99",
    year = "2017",
    issn = "0022-5096",
    doi = "https://doi.org/10.1016/j.jmps.2017.03.015",
    url = "http://www.sciencedirect.com/science/article/pii/S0022509616308341"
}

Description
Crack Surface Energy
Elastic Energy
Evolution Equation ( Allen - Cahn )
Example Input File Syntax
Input Parameters
Input Files
References

Application Usage

Physics and Syntax

Application Development

Framework Development

MOOSEDocs

Infrastructure

Questions

Information and Tools