Princpal stress and incompressiblility

[sanshi2000]

I am doing a stress analysis at crack tip in a incompressible material. So I use modified neohookean mateial to specified the material. Unfortunately, I found the principal stress  at crack tip is negative, if the \nu is 0.4999, which indicate the material incompressible, so that means the tip is compressed which seem opposite to realily. In order to find reasons, I have also tried a linear elasitic material and solve a plane strain problem. The problem is similar, when the \nu is small as 0.49, the stress concentration is significant at tip as expected, but terrible when \nu reach 0.4999.

Anyone has idea on the reasons. The input files and principle stress ploting with different \nu are attached. The Inhkf is the finite deformation, while Inhkl is for linear plane strain problem.

Best,
Lei

[FEAP_Admin]

Why do you feel the principal stresses should be positive?  Incompressibility only means that the volumetric strains are zero:  trace(epsilon) = 0 [small strain case] or det(F) = 1 [finite strain case].

[sanshi2000]

Sorry for my confusing problem. I know the meaning of incompressibiliy as you said. As I know the Possion's ration \nu=0.5 for incompressible material, so I used a \nu=0.49~0.4999 in the neohookean model to approach the incompressible material. With the \nu is increasing and closer to 0.5, the deformation or displacement will conserve, so the big value for \nu (like 0.4999 or bigger) seem suitable for the incompressibility of material, which is named penalty method for constraint I guess.

At same time, with increasing the \nu (or correspondingly the bulk modulus K) from 0.49~0.49999, the principal stress at crack tip is always changing from positive to negative and proportional to the change of \nu after \nu is beyond some value. Should the stress be like deformation conserve when \nu is beyond some value?

In practice, I guess the crack will propagate under pressure on the crack surface. Associated with propagating, the stress at crack tip should be positive but I got negative in computing. A negative stress meaning the material is compressed at crack tip, so the crack can not propagate. That makes me lost.

I realized that there probably is a numerical computing difficulty when \nu is close to 0.5 or K lead to infinity without knowing the details. Just from my calculation, I learnt that the big \nu is OK for calculation of displacements but seem not for stress. So anyone share some idea or just indicate some relative literature on that must very very appreciated.

Thanks very much for your kind replying. 

Best,
Lei   

[Prof. R.L. Taylor]

You are using a displacement element to perform nearly incompressible behavior.  Have you tried using the 'mixed' element that is designed to work for this case?   Elements must be quadrilateral or brick shape and can be bi-linear or
tri-linear for linear elements or bi-quadratic or tri-quadratic for quadratic elements.  Latest version may support cubic also.

[sanshi2000]

Thanks for your kind tip. I have used the 'mixed' option in the input file 'Inhkf' attached, but the principal stress at crack tip is still negative, opposite to the fact. At the same time, I have also try 9-node biquadratic quadrilateral element. No change occurred. Thanks very much for your kind considering.

Best,
Lei

[sanshi2000]

Now I guess the reason is that the coarse mesh can not recover the stress concentrate at the crack tip. I was just wondering whether I can generate nonuniform mesh in feap with what command. I want to make the mesh is very fine when close to crack tip and coarse far away from tip. Could you recommend some command or software?

By the way, you said I must use quadrilateral elements. Could you please indicator the reason as well? What is the disadvantage of triangle element?

Best,
Lei 

[FEAP_Admin]

Elements do not have to be quads; I think Prof. Taylor is referring to mixed elements.

In general, triangles are too stiff for mechanical problems and especially incompressibility.

If you want to grade your mesh you use the block generator move the mid side nodes of the
master block closer to one side -- but not closer than 1/4 of the side length (otherwise a singularity
occurs in the mapping function).  If you need a very carefully graded mesh then you should generate
the mesh externally and just read it into FEAP.

[sanshi2000]

Thanks a lot. Very helpful.   

Best,
Lei

[sanshi2000]

Hi. Any replying is reaaaaaly appreciated. 

1. I have been seeking the effect of \nu in modified neohookean material on principal stress around crack tip? Why does the stress become to negative if I change the \nu from 0.49 to 0.499 at line 20 in the attached file 'Inhkf'? I thought the material property is similar for \nu=0.49 or \nu=0.499, but the fem results is opposite. Any ideas on explaining that?

2. For the singularity of stress at tip, I realized that the fem with graded mesh or xfem are two right ways. For the xfem, I realized that more enriched terms are needed as shown in the equation (6) at p4 of the paper attached. I am just wondering whether there have been some options on dealing with that in feap if I know the form of \phi in (6)? Otherwise, what kind of change should I do in feap for considering the enriched displacement approach?

3. Also the paper indicates that 'quadrilateral quarter-point element' at Figure2 can remove one kind of singularity. So is there this kind of elements in feap?

4. Following your kind suggestion, I'd better use quadrilateral elements. But I realized that the quality of graded quadrilateral mesh is more difficult to control, so what difficulty am I going to meet if I use graded triangle or tetrahedral mesh?

5. Is there some examples in feap to consider the stress field around crack tips?

Thanks a lot.

[FEAP_Admin]

One thing to think about is do you really need to resolve the crack tip stresses super accurately?  Or in only an average sense.  If for example you just want crack tip driving forces, then J-integral theory will tell you that the having the local field super precise may be unnecessary.

Notwithstanding, if you want to try quarter point elements it is easy.  Place Quadratic elements around the crack tip.  Then for the elements that have a node on the crack tip move their midside nodes to the quarter points.  This will work with 6-node triangles and 9-node quads; with the 9-node quad move the center point to quarter-point too.

Also, at one time there was a Jintegral plot command in FEAP to plot driving forces on nodes [plot,jint].   It has not been used in a long time so it may no longer work but you could try this to see how good your meshes are (by using a standard problem with a known crack tip driving force).

[sanshi2000]

Thanks. I thought innocently the crack would growth if 'maximum principle stress at tip is beyond some value', which seem to be wrong and should be replaced by some average sense. Thanks very much for your kind advise on J-integration. But based on my initial reading, I realized that the J-integration, used to calculate energy released rate while the load is fixed, probably fails because the load in my problem is applied onto crack surfaces and also new generated crack surfaces.   

Also any average indicator of stress should depends on the stresses calculated in feap, so I still need to seek precise stress distribution by considering the problems 1,2,4,5 in previous post. Any ideas are reaaaaly appreciated.

Most importantly, is there some option on enriched displacement approach (xfem) in feap? Or is it easier to add them by feap user subroutines?

Thank your very much.

[Prof. S. Govindjee]

J-integrals apply to your case too.  The more general idea is to compute the energy release rate at the crack tip.  Have a look at the Tada Paris Shih handbook on the stress analysis of cracks.

[sanshi2000]

Thanks very much. I just employed a non-uniform triangular mesh and got a positive stress at crack tip.