FEAP User Forum
FEAP => Programming => Topic started by: Sebastian Castro on December 12, 2014, 12:22:29 PM
-
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 (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
(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:
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
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
-
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.
-
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.
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
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
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.
-
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.
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.
-
can you post your user modules?
-
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
-
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.
-
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 :P).
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?
-
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.
-
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.