Crack tip(s) definition

[tuanhie...@gmail.com]

Hello,

I am trying to reproduce the code_aster test case sslp115 (https://www.code-aster.org/V2/doc/v12/en/man_v/v3/v3.02.115.pdf) in MoFEM. In this test case, code_aster defines only one crack tip in the center of disk, while MoFEM considers automatically 2 crack tips. In consequence, the G calculated by code_aster at a point on crack tip is around 2 times bigger than MoFEM. Could you confirm me this result is logic? Can I define the crack tip in MoFEM to have a coherent result as code_aster?

Best regards,

Pham

[Lukasz Kaczmraczyk]

Hello,

Do you model it as 3D problem? Note that in 2D problem and 3D problem you have to consider the thickness of your model. Note in 2D model G represents release energy through the thickness, in 3D model it is release energy at the point of crack front. So results will differ depending on the thickness which is used in Code_Aster and MoFEM.    

Can you send me to mesh file, and how you run code. I will look at it.

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

Ok, now I get what you asking for. Yes you right to compare results between 2d and 3d model you have to sum integrate G values from 3d model over thickness, in that case  sum them from tips.

[tuanhie...@gmail.com]

Hello,

I am sorry that my question was not clear. I talk about 3D simulation (the 2D disk is extruded with a given thickness). In sslp115 test case of code_aster, only the crack front in the centre of disk (the black point in the attached figure) is considers. On the other hand, 2 crack frond are present in MoFEM simulation. The J calculated by using code_aster is 2 times bigger than J calculated for each point on the crack front.

Here are the command lines:

read_med -med_file sslp115.med -meshsets_config sslp115.config
read_med -med_file $dirs_files/sslp115_RS.rmed  

 

mpirun --allow-run-as-root -np 2 ./crack_propagation -my_file sslp115.h5m  \

-snes_monitor -ksp_monitor -ksp_type fgmres -my_max_post_proc_ref_level 1 -ksp_atol 1e-12 -ksp_rtol 0 \

-my_order 2 -mofem_mg_verbose 1 -mofem_mg_coarse_order 1 -mofem_mg_levels 2 \

-mg_coarse_ksp_type preonly -mg_coarse_pc_type lu -mg_coarse_pc_factor_mat_solver_package mumps \

-pc_mg_smoothup 10 -pc_mg_smoothdown 10 -pc_mg_galerkin true -pc_mg_type multiplicative \

-pc_type lu -pc_factor_mat_solver_package mumps \

-my_ref 0 -my_geom_order 0 -my_ref_order 0 \

-my_residual_stress_block 1 \

-my_mwls_approx_file sslp115_RS.h5m -my_internal_stress_name SIEFELGA \

-mwls_dm 1.25 -mwls_number_of_base_functions 4  | tee log_sslp115

 

[Lukasz Kaczmraczyk]

Hi,

That makes sense.  Calculated release energy expresses energy release due to the increased crack area. In the 3D model it is energy at the point of crack front, where for 2D model is at crack tip.

In case that are two crack front nodes in 3D model, on both sides, and the edge is length one, the summation of Gs gives G from 2D model  — that why G from code_aster is two times bigger.

L.

[Lukasz Kaczmraczyk]

Hello,

Those results are not good. If you apply Dirichlet boundary condition on the surface, to the group you should add faces, not nodes. Otherwise, results will fix displacements only at the nodes, not on the edges and faces in between. Look at attached file. 

Moreover, in the command line, you specify parameters for the multi-grid solver, but you have to indicate that you like to use it in the eleastic solver, you can do that as follows,

 

./crack_propagation -my_file sslp115.h5m  \
-snes_monitor -ksp_monitor -ksp_type fgmres -my_max_post_proc_ref_level 1 -ksp_atol 1e-12 -ksp_rtol 0 \
-my_order 2 -mofem_mg_verbose 1 -mofem_mg_coarse_order 1 -mofem_mg_levels 2 \

-elastic_pc_type fgmres \
-elastic_pc_type mg \
-elastic_mg_coarse_ksp_type cg \
-elastic_pc_mg_smoothup 4 \
-elastic_pc_mg_smoothdown 4 \
-elastic_mg_coarse_pc_type lu \

-my_ref 0 -my_geom_order 0 -my_ref_order 0 \
-my_residual_stress_block 1 \
-my_mwls_approx_file sslp115_RS.h5m -my_internal_stress_name SIEFELGA \

-mwls_dm 1.25 -mwls_number_of_base_functions 10

Using MG solver you can solve much bigger problem. Consider use -my_ref_order 0, that improves results, but you know about that.

I put at the end 

-mwls_number_of_base_functions 10

that uses 2nd order polynomials, i.e. 10 base functions to approximate internal stresses. You can consider using 

-mwls_dm 2

if your mesh for the thermal problem is coarse, or you using lower order elements to calculate internal stresses.

Note, you do not have to calulate J, is is here

griffith force at ent     1    coords      22.98          0          1        g1 2.543889e+01 / 2.543901e+01

And J is

J =2.543901e+01

Kind regards,

Lukasz

[tuanhie...@gmail.com]

Hi,

Thanks for very kind explanation. In fact, when I tried to apply the Dirichlet conditions on surface in this test case and this one 

https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!topic/mofem-group/dciCgI0saQk

by juste changing NODESET to SIDESET (the others parameters stayed exactly the same), the calculations took a lot of time, before the maximum number of iteration was reached. The use of elastic solver did not help me much.

Regards,

Pham

[Lukasz Kaczmraczyk]

Hello,

This can not be done NODESET to SIDESET, it should be always NODESET for Dirichlet. This have to be done in Salome, you need add to group/family faces not nodes.

The rule is like that:
Fix at point: add nodes to group.

Fix in edges: add egdes, no need to add nodes

Fix in surface: add faces, i.e. triangles, no need to add nodes or edges to group.

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

BTW if you set Dirichelt BC in SIDESET you are lucky that it converge to something at all. Essentially you run problem without kinematic constrains.