J calculation - force magnitude

[tuanhie...@gmail.com]

Dear all,

I am trying to run the code_aster's sslv154 test case (https://www.code-aster.org/V2/doc/v11/fr/man_v/v3/v3.04.154.pdf). The inclined forces are applied on (NODESET) surfaces with 2 following methods:

- Either directly on surface groups (sslv154.config)

- Or on node groups based on surface groups created by using Salome (sslv154.configNode).

The force magnitude stays the same. I obtain very different J for 2 methods, and the MoFEM's results are different from code_aster reference results. Could you please give me some advice to calculate correctly the force?

I send you in the attached file the mesh file and 2 config file corresponding to 2 methods.

Thanks,

[Lukasz Kaczmraczyk]

Hello,

Problem is with the forces, in your group for nodal forces you have only nodes. As you know in mofem we heave hierarchical approximation, and DOFs are not only on nodes but can be on edges, fasces and volumes. So if you apply forces to surface, you should have in the group faces on that surface. If you apply forces on the edges, you should have in the group edges. If you have forces on all nodes only, you apply forces on the nodes, but there is no force on edges or faces. This is a case how you would apply force by a set of tiny needles. See the attached figure. To fix this, to the group with the forces, add faces.

For the problem with pressure, the problem is not restrained and you have some rigid rotation.

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

Just to let you know, to get a figure, I solved the elastic problem, and not introduced crack. This is good enough to check boundary conditions.

Just to let you know, to get a figure, I solved the elastic problem, and not introduced crack. This is good enough to check boundary conditions.

./read_med -med_file sslv154.med -meshsets_config sslv154.configNode
./elasticity -my_file out.h5m -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monito

The elasticity problem is located in users modules, in sub-directory

basic_finite_elements/elasticity/elasticity

KInd regards,

Lukasz

[tuanhie...@gmail.com]

Hello,

Thanks a lot for your answer, that helps me understand better now about the applied force. For a total applied force of F=1 kN on a surface S (n1 nodes, n2 elements), how can I calculate the value of force_magnitude in config file?

My test case has been performed with convergence, so I do not really get your remark about the rigid rotation...

Best regards,

Pham

[Lukasz Kaczmraczyk]

Hello,

To apply an overall force of 1kN on the surface, assuming that force is uniform, you have to apply force traction = 1kN/Surface Area. To apply an overall force of 1kN on the surface, assuming that force is uniform, you have to apply force traction = 1kN/Surface area. You can apply force on the edges, or nodes only, then you create group not with the faces, but edges or nodes, and then set force nodeset in the config file. 

I could be wrong; when I plotted displacements for a case with pressure, rotations looks extremely large. Usually, this is a case when a problem is not constrained correctly, however it could be in this case as it is.

Regards,

Lukasz

[tuanhie...@gmail.com]

Hello,

It means that the following input:

[block_13]

# Applying force on surface Up

id=1013

add=NODESET

force_magnitude=1000

force_fx=1

force_fz=1

gives me 1kN force in x direction and 1kN in z direction? Or the magnitude of force vector is equal to 1kN?

Regards,

Pham

[tuanhie...@gmail.com]

Hi, 

By using your method of force calculation and the data of sslv154 test case 

[block_5]

# Applying force on Left

id=1005

add=NODESET

force_magnitude=500000

force_fx=1

force_fz=1

[block_6]

# Applyng pressure on Low

id=1006

add=NODESET

force_magnitude=500000

force_fx=1

force_fz=1

[block_11]

# Applying pressure on Right

id=1011

add=NODESET

force_magnitude=500000

force_fx=-1

force_fz=-1

[block_13]

# Applying pressure on Up

id=1013

add=NODESET

force_magnitude=500000

force_fx=-1

force_fz=-1

I obtain 

J = sqrt(G1^2 + G2^2 + G3^2) 

20% bigger than the reference result. Have you an idea about this big difference?

Regards,

Pham

[Lukasz Kaczmraczyk]

Hello,

Do the results change if you use 

-my_ref_order with 1, 2,...

or higher number.

Have you checked if you apply Dirichlet boundary conditions to faces not nodes, see my response here

https://groups.google.com/d/msg/mofem-group/t6A8jhKKLjY/CFpDrNNtCQAJ

I do not see that you apply any Dirichlet boundary conditions, you always should have them. At least you can constrain rigid translation and rotation.

Kind regards,

Lukasz

[tuanhie...@gmail.com]

Hi,

The Dirichlet BC are applied on face as well. 

Regards,

Pham

[Lukasz Kaczmraczyk]

That looks OK. 

If you do not introduce crack is homogenous stress state mixture tension and sheer. I have not to check values.

Regards,

Lukasz

[tuanhie...@gmail.com]

Thanks for your answer. 

Do you mean that the issue of obtained values (different from code_aster) can relate to the mixture traction-shear stress? Almost test case with mode I crack look so good, but mixed mode cracks give me a lot of hesitation in MoFEM. :( 

[Lukasz Kaczmraczyk]

Hello,

Mofem is not designed to calculate stress intensity factors, in particular in mix modes. G1 represents crack release energy if the crack propagate in the current crack plane. You verified this when mode-I was calculated. 

Note, that we build mofem to analyse the propagation of the crack so that we do not have to calculate stress intensity factors in mix mode since our crack always rotates to maximise dissipation of energy, from the theory it can be proven that this is mode I. 

You calculating g1, g2 and g3, when a crack propagates, you can find out that components g2, g3, are negligible, and J=g1. This is in the root of the theory, and the direct consequence of maximal energy dissipation, i.e. crack release energy. 

We can change the code,  that KII, and KIII are calculated, but this will be no use to crack propagation algorithm. 

Kind regards,

Lukasz

[Lukasz Kaczmraczyk]

Note that this is not the drawback, but the biggest advantage of the theory.  Since crack always rotates to mode one, we do not need a sophisticated theory for crack initialisation mix mode, we use directly Griffith theory.  No need to heuristic algorithms choosing a crack direction, maximal hoop stress, principal stress, or any equation using K1, K2, or K3. No need for additional parameters, etc. 

[Lukasz Kaczmraczyk]

Here is case of propagation in mix mode,

http://mofem.eng.gla.ac.uk/mofem/html/frac_mech_simple_examples.html

whatever you do, crack always rotate itself to mode I, like in experiment.

[tuanhie...@gmail.com]

Thanks a lot for your answer. Evidently, the MoFEM's approach shows advantages in theoritical and numerical points of view, that's why I'd love to use it for my current and future crack studies.

I means in my previous post, that I obtained different J for mixed mode crack by using either MOFEM or Aster (as you said, the K calculation is not the objectif of MoFEM).

[Lukasz Kaczmraczyk]

Yes,

If you look at the code, you can notice that we aim higher than classical linear fracture mechanics, we dropped idea of stress intensity factors, which make sense in the case of large strains or heterogenous materials, i.e. varying material properties. 

As results, MoFEM can propagate a crack in soft material, which exhibit brittle behaviour, like gels. You can notice that all equations hold, and in fact, we have done such calculations. Current version of code can handle this. But for such case, K1, K2, or K3, can't be calculated, because those are standing next to the analytical solution, which can not be obtained. So instead of dealing with mix mode, we have implicit solution scheme rotating crack to mode I, which is unique, because we could calculate implicitly not only crack direction but as well crack area extension, exploiting equilibrium equation at crack front.

So at the end of load step crack front is mode I, that is not the case for classical approach, where crack propagation is explicit, and tiny error leads to wrong crack path. In the case of mofem, you could switch off singularity at the crack front or have very coarse mesh, and the crack path will not change too much. 

Regards,

Lukasz

[Lukasz Kaczmraczyk]

To add more.

If we find that is demand for code which can calculate K1, K2 and K3, we can consider adding this to fracture module. In fact if you like to program, you can help us do this.

Regards,

Lukasz