GriffithForceElement parameters

[D-Nguyen]

Hello,

I really want to understand your code, not only the theoretical part but also your implementation. I only found a limited theoretical aspest from the papers as well as code instructions but  sometimes it is still difficult to understand for me.

In the GriffithForceElement, BlockData has 4 component, but the roles of BlockData::E and of BlockData:: r are not clear to me. It seems that BlockData::E is not the Young modulus as its instruction. By default, BlockData::E =BlockData::gc^3 if the value of BlockData::E  is not set by -griffith_E. 

It seems that BlockData::E  afffects heavily the convergence, a larger value increase quickly the SNES function norm and the numer of iteration of propagation solver is too big. If it is large enough,  the diverge is observed.

Moreover, a larger value of BlockData:: r (>1) seems to improve a bit the convergence.  

BlockData:: r  is used in the approximation of the max-operator but BlockData::E is not clear to me.

Can you explain a bit about these two parameters? and how to correctly choose them?

Best regards,

Van Dung

[Lukasz Kaczmraczyk]

As you know Griffith criterion is not smooth, bases on it crack can or not propagating. That lead to unilateral constrains, very similar which you have in contact mechanics (bodies are in contact or not). That why we have used ideas from contact mechanics to enforce constrains at the crack front, on particular Popp papers (https://doi.org/10.1002/nme.2614) with some additional trick.

Please see attached PDF with show how it works. E acts alike g_c on the slides. This parameter has no impact on results but may influence the convergence rate. If you make it wrong you will have problems with convergence. In future, we have to add a way to estimate such that the user does not have to do it.

"r" parameter is added to regularise unilateral constrains, and improve convergence. However, in practice, we find out that we do not have to use it and always is set always to one. If you set it big it will improve convergence but will deteriorate accuracy.

We have not yet published this, a bit of work with the formulation is needed, but from the slides, you will get the idea. BTW this is very cool stuff, that we can propagate crack as it is driven by physics, very unique to MoFEM.

 

[Lukasz Kaczmraczyk]

If you look at this log file:

 1 SNES Function norm 2.441726967660e-05 [

LAMBDA_ARC_LENGTH 1.767176697813e-06,

SPATIAL_POSITION 2.300004794645e-05,

MESH_NODE_POSITIONS 5.945526210446e-06,

LAMBDA_SURFACE102 5.246020546110e-28,

LAMBDA_SURFACE200 1.050104092973e-13,

LAMBDA_CRACKFRONT_AREA 5.359912081408e-06,

LAMBDA_CRACKFRONT_AREA_TANGENT 8.093663499978e-12,

LAMBDA_EDGE 0.000000000000e+00,

LAMBDA_BOTH_SIDES 4.767699904181e-22]

delta vec nrm  = 3.82725726e-06

Area active constrains on entity 8521 diag = +1.00000e+00 delta = +9.22786e-07 lambda = -1.29294e-06

Area active constrains on entity 8510 diag = +1.00000e+00 delta = +4.45693e-07 lambda = -6.25759e-07

Area active constrains on entity 8498 diag = +1.00000e+00 delta = +5.89808e-07 lambda = -8.20307e-07

Area active constrains on entity 8427 diag = +1.00000e+00 delta = +6.26085e-07 lambda = -8.76319e-07

Area active constrains on entity 8425 diag = +1.00000e+00 delta = +5.43737e-07 lambda = -7.59464e-07

Area active constrains on entity 8423 diag = +1.00000e+00 delta = +3.44044e-07 lambda = -4.83134e-07

Area active constrains on entity 8411 diag = +1.00000e+00 delta = +5.21534e-07 lambda = -7.20527e-07

Area active constrains on entity 8408 diag = +1.00000e+00 delta = +4.05463e-07 lambda = -5.67395e-07

Area active constrains on entity 8351 diag = +0.00000e+00 delta = -1.39411e-06 lambda = +1.93273e-06

Area active constrains on entity 7397 diag = +1.00000e+00 delta = +4.45790e-07 lambda = -6.37910e-07

Area active constrains on entity 7396 diag = +1.00000e+00 delta = +5.16986e-07 lambda = -7.31198e-07

Area active constrains on entity 7395 diag = +1.00000e+00 delta = +4.31120e-07 lambda = -6.03880e-07

The LAMBDA_CRACKFRONT_AREA 5.359912081408e-06 residual ac crack front Lagrange multipliers enforcing unilateral constrains, and following lines shows unilateral constraint at each node. If diag is 1, then the node is free to move, i.e. propagate the crack. If not it is an additional force in the magnitude if lambda, to prevent the reduction of crack area. All are enforced implicitly.

[D-Nguyen]

Hello,

It is clear to me.Thanks a lot for your kind explanation ;)

Best regards,

Van Dung