more cleavage plane in ACInterfaceCleavageFracture

[wsfsz123]

Hi, all,

The kernel modification of ACInterfaceCleavageFracture is used to simulate the anisotropic fracture with transverse isotropic (one cleavage plane), then if I want to use it to simulate the cubic symmetric with two cleavage plane or more, like three cleavage plane, refer to "S. Teichtmeister, Phase field modeling of fracture in anisotropic brittle solids, 2017". What should I do?

Thanks,
Shizhe

[aeslaughter]

@dschwen @laagesen @SudiptaBiswas Do you have an answer for this, it seems to have slipped through without a response.

[jiangwen84]

Sorry for the late reply. I think one way is to introduce additional phase field damage parameters to represent fractures on different cleavage planes, and enforce the sum of all damage parameters to be 1.

[jiangwen84]

If you want to implement Teichtmeister's approach, I think we should generalize ACInterfaceCleavageFracture to include more than one preferred planes in anisotropic metro A. We do not have that capability yet, but it should be straightforward to do. Let me know if you would like to contribute, and we can provide you some guidance.

[wsfsz123]

Thank you! I tried, but A is a second-order tensor, can only represent one cleavage plane, to introduce bi-cleavage plane, it should include Laplace (second-order) of c to the fracture density function, then Allen-Cahn is not suitable anymore, does any evolution method works for it on MOOSE?

[jiangwen84]

Are you referring to the third term of equation (38) in Teichtmeister's paper?

[wsfsz123]

Thank you, I quit this method, it seems not suitable for my problem. I am now trying to modify Gc(theta) directly (see attached figure, the reference to https://doi.org/10.1016/j.jmps.2020.104253), it is similar to ACInterfaceKobayashi kernels, but I have encountered the convergence problem even at the first step, attached is my input file. Besides, may I know how to choose the parameters or methods in executioners to improve the calculation efficiency on MOOSE?

Thanks,
Shizhe

[crack_Kobayashi.txt](https://github.com/idaholab/moose/files/5693543/crack_Kobayashi.txt)

[jiangwen84]

Let us set use_current_history_variable = false first, and see if it improves the convergence. To further simplify the model, you can also set decomposition_type = none, to reduce nonlinearity.

[wsfsz123]

I tried and failed. since gc is constant before, now I need to use gc (grad(c)), do I need to modify the code in ComputeLinearElasticPFFractureStress.C, and how?

[jiangwen84]

Yes. If your Gc is a function of normal (gradient of c), you need to account for the derivative of G_c respect to c. Current implementation assumes Gc does not depend on c. To make the change, you have to make Gc as a derivative material, and add its derivative (g_d in eqn 23 in your referred paper) to the driving force term.

[wsfsz123]

The g_d is calculated in InterfaceOrientationMaterial (defined as _deps), so I think it could be used directly (ACInterfaceKobayashi is used for the gradient part of the fracture energy), my problem is: in the non-gradient term (f_elastic+Gc/2l*c^2), do I need to modify the calculation? I am thinking that Gc depends on the gradient of c, not c, so the partial gradient of c for this term shouldn't consider the gradient of Gc. But during my simulation, it has a convergence problem even when I only consider the non-gradient term (no ACinterface term).

[jiangwen84]

First, I think you should remove ACInterfaceKobayashi2. ACInterfaceKobayashi2 is the double-well potential term, and should not be considered in phase field fracture model.

Secondly, you should definitely use chain rule to include the derivative of Gc respect to c through \nabla c. See the second term on the right side of eqn 14 , and eon 16 in Shahed Rezaei's paper.

Lastly, you should make sure you use a consistent Mobility in your input file. For phase field fracture model, the M = 1/(gc*visco). I saw you defined another M that was used in ACInterfaceKobayashi1

[wsfsz123]

I have further questions about the parameter type in MOOSE.

such as the following code in the material section, where c is the order parameter, then, my question is must gc_prop be a constant? if I have defined gc_prop to be a field that depends on the location (x,y), is it still works for this expression or should I put the gc_prop in "args='? In the other words, the parameter in material_property_names should be constant or any type?

[./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
[../]

Thanks,
Shizhe

[jiangwen84]

The gc_prop can be a function material that depends on (x,y,z,t). The args is the coupled variable name, and for your case it should not include gc_prop.

[jiangwen84]

Currently, gc_prop cannot be a function of damage parameter because we do not account for the derivative of gc_prop respect to the damage variable in our material class.

If that is something you need in the future, It should be straightforward to modify the codes to make gc_prop that depends on damage variable.

[wsfsz123]

Then, do I need to modify the definition of gc_prop in ComputePFFractureStressBase.C : _gc(getMaterialProperty("gc_prop"))? Besides, I can't find any code of *.i have use the gc_prop directly except the definition in ComputePFFractureStressBase.C.

[jiangwen84]

If you use a function material to define gc_prop, you do not need to change the source code.

For regression tests, we usually define gc_prop as GenericConstantMaterial

[./pfbulkmat]
type = GenericConstantMaterial
prop_names = 'gc_prop l visco'
prop_values = '1e-3 0.04 1e-4'
[../]

[wsfsz123]

Thank you! and how to solve the negative damage order parameter problem? could I restrict the range of the order parameter to be [0,1]?

[jiangwen84]

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