Fracture mechanic calculation results understanding

[François EDF Energy]

Hi all,

I launched a fracture mechanic analysis with this command line :

mpirun -np 4 ./crack_propagation 
\ -my_file MyFile.h5m 
\ -my_order 2 
\ -mofem_mg_verbose 1 -mofem_mg_coarse_order 1 -mofem_mg_levels 2 
\ -my_ref 1 -my_geom_order 1 -my_ref_order 1 -material HOOKE
\ -nb_load_steps 1 -load_scale 1 2>&1 | tee log

Here is the results I obtained :

Nonlinear elastic_ solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
solution fnorm  = 6.14040397e+05
griffith force at ent    54    coords     -140.3        117          0        g1 5.076345e-02 / 5.131979e-02 ( 9.891594e-01 )        g2 7.536088e-03 ( 1.468457e-01 )        g3 0.000000e+00 ( 0.000000e+00 )
griffith force at ent    99    coords     -147.4        124         10        g1 7.164394e-01 / 1.136281e+00 ( 6.305124e-01 )        g2 8.819578e-01 ( 7.761792e-01 )        g3 0.000000e+00 ( 0.000000e+00 )
griffith force at ent  1190    coords     -145.8      122.5      9.749        g1 -1.790881e-01 / 2.221492e-01 ( -8.061613e-01 )        g2 1.273247e-01 ( 5.731493e-01 )        g3 -3.265197e-02 ( -1.469821e-01 )
griffith force at ent  1191    coords     -144.3        121       9.01        g1 2.267564e-01 / 2.683138e-01 ( 8.451167e-01 )        g2 1.433193e-01 ( 5.341482e-01 )        g3 5.776399e-03 ( 2.152852e-02 )
griffith force at ent  1192    coords       -143      119.6      7.818        g1 1.336679e-03 / 1.150214e-01 ( 1.162114e-02 )        g2 1.114684e-01 ( 9.691099e-01 )        g3 2.833611e-02 ( 2.463552e-01 )
griffith force at ent  1193    coords     -141.8      118.5      6.235        g1 4.603272e-01 / 7.969787e-01 ( 5.775903e-01 )        g2 6.304077e-01 ( 7.909970e-01 )        g3 1.608105e-01 ( 2.017751e-01 )
griffith force at ent  1194    coords       -141      117.7      4.339        g1 1.606286e-01 / 2.535961e-01 ( 6.334032e-01 )        g2 1.955429e-01 ( 7.710802e-01 )        g3 1.650468e-02 ( 6.508253e-02 )
griffith force at ent  1195    coords     -140.5      117.1      2.225        g1 -3.475661e-03 / 2.614183e-01 ( -1.329540e-02 )        g2 2.542548e-01 ( 9.725974e-01 )        g3 -6.067915e-02 ( -2.321151e-01 )
griffith force at ent  2318    coords     -146.6      123.2      9.875        g1 3.346385e-01 / 3.424017e-01 ( 9.773274e-01 )        g2 7.050118e-02 ( 2.059020e-01 )        g3 1.689762e-02 ( 4.935030e-02 )
griffith force at ent  2319    coords       -145      121.7      9.379        g1 3.198514e-01 / 4.797413e-01 ( 6.667163e-01 )        g2 3.553027e-01 ( 7.406132e-01 )        g3 4.008485e-02 ( 8.355513e-02 )
griffith force at ent  2320    coords     -143.6      120.3      8.414        g1 5.545587e-02 / 6.617463e-02 ( 8.380232e-01 )        g2 3.496400e-02 ( 5.283596e-01 )        g3 9.013655e-03 ( 1.362101e-01 )
griffith force at ent  2321    coords     -142.4      119.1      7.027        g1 -2.803904e-01 / 4.487019e-01 ( -6.248923e-01 )        g2 3.438054e-01 ( 7.662223e-01 )        g3 -6.717491e-02 ( -1.497094e-01 )
griffith force at ent  2322    coords     -141.4      118.1      5.287        g1 -2.536785e-01 / 5.678011e-01 ( -4.467735e-01 )        g2 4.931517e-01 ( 8.685289e-01 )        g3 -1.218472e-01 ( -2.145949e-01 )
griffith force at ent  2323    coords     -140.7      117.4      3.282        g1 -8.060132e-02 / 1.219186e-01 ( -6.611077e-01 )        g2 9.050557e-02 ( 7.423443e-01 )        g3 1.327824e-02 ( 1.089107e-01 )
griffith force at ent  2324    coords     -140.4        117      1.113        g1 -2.816858e-02 / 2.024589e-01 ( -1.391324e-01 )        g2 1.951563e-01 ( 9.639307e-01 )        g3 -4.593622e-02 ( -2.268916e-01 )
max g1 = 7.1644e-01
Elastic energy 6.5673e+04
Done rank = 0
Done rank = 1
Done rank = 2
Done rank = 3

Why there are two values for g1 ?
And also why there is a value in brackets fro g1, g2 and g3 ?

Regards,
François

[Lukasz Kaczmraczyk]

Hello,

The first value is the real release of energy for mode I. Second is the magnitude of the material force, let call it J, from J-integral. If those two are close to each other, it means that you have pure mode I at the node. If values are different, it is mix mode. Values of in brackets are relative values to J, i.e. (g1/J), (g2/J) and (g3/J). You have to have a case when you have mix-mode. So you have to have K-II or K-III not equal to zero. When you will run crack propagation, it will mean that crack front will rotate itself to be in pure mode I. And g2 will be very close to zero.

However something is wrong with those results, g1 should be positive. It could be discretization error. I will suggest increasing -my_ref_order and see what will happen. 

What problem is it?

Kind regards,

Lukasz

[François EDF Energy]

Hi Lukasz,

Thank you for this very clear explanation.
It is the test case 1.1, the Brick slice with a slot. I'm working on this one at the moment because it's the most simple.

What is the difference between -my_order (which I believe is the order of the elements) and -my_geom_order ?
And between -my_ref and -my_ref_order ? (I believe that if -my_ref_order = x, it add x order to -my_order value)

Can you confirm that to define the crack in my .config file, I only have to declare the crack surface ? May be something else, which I don't yet known, has to be declared.

[block_4]        # Applying crack on Crack_Init_Surface
id=200          # ID must be 200
add=SIDESETEnter

regards,
Francois

[Lukasz Kaczmraczyk]

Hello,

Correct you need only following part to define crack surface, however block has to have unique number.

[block_4]

# Crack surface

id=200   # ID must be 200

add=SIDESET

-my_order set approx. order for every element, it can be costly, usually you can have 2 (sometimes 1). This is global p-adaptivity.

-my_ref_order set approx. order in vicinity of crack front, for example if you have -my_order 2 -my_ref_order 1, approximation order at crack front is 2+1. This is very efficient way to increase accuracy of the solution. Try values 2, 3, and more if needed. This is local p-adaptivity.
-my_geom_order is geometry order, so if you have 10-node Tet you can set it to 2. For linear tetrahedrons is 1. This only controls approximation of geometry not a quality of approximation of displacements.

-my_ref refine mesh at crack front. That way has slow convergence comparing to -my_ref_order, however for some cases if you get negative g1, you have to do both -my_ref and -my_ref_order. This is local h-adaptivity.

Look at that plot where some combinations of local/global p-adaptivity and adaptivity is tested.

http://doi.org/10.5281/zenodo.569099

Kind regards,

Lukasz

[François EDF Energy]

Hello,

I'm testing different combinations of all these parameters.
I try to find the best set of parameters for my cases.

So far this set it's the only one to give positive G1 results.

mpirun -np 8 ./crack_propagation -my_file Brick_Slice_Slot.h5m  
-my_order 1

-mofem_mg_verbose 1
-mofem_mg_coarse_order 1 -mofem_mg_levels 2

-my_ref 0
-my_geom_order 1
-my_ref_order 3

-material HOOKE -nb_load_steps 1 -load_scale 1  

2>&1 | tee log_test_7





griffith force at ent    54    coords     -140.3        117          0        g1 2.084589e+00 / 6.549634e+00 ( 3.182757e-01 )        g2 6.209041e+00 ( 9.479982e-01 )        g3 0.000000e+00 ( 0.000000e+00 )
griffith force at ent    99    coords     -147.4        124         10        g1 3.055327e+01 / 3.288489e+01 ( 9.290977e-01 )        g2 1.216196e+01 ( 3.698343e-01 )        g3 0.000000e+00 ( 0.000000e+00 )
griffith force at ent  1190    coords     -145.8      122.5      9.749        g1 1.468947e+01 / 1.521805e+01 ( 9.652661e-01 )        g2 3.858733e+00 ( 2.535628e-01 )        g3 9.585190e-01 ( 6.298565e-02 )
griffith force at ent  1191    coords     -144.3        121       9.01        g1 1.008093e+01 / 1.225339e+01 ( 8.227053e-01 )        g2 6.760906e+00 ( 5.517582e-01 )        g3 1.676470e+00 ( 1.368168e-01 )
griffith force at ent  1192    coords       -143      119.6      7.818        g1 3.327441e-01 / 4.798489e+00 ( 6.934353e-02 )        g2 4.746353e+00 ( 9.891350e-01 )        g3 -6.220185e-01 ( -1.296280e-01 )
griffith force at ent  1193    coords     -141.8      118.5      6.235        g1 1.652412e+00 / 8.584084e+00 ( 1.924972e-01 )        g2 8.172831e+00 ( 9.520912e-01 )        g3 2.039820e+00 ( 2.376281e-01 )
griffith force at ent  1194    coords       -141      117.7      4.339        g1 1.180481e+00 / 6.447072e+00 ( 1.831034e-01 )        g2 6.286556e+00 ( 9.751025e-01 )        g3 -8.064833e-01 ( -1.250930e-01 )
griffith force at ent  1195    coords     -140.5      117.1      2.225        g1 8.030654e+00 / 1.789031e+01 ( 4.488828e-01 )        g2 1.552106e+01 ( 8.675678e-01 )        g3 3.829962e+00 ( 2.140802e-01 )
max g1 = 3.0553e+01
Elastic energy 1.1231e+05

Done rank = 0
Done rank = 1
Done rank = 2
Done rank = 3

Done rank = 4
Done rank = 5
Done rank = 6
Done rank = 7

So far I only run G-calculation analyses. What is the command line to run a crack propagation analysis on the last version ?

Regards,
Francois



 

[Lukasz Kaczmraczyk]

How you introduced initial crack? 

For crack propagation make mesh without crack and separately surface for initial notch. MoFEM cut the mesh and place cack automatically.  Pleas give me some time, I have to make some proper description how to do that.

BTW this is good test for crack propagation. 

Kind regards,

Lukasz

[François EDF Energy]

Here is the procedure I use :

The crack initialisation is a surface that is a part of my mesh.

I define this crack initialisation surface as a group when I create the .med mesh.

Then, when I write the .config file, I associate the number of the block that define the crack with the group's number of the crack initialisation surface provide with the command line

 ../tools/read_med -med_file mesh.med

Once this is done I create the .h5m file from my .med mesh and my .config file with

../tools/read_med -mesh_file mesh.med -meshsets_config MyFile.config

And finally launch the calculation.

It's the wrong way to do it ? Should I create the initial crack separately from my mesh ?

King regards,
Francois

[Lukasz Kaczmraczyk]

That is the correct procedure. However now when you do model, you cat geometry to set crack surface and then mesh. If you do model for crack propagation, you mesh cut surface independently from the model. This surface will grow and change while crack is propagating and MoFEM will do mesh cutting on the fly. 

Making a model for crack propagation as you can see is simpler, more like for xfem or z-crack, that you do not have to mesh model with the crack. 

This is nowhere described correctly; we will have to do proper documentation. Could you try to make model like that and send me file. I will look at it when is correct you will be able to kick-start calculations.

See example. 

1) make a model of the brick

2) mesh brick

3) create initial cutting surface, for example circle

4) mesh surface of circle
5) save surface of the circle in sideset, for example sideset 400

See example that model and surface is mesh independently. And model is cut on the fly.

[François EDF Energy]

Hi Lukasz,

For the step 5, how to save the circle in the sideset ?
I'm supposed that it has to be done with the .config by writing:

[block_4]

# Crack surface

id=400  

add=SIDESET

The thing is that my .config file is associate with the block number of my model without a crack (Brick_Slice_Slot_No_Crack.med enclosed) and not with my crack meshed (Crack.med)
 

[François EDF Energy]

Hi,

I've run a crack propagation analysis (brick_slice_slot).
After 11 load steps a G value equal to -inf appears.

griffith force at ent  1034    coords    -0.1474      0.124    0.02559        g1 -inf / 1.549025e-04 ( -inf )        g2 inf ( inf )        g3 0.000000e+00 ( 0.000000e+00 )     relative error (gc-g)/gc inf / -1.064463e-01
griffith force at ent  6039    coords    -0.1454     0.1211    0.01822        g1 -inf / 1.422999e-04 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / -1.642807e-02
griffith force at ent  6370    coords    -0.1428      0.119   0.006503        g1 -inf / 1.447768e-04 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / -3.412014e-02
griffith force at ent  6980    coords    -0.1443     0.1202    0.01388        g1 -inf / 1.435605e-04 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / -2.543184e-02
griffith force at ent  7030    coords    -0.1463     0.1224    0.02157        g1 -inf / 1.518306e-04 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / -8.450394e-02
griffith force at ent  7031    coords    -0.1435     0.1193   0.009794        g1 -inf / 1.420280e-04 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / -1.448580e-02
griffith force at ent  7032    coords    -0.1408     0.1172   0.002227        g1 -inf / 9.377892e-05 ( -inf )        g2 inf ( inf )        g3 -inf ( -inf )     relative error (gc-g)/gc inf / 3.301506e-01
griffith force at ent  7034    coords    -0.1403      0.117          0        g1 -inf / 7.079844e-05 ( -inf )        g2 inf ( inf )        g3 0.000000e+00 ( 0.000000e+00 )     relative error (gc-g)/gc inf / 4.942968e-01
max g1 = 0.0000e+00 max J = 1.5490e-04

I decided to restart the calculation from load step 9 without modifying anything, and everything worked fine.

But I'm wondering why such a result comes out, may be some of you have an explanation about it.

The log file is enclosed if you need to have a look at it.

Regards,
Francois

[Lukasz Kaczmraczyk]

Hello,

Thanks for this, this type of error is very new to me, could you send me on my email mesh file and parameters file. I will try to reproduce this. 

In some sense, this is kind of the errors which we should fix over the time, while you run those tests.

Kind regards,

Lukasz

[François EDF Energy]

It's done.

[François EDF Energy]

Hi,

The error appears again,

In the log file have a look at the lines 28776 and 171547.

Regards,
Francois

[Lukasz Kaczmraczyk]

Hello,

What you using to compress those files?

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

Can you send parameter file as well.

[François EDF Energy]

I'm using this command

tar -czvf

[Lukasz Kaczmraczyk]

Thanks,

I can restart your analysis from step 77 and is looking ok. It is a bug in the code which you reported before. I will try to look at it and will make an update. 

I am correct that restart work for you, and there is no problem after restart, at least for some time?

Kind regards,

Lukasz

[François EDF Energy]

yes you're right, when I restarted everything worked fine then.

For your information this error occurred at least one time, sometimes two for each crack propagation analysis so far (Brick_Slice_Slot and Full_Brick_Slot). 

Kind regards,
Francois

[Lukasz Kaczmraczyk]

This is an annoying problem, related to updating variables on partitioned meshes.