Incompressible user material

[Sebastian Castro]

Hi,

I'm writting an incompressible user material. First of all, I wrote an compressible hyperelastic material (using the umati and umatl functions) and it works. Then, I want to add incompressibility using the augmented lagrangian option (AUGM). I found a topic from 2012 (http://feap.berkeley.edu/forum/index.php?topic=91.0) where the same question was formulated.

Based on that topic, I wrote a function which is called by umodelf with the inputs

Code: [Select]

(f,finv,df,detf,be,ta,d,ud,hn,h1,nh,ntm,ii,istrt, sig,dd, xlamd,ha, isw)
This function returns the cauchy stresses and the spatial tangent moduli, adding the incompressibility term, as the mixed formulation says (I took as reference the book Nonlinear Solid Mechanics A Continuum Approach for Engineering by Gerhard A. Holzapfel), i.e. using the deviatoric part of the gradient deformation. However, it's not working . Now I list all the questions that I have:

• How do I choose the AUGM's parameters? If I run the script attached at the end of this post, the first augment iteration is done correctly (the solution converges). However, in the second iteration, when I solve the problem with UTAN, the energy does not converge (it says "indefinite matrix").
• Do I have to compute something else than the cauchy stresses and the tangent spatial moduli? I don't know if the AUGM function uses something that I'm not calculating and saving. For example, in neoh3f you compute the energy density, it is necessary for an user material?
• I have two type of elements, a solid element with my user material and a pressure element. Could that affect in the augmented lagrangian iteration? I don't believe.
• How do I determine a value for d(21) if I don't have an elastic modulus and either a Poisson's ratio in my material?

Any help is really appreciated.

Thanks

Sebastian

EDIT

I was reading the modlfd.f file and I noticed that there is a difference when FEAP calls umodel and umodelf:

Code: [Select]

c     U s e r    M o d e l    I n t e r f a c e                                                                                                                                                                                                                       

      elseif(umat.le.100) then

elamd = xlamd
eha   = ha
call umodel(umat,f,detf(1),ta,d(1),d(uprm+1),hn(1),hn1(1),nh,
     &              ii,istrt, sig,aa, isw)
xlamd = elamd
ha    = eha

c     U s e r    F i n i t e    M o d e l    I n t e r f a c e                                                                                                                                                                                                         

      else

umat = umat - 100
call umodelf(umat,f,finv,df,detf(1),b2,ta,d(1),d(uprm+1),
     &               hn(1),hn1(1),nh,ntm,ii,istrt, sig,aa,xlamd,ha,isw)

      endif
I added the same two lines before and after umodelf (i.e. "elamd = xlamd, eha   = ha" and "   xlamd = elamd,   ha    = eha") as in the umodel case and now it doesn't crush at the second iteration. It was a bug or those lines shouldn't be there?

Bye

Code: [Select]

FEAP ** 3D bar
0,0,0,3,3,8

PARAmeters
lx = 10.0
ly =  1.0
lz =  1.0
mx = 50
my = 5
mz = 5

BLOCk 1 ! Beam
  CARTesian,mx,my,mz,0,0,1,10
  1 0.0 0.0 0.0
  2 lx  0.0 0.0
  3 lx  ly  0.0
  4 0.0 ly  0.0
  5 0.0 0.0 lz
  6 lx  0.0 lz
  7 lx  ly  lz
  8 0.0 ly  lz

BLOCk 2   ! Pressure Elements
  CARTesian,mx,my,0,0,2,0
  1 0.0 0.0 0.0
  2 lx  0.0 0.0
  3 lx  ly  0.0
  4 0.0 ly  0.0

EBOU
  1, 0., 1, 1, 1

EDIS
  1, 0., 0., 0., 0.

MATErial,1
  SOLID
  FINIte
  MIXE
  UCON gucc 2.0 8.0 2.0 4.0 1e12
  AUGM ON

MATErial,2
  PRESsure
  LOAD 0.0000004 ! very low load
  FOLLower
  AUGM OFF

END

TIE
OPTI

BATCh
  LOOP augment 3
    AUGM,,1e7               ! I don't know how it works really
    LOOP,,100
      UTANG,,1
    NEXT
  NEXT augment
  DISP,node
  STRE,node
  PVIE,time
END

STOP


[FEAP_Admin]

Some partial answers.

(1) You do not need to use umodelf if you need augmentation.  The lines around the call to umodel take the augmentation parameters and copy them to elamd and eha which are stored in the common block elaugm.h.  You can reference elaugm.h in your user model to have access to these parameters.  They will needed to compute the pressure properly (have a look at neoh3f.f).

(2) You need to compute the pressure and the pressure tangent taking into account the augmentation; see neoh3f.f and fengy3.f (which computes among other things the derivatives hp and hpp).

(3) should be ok -- but not entirely certain.  Replace the pressure element with a simple dead load to check.

(4) d(21) is a computed quantity IF you use FEAP's built-in machinery for property input.  Since you do not have E or nu you have two choices.  (a) fake it by entering values of E and nu that result in the desired K value or (b) add lines in you property input file to input K and store it in ud( ).

Following what you see in neoh3f.f indicates all the items you need to compute in your material model.  Note that FEAP use h(J) = J-1 as the augmentation function and the Uzawa update for the lagrange multiplier.  See the paper of Simo and Taylor 1991 in CMAME but note that FEAP does not require the split energy format the user model is passed the mixed deformation gradient F_mixed = F*(theta/J)^1/3.   Also assume augf = 1, the default, when trying to follow the code.

Note some of these thing were not originally intended for users to access so there may still be hiccups.

[Sebastian Castro]

Thank you for your answer. I have been reading the codes of neoh3f, modlfd and others in order to understand well what are the inputs.

This is an idea of what I'm doing in my material.

Code: [Select]

      subroutine umatl3( f , detf,td,d,ud,hn,h1,nh,ii,istrt, sig,dd,isw)

      implicit none

      include   'pmod2d.h'      ! etype = 1 (displacement method)                                                                                                                                                                                                       
                                !       = 2 (mixed method)                                                                                                                                                                                                             
      include   'elaugm.h'      ! elamd, eha                                                                                                                                                                                                                           
      include   'elengy.h'      ! estore                                                                                                                                                                                                                               
      include   'pconstant.h'   ! program constants                                                                                                                                                                                                                     

      if (isw.eq.14) then       ! initialize history terms                                                                                                                                                                                                             

      else

! This function computes the cauchy stresses, the spatial tangent moduli                                                                                                                                                                                       
        ! and the energy density                                                                                                                                                                                                                                       

         call guccione(f,detf,d,ud, sige,dde,estore)

         if(etype.eq.1) then

            ! Write sig and dd considering sige and dde,
            ! in this case sig = sige                                                                                                                                                                                                                                   
            !              dd  = dde

         elseif(etype.eq.2) then

c     Augementd Lagrangian parameters                                                                                                                                                                                                                                   
            ha    = eha
            xlamd = elamd

c     Compute pressure and volumetric moduli                                                                                                                                                                                                                           
            jsw = nint(d(170))                                                                                                                                                                                                           
            call fengy3(d,detf(3), u,up,upp,ha,hp,hpp,jsw)

c     Pressure and tanget (not mixed pressure)                                                                                                                                                                                                                         
            press =  up  + xlamd * hp
            upp   = (upp + xlamd * hpp ) * detf(1)

            ! Write sig and dd, following the paper of Simo and Taylor 1991,                                                                                                                                                                                       
            ! considering sige, dde and the pressure field                                                                                                                                                                                                                           

c     Correct stored energy                                                                                                                                                                                                                                             
            estore = estore + u

         endif                  ! etype                                                                                                                                                                                                                                 
      endif                     ! isw case                                                                                                                                                                                                                             
      end                       ! program

Now, in the input file, I define the materials as follows

Code: [Select]

MATErial,1
  SOLID
  FINIte
  INCOmpressible
  MIXEd
  UCON gucc 2.0 8.0 2.0 4.0 1
  AUGM ON

MATErial,2
  PRESsure
  LOAD 1 1
  FOLLower
  AUGM OFF

And the batch mode

Code: [Select]

BATCh
  PROP,,1
  dt,,0.1
  LOOP,,2
    TIME
    LOOP augment 3
      AUGM,,1.0
      LOOP,,100
UTANG,,1
      NEXT
    NEXT augment
    DISP,node
    STRE,node
    PVIE,time
  NEXT
END
1 0 0.0 100 0.0 0.000004 0.0 0.0

However, I'm still having problems

• Is the general idea of the user material right? Or am I forgetting something? Because when I run my example, I noticed that xlamd is always equal to zero.
• In the batch mode, when I use AUGM, is it well used? I wonder if the problem with xlamd has a relation with this.
• Is there a function in FEAP to output the volume (or the determinant of the gradient deformation) of an element? I know I can write my own function, but I'd like to know if FEAP has one.

Thanks again for your time, I really appreciate it.

[Sebastian Castro]

Well, I found a mistake in my code. I forgot to add two lines after to call fengy3. Now both xlamd and ha are changing along the solution.

Code: [Select]

      subroutine umatl3( f , detf,td,d,ud,hn,h1,nh,ii,istrt, sig,dd,isw)

      implicit none

      include   'pmod2d.h'      ! etype = 1 (displacement method)                                                                                                                                                                                                       
                                !       = 2 (mixed method)                                                                                                                                                                                                             
      include   'elaugm.h'      ! elamd, eha                                                                                                                                                                                                                           
      include   'elengy.h'      ! estore                                                                                                                                                                                                                               
      include   'pconstant.h'   ! program constants                                                                                                                                                                                                                     

      if (isw.eq.14) then       ! initialize history terms                                                                                                                                                                                                             

      else

! This function computes the cauchy stresses, the spatial tangent moduli                                                                                                                                                                                       
        ! and the energy density                                                                                                                                                                                                                                       

         call guccione(f,detf,d,ud, sige,dde,estore)

         if(etype.eq.1) then

            ! Write sig and dd considering sige and dde,
            ! in this case sig = sige                                                                                                                                                                                                                                   
            !              dd  = dde

         elseif(etype.eq.2) then

c     Augmented Lagrangian parameters                                                                                                                                                                                                                                   
            ha    = eha
            xlamd = elamd

c     Compute pressure and volumetric moduli                                                                                                                                                                                                                           
            jsw = nint(d(170))                                                                                                                                                                                                           
            call fengy3(d,detf(3), u,up,upp,ha,hp,hpp,jsw)

c     Save current Augmented Lagrangian parameters                                                                                                                                                                                                                                   
            eha   = ha
            elamd = xlamd

c     Pressure and tanget (not mixed pressure)                                                                                                                                                                                                                         
            press =  up  + xlamd * hp
            upp   = (upp + xlamd * hpp ) * detf(1)

            ! Write sig and dd, following the paper of Simo and Taylor 1991,                                                                                                                                                                                       
            ! considering sige, dde and the pressure field                                                                                                                                                                                                                           

c     Correct stored energy                                                                                                                                                                                                                                             
            estore = estore + u

         endif                  ! etype                                                                                                                                                                                                                                 
      endif                     ! isw case                                                                                                                                                                                                                             
      end                       ! program
Now my problem is with the AUGM option. If I run my script, the first augment loop converges, but in the second step, when is solving the system with a new xlamd, it diverges. I have been playing with the inputs but I don't get any successful result.

Thank you.

[Prof. R.L. Taylor]

can you post your user modules? 

[Sebastian Castro]

I'm attaching my user material and the input file.

I found where was the problem. I was using a bulk modulus equal to 1 because I though that it shouldn't affect in my case, but I noticed that for the incompressible case you set the bulk modulus as 0. I changed the value to zero and now I don't have that problem. However, setting the bulk modulus as zero, xlamd remains always as zero too. I'm curious about this, why does not this parameter change? I changed the material of my script for a modified Neohookean (E = 2000, nu = 0.499) with the same options for AUGM and xlamd always had a value of zero (I'm setting the INCOmpressible option to be sure).

So, I don't understand why is this happening. Where do you add the pressure field if xlamd is always zero?

Thanks

[Prof. R.L. Taylor]

The idea of augmented Lagrangian is to do a slightly compressible material and then enforce incompressibility as an iterative constraint.  Thus, you should not use the INCOmpressible option.  This was inserted to allow us to do a u-p formulation (nodal p) where one can do the incompressible option.  If you remove the INCO then you can check if augmenting does enforce the incompressibility correctly.  You do still nedd the MIXEd option.

[Sebastian Castro]

Thanks for your explanation. I removed the INCO option and now it's working. I checked the value of J in every time step and it was equal to 1 or almost (0.99999999999989 ).

A few questions:

• Does the election of U(J) affect in the displacement field when J=1?
• How can I get the pressure field as an output? Should I save the press variable in my material and write a file later?
• Does the pressure element work in parFEAP? At least it's not working for me. (This may be in other subforum I think)

Thank you very much!

EDIT

About question 3, I noticed the problem wasn't the pressure element but the mixed formulation. If I use MIXEd, even with the modified Neohookean model, the displacement is always zero. I tried with EFORce instead the pressure element and the result was the same, but when I quit the mixed option it worked well with both cases (eforce and pressure element). Can't I use mixed formulation in parFEAP?

[Prof. R.L. Taylor]

1.) No
2.) The projected Lagrange multiplier after you converge J=1 is the pressure
3.) What do you mean it does not work?

 Look at how augmented Lagrangian is performed in fld2d2.f  If you are projecting at the user material model level  that is not consistent with the mixed model.  There is a projector on the stresses (phi) for each element, not each Gauss point.

[Sebastian Castro]

Finally I found what was happening. I had defined augf = 1 in the serial case and it was working well, but in the parallel case I had to change the augf value. With the new value (1e-5 in my case) both cases are solved fine. I used the modified Neohookean model to verify this.

Thank you for all your help professor Taylor and FEAP_Admin.