Variation of Modulus of Elasticity, E with mean stress using USDFLD in ABAQUS 1

[johniftekhar]

Hi,
I am new here. I wrote a USDFLD subroutine to change young's modulus based on the mean stress. First, I tried to keep the E value constant with field(1) then I am tring to vary the E with mean stress. In usdfld file, I gave constant field(1)=60000, but when it is running input file, initially it taking E value less than 60000. Why it is taking different E value as I keep the field(1) value constant.

Input file:

*HEADING
Mohr-Coulomb, CAX4 element, Triaxial Compression
E=A(p)^N is used using usdfld fortran file
Inputs in USDFLD.for are C, P_ref,N, and minimum p
*NODE,NSET=ALLN
1, 0., 0.
2, 1., 0.
3, 1., 1.
4, 0., 1.
*ELEMENT,TYPE=CAX4,ELSET=ALLE
1,1,2,3,4
*SOLID SECTION,ELSET=ALLE,MATERIAL=ALLE
*MATERIAL,NAME=ALLE
*Density
1.77,
*ELASTIC, TYPE=ISOTROPIC, DEPENDENCIES=1
100.,0.3,0.,100.
300000.,0.3,0.,300000.
*USER DEFINED FIELD
*DEPVAR
5
***ELASTIC
**160000.,0.3
***Mohr Coulomb
**44.,16.
***Mohr Coulomb Hardening
** 0.5,0.
*BOUNDARY
1,PINNED
2,2
4,1
*INITIAL CONDITIONS,TYPE=STRESS
ALLE,-60.,-60.,-60.
*INITIAL CONDITIONS,TYPE=Field, VARIABLE=1
ALLN,70000.
*STEP,INC=1000
*GEOSTATIC
1.,1.
*DLOAD
1,P2,60.
1,P3,60.
*EL PRINT
S,
SINV,
E,
PE,
SDV,
FV
**EE,
*NODE PRINT, TOTALS=yes
U2, RF2
*EL FILE,FREQ=10
S,E,PE
*END STEP
*STEP,INC=20000
*Static
0.0001, 1, 1e-15, 1.
***STATIC,DIRECT
**1.,20.
*BOUNDARY
3,2,,-.2
4,2,,-.2
*Output, field, variable=PRESELECT
*Element Output, directions=YES
LE, PE, PEEQ, S, SDV, FV
*END STEP


usdfld file:


SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT,
1 TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,
2 KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO,
3 LACCFLA)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME,ORNAME
CHARACTER*3 FLGRAY(15)
DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3),
1 T(3,3),TIME(2)
DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),
1 COORD(*)

C Get values of current stress:

CALL GETVRM('S',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,
1 MATLAYO,LACCFLA)
C ARRAY(1)=S11, ARRAY(2)=S22, ARRAY(3)=S33,ARRAY(4)=S12,ARRAY(4)=S13,ARRAY(4)=S23
c write(*,*)array(1),array(2),array(3)

c if (time(1) .GE. 1.0)then
C MODULUS OF ELASTICITY DEFINED AS E=C*(P/P_REF)^AN
!INPUT
C=100000.0D0
P_REF=100.0D0
AN=0.5D0
P_min=1.0D0 !kPa, for numerical issue
P_max=1000.0D0 !kPa, for numerical issue
!END OF INPUT FOR ELASTICITY

P_current=(array(1)+array(2)+array(3))/(-3.0D0)
IF (P_current .LE. P_min)P_current=P_min
IF (P_current .GE. P_max)P_current=P_max
E_current = C*(P_current/P_REF)**AN
if (time(2) .ge. 1.0)then
write(6,*) time(1),time(2), field(1)
field(1) = 60000.*time(2) !E_current
else
field(1) = 60000.
endif
sig22= array(2)
c q=array(2)-array(3)
c statev(1)=p_current
c statev(2)= E_current
C Get value of current strain:

CALL GETVRM('E',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,
1 MATLAYO,LACCFLA)
C ARRAY(1)=E11, ARRAY(2)=E22, ARRAY(3)=E33,ARRAY(4)=E12,ARRAY(4)=E13,ARRAY(4)=E23
strain22=array(2)
c write(6,*)'test',P_current, field(1)
c write(6,*)'test',vertical_strain, q, field(1)
C If error, write comment to .DAT file:

if (time(2) .ge. 0.1)then

E_modulus=(sig22-statev(1))/(strain22-statev(2))

print*,field(1),E_modulus,time(2)
endif


statev(1)=sig22
statev(2)=strain22

IF(JRCD.NE.0)THEN
WRITE(6,*) 'REQUEST ERROR IN USDFLD FOR ELEMENT NUMBER ',
1 NOEL,'INTEGRATION POINT NUMBER ',NPT
ENDIF


C
RETURN
END

 

[YhTian76]

Hi, John
Have not run your code and input file. Some comments for you to consider:
1. I guess your initial E should be 100. Check it
2. I ever used USDFLD and I found the FV is "explicit" - it becomes effective at next increment when you give it a new value.
3. Try to add 60000.,0.3,0.,0. after your *elastic to whether you can get what you want or not.
yt

 

[akabarten]

Hi,

I ran your model and in first increment I got FV1 value of 70000 (value set with *INITIAL CONDITIONS keyword) and in other increments Abaqus keeps constant value of 60000.

initially it taking E value less than 60000

What value did you get?

Regards,
Bartosz

 

[IceBreakerSours]

Here's what I got by running the code:

FV1:

60000
60000
60006
60012
60021
60034.5
60054.8
60085.1
60130.7
60199
60301.5
60455.3
60686
61032
61551
62329.4
63497.2
65248.7
67876.1
71817.1
77728.7
86596
99897.1
119849

http://www.eng-tips.com/faqs.cfm?fid=376

 

[IceBreakerSours]

The first set of values for E is actually close to "zero" and not just less than what you want it to be.

http://www.eng-tips.com/faqs.cfm?fid=376

 

[akabarten]

Hi,

Thanks for your updated and sorry for my late response.
I found some time to take a look for your model and the subroutine once again.

Please find my comments.

Why it is taking different E value as I keep the field(1) value constant.

You keep FV1 constant only for first step. You have two steps, first *GEOSTATIC with time 1.0 and second *STATIC with 1.0 time.
To assign field(1) variable you are using:

if (time(2) .ge. 1.0)then
write(6,*) time(1),time(2), field(1)
field(1) = 60000.*time(2) !E_current
else
field(1) = 60000. 
endif

The is statement is true for second step and false for first one.
It means for first step you set constant value but for second you set field(1) = 60000.*time(2).

The Results you post come from this equation (field(1) = 60000.*time(2)).


Regards,
Bartosz

 

[IceBreakerSours]

Hi Bartosz,

I may be wrong but I think johniftekhar's concern was the small initial value of the Young's modulus. And, in my opinion, its only because the first iteration of the first increment has no "previous" incremental strain and stress values which is why the Young's modulus comes out as nearly zero. Is that right? Also, I "think" ABAQUS calls the subroutine twice in the first iteration of every increment which explains values for, say Modulus, showing up 8 times for a 4 noded element. Job diagnostics show that every increment was completed in one iteration.

I'd appreciate your comments.

http://www.eng-tips.com/faqs.cfm?fid=376

 

[akabarten]

Hi,

johniftekhar's concern was the small initial value of the Young's modulus

I am not sure, as I see the first value is 60000.

first iteration of the first increment has no "previous" incremental strain and stress values which is why the Young's modulus comes out as nearly zero

I agree with you, in first increment the FV values should be 0.0 ... but not when you set initial conditions and in the inputdeck we have

*INITIAL CONDITIONS,TYPE=Field, VARIABLE=1
ALLN,70000.

I did following, remove if condition and set only

field(1)=60000

And I got FV=70000 for first subroutine call and 60000 for all other calls, both single and double precision.

I "think" ABAQUS calls the subroutine twice in the first iteration

I think the same.

Regards,
Bartosz

 

[IceBreakerSours]

I see the first value is 60000.

That's odd because the first 8 values for modulus that I see in the command prompt are of the order of e-10. The following piece of code displays it:

if (time(2) .ge. 0.1d0) then
  E_modulus=(sig22-statev(1))/(strain22-statev(2))
  print*,field(1),E_modulus,time(2)
endif

Also, the field(1) values start from 60000 in the command prompt. On the other hand, the field(1) values in the .dat file stays constant at 70000.

http://www.eng-tips.com/faqs.cfm?fid=376

 

[akabarten]

Hi,

I see in the command prompt are of the order of e-10.
Also, the field(1) values start from 60000 in the command prompt.

I am a little bit confus, so what value do you have?

On the other hand, the field(1) values in the .dat file stays constant at 70000.

It because the line

write(6,*) time(1),time(2), field(1)

is before you set field(1) value.

The field array does not have values from previous increment.
When Abaqus passes the array into the subroutine it is always initialize to 0.0 or value from initial condition.
The array should be only write not read. If for some reason we need to know field values from previous increment we have to save them into STATEV array.

Regards,
Bartosz

 

[IceBreakerSours]

Here are the first 8 values of field(1), E_modulus, time(2) that show up in the command prompt:

field (1) E_modulus time(2)
60000.0000000000 -8.331113576788606E-009 1.00000000000000
60000.0000000000 -8.704148513061304E-010 1.00000000000000
60000.0000000000 -3.854694341498849E-009 1.00000000000000
60000.0000000000 3.606004383982206E-009 1.00000000000000
60000.0000000000 -8.331113576788606E-009 1.00000000000000
60000.0000000000 -8.704148513061304E-010 1.00000000000000
60000.0000000000 -3.854694341498849E-009 1.00000000000000
60000.0000000000 3.606004383982206E-009 1.00000000000000

When I switch the lines as per your suggestion, here's what I get:

60000.0000000000 55678.1938255410 1.00000000000000
60000.0000000000 73324.8836068822 1.00000000000000
60000.0000000000 57955.0168574356 1.00000000000000
60000.0000000000 50547.7375126345 1.00000000000000
60000.0000000000 55678.1938255410 1.00000000000000
60000.0000000000 73324.8836068822 1.00000000000000
60000.0000000000 57955.0168574356 1.00000000000000
60000.0000000000 50547.7375126345 1.00000000000000

Also, in the .dat file, field(1) starts from 60000. Shouldn't it have started from 70000 in the first step? Perhaps, there is something wrong with what I did. INP/subroutine are attached.

http://www.eng-tips.com/faqs.cfm?fid=376

 

[akabarten]

Hi,

Shouldn't it have started from 70000 in the first step?

You are right however you do not write field value for first step.
First step has only one increment with value of 1.0.

So total time (time(2) variable) is changing as follow:
0.0
1.0
1.0 + delta t from automatic incrementation for step 2
...
2.0

We have two condition to control write and print statement:
- first

IF (time(2) .ge. 1.0) then
field(1) = 60000.*time(2) 
write(6,*) time(1),time(2), field(1)
end if

- second

if (time(2) .ge. 0.1d0) then
print*,field(1),E_modulus,time(2)
endif

For time(2) = 0.0 both conditions are not true (we do not write/print field value however it's 70000).
For time(2) = 1.0 first condition is true so we assign field(1)=60000 befor second condition.

As you can see there is no mid-time value between 0.0 and 1.0.

Regards,
Bartosz