Problems about outputing the state variables from USDFLD subroutine...help emergency

[JerryHe1]

Dear members:

I was trying to model the distribution of the temperature and degree of cure for a material. And as you can see from my codes,I used field(1) for the temperature and Fiedl(2) for the degree of cure,and each of them were stored in Statev(1) and Statev(2) respectively.But what confused me a lot is that the output for the SDV(2)--degree of cure, was always zero and didn't change with time.It seems like ABAQUS didn't call the subroutine at all. So can you help me check what's the problems with my subroutines? Thank you so much!!

Here is the codes:

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(*)
Real T_g,C_p0,C_p1,C_p,KK_1,KK_0,KK,E_1,k_1,k_2
PARAMETER(K_0=0.188)
PARAMETER(m=1.06,n=1.79,A_1=2.6e11,A_2=1.19e3,
1 E_1=1.24e4-273.15,E_2=5.7e3-273.15)
C
CALL GETVRM('TEMP',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,LACCFLA)
FIELD(1)=ARRAY(1)
STATEV(1)=FIELD(1)
IF (KINC.EQ.1) THEN
FIELD(2)=1E-4
FIELD(3)=0.0
ELSE

END IF
STATEV(2)=FIELD(2)
STATEV(3)=FIELD(3)
C CALCULATE THE GLASS TRANSFORMING TEMPERATURE (unit:C)
T_g=0.44*STATEV(2)*178/(1.0-0.56*STATEV(2))-42.
C
C
C CALCULATE THE SPECIFIC HEAT CAPACITY (unit:J/g*C)
C_p0=1.8500+2.625*1.0e-3*STATEV(1)
IF(STATEV(1).LT.T_g) THEN
C_p1=1.3125+4.437*1.0e-3*STATEV(1)
ELSE
C_p1=1.8500+2.265*e-3*STATEV(1)
END IF
C_p=(1-STATEV(2))*C_p0+STATEV(2)*C_p1
C
C CALCULATE THE COEFFICIENT OF HEAT CONDUCTION (unit: W/m*C)
KK_1=(3555.529*e-4)-2.727*e-4*STATEV(1)
KK_0=0.188
KK=(1-STATEV(2))*KK_0+STATEV(2)*KK_1
C
C THE KINETIC MODEL FOR CURING REACTION
k_1=A_1*EXP(-E_1/STATEV(1))
k_2=A_2*EXP(-E_2/STATEV(1))
STATEV(3)=(k_1+k_2*STATEV(2)**m)*(1-STATEV(2))**n
STATEV(2)=STATEV(2)+STATEV(3)*DTIME

RETURN
END


SUBROUTINE UMATHT(U,DUDT,DUDG,FLUX,DFDT,DFDG,
1 STATEV,TEMP,DTEMP,DTEMDX,TIME,DTIME,PREDEF,DPRED,
2 CMNAME,NTGRD,NSTATV,PROPS,NPROPS,COORDS,PNEWDT,
3 NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION DUDG(NTGRD),FLUX(NTGRD),DFDT(NTGRD),
1 DFDG(NTGRD,NTGRD),STATEV(NSTATV),DTEMDX(NTGRD),
2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3)
C
Real T_g,C_p0,C_p1,C_p,KK_1,KK_0,KK,E_1,k_1,k_2
PARAMETER(K_0=0.188)
PARAMETER(m=1.06,n=1.79,A_1=2.6e11,A_2=1.19e3,
1 E_1=1.24e4-273.15,E_2=5.7e3-273.15)
C
C CALCULATE THE GLASS TRANSFORMING TEMPERATURE (unit:C)
T_g=0.44*STATEV(2)*178/(1.0-0.56*STATEV(2))-42.
C
C
C CALCULATE THE SPECIFIC HEAT CAPACITY (unit:J/g*C)
C_p0=1.8500+2.625*1.0e-3*STATEV(1)
IF(STATEV(1).LT.T_g) THEN
C_p1=1.3125+4.437*1.0e-3*STATEV(1)
ELSE
C_p1=1.8500+2.265*e-3*STATEV(1)
END IF
C_p=(1-STATEV(2))*C_p0+STATEV(2)*C_p1
C
C CALCULATE THE COEFFICIENT OF HEAT CONDUCTION (unit: W/m*C)
KK_1=(3555.529*e-4)-2.727*e-4*STATEV(1)
KK_0=0.188
KK=(1-STATEV(2))*KK_0+STATEV(2)*KK_1
C
C THE KINETIC MODEL FOR CURING REACTION
k_1=A_1*EXP(-E_1/STATEV(1))
k_2=A_2*EXP(-E_2/STATEV(1))
STATEV(3)=(k_1+k_2*STATEV(2)**m)*(1-STATEV(2))**n
STATEV(2)=STATEV(2)+STATEV(3)*DTIME
C
DUDT=C_p
DU = DUDT*DTEMP
U = U+DU
C
DO I=1, NTGRD
FLUX(I) = -KK_1*DTEMDX(I)
DFDG(I,I) = -KK_1
ENDDO
C
RETURN
END

 

[akabarten]

Hi,

Does your *MATERIAL definition includes *USER DEFINED FIELD keyword to active field variables? If not USDFLD will be never call.
Did you allocate space for state dependent variables with *DEPVAR? If not state variable may be incorrect stored in the memory during calculation.
To check does Abaqus call USDLFD add following line to the subroutine:

write(6,*) '!!! USDFLD subroutine CALLED !!!'

After simulation see does the line exist in *.dat file, if not it is proof your USDFLD was never used.

You wrote STATEV(2) is always zero, what about others? Are they ok?

Another point is this part of yout code:

IF (KINC.EQ.1) THEN
 FIELD(2)=1E-4
 FIELD(3)=0.0
ELSE

END IF
 STATEV(2)=FIELD(2)
 STATEV(3)=FIELD(3)

It will set to unknown value STATEV(2) & STATEV(3) for all increment except first one since FIELD variables are not define previously.

Regards,
Bartosz

VIM filetype plugin for Abaqus

https://github.com/gradzikb/vim-abaqus

 

[JerryHe1]

Hi,Bartosz

Thank you for your help.I double check the input file and both the *USER DEFINED FIELD and *DEPVAR are there.I used the codes you showed me to check whether the USDFLD is called or not.And it turns out that all the subroutines are called by ABAQUS.So I think the real problem lies in the codes itself such as the following:

IF (KINC.EQ.1) THEN
FIELD(2)=1E-4
FIELD(3)=0.0
ELSE

END IF
STATEV(2)=FIELD(2)
STATEV(3)=FIELD(3)

My original idea was that for the first increment, the FIELD(2)and FIELD(3) are given initial values and then stored into STATEV(2) and STATEV(3) in USDFLD subroutine.After that,the USDFLD transfers the two STATEV() into another subroutine UMATHT,in which the STATEV() will be updated by this way:

STATEV(3)=(k_1+k_2*STATEV(2)**m)*(1-STATEV(2))**n
STATEV(2)=STATEV(2)+STATEV(3)*DTIME

So I think none of these STATEV() will be zero when outputing them.Is there something wrong with my idea or my codes?

Thanks again,
Ge

 

[akabarten]

Hi,

Is all STATEV are zero or only STATEV(2)?
Could you share your subroutine and model so I can run it by myself?

Regards,
Bartosz

VIM filetype plugin for Abaqus

https://github.com/gradzikb/vim-abaqus

 

[JerryHe1]

Hi,Bartosz

Here is my subroutine and input file,thank you so much for your help

Regards,
Ge

 

[JerryHe1]

and here is the input file.

Thanks

 

[akabarten]

Hi,

Thanks for the file, I ran them and here are my comments.

Following code is incorrect

**
IF (KINC.EQ.1) THEN
 FIELD(2)=1E-4
 FIELD(3)=0.0
 ELSE

 END IF
 STATEV(2)=FIELD(2)
 STATEV(3)=FIELD(3)
**

With this code STATEV(2) & STATEV(3) is always 0.0 expect 1st increment.
For increments higher than 1 FORTRAN does not execute IF and FIELD variables are not defined, they will be set to numerical garbage.
Very often set to 0.0 like in this case but it is not a rule. After this value of 0.0 is stored in your STATEV.

In your case there is no need to use FIELD variables at all, just use:

**
IF (KINC.EQ.1) THEN
 STATEV(2)=1E-4
 STATEV(3)=0.0
END IF
**

So I think none of these STATEV() will be zero when outputing them

And they are not, they are just very close to zero.
I took your eqation and make hand calculations:
1. For wrong SVARS(2)=0.0 and initial temperature SVARS(1)=20.0 I got 1.2e-212, with the highest temperature SVARS(1)=150.0 it is 2.0e-24
2. With correct SVARS(2)=1.0e-04 I got 9.8e-120 and 1.3e-17.

Values e-120 or e-242 is just zero for Abaqus, values of e-24 and higher is not zero and these values can be find in results.

Are you really expect such small values? Maybe the equation is incorrect or parameters used in it are wrong?

Regards,
Bartosz







VIM filetype plugin for Abaqus

https://github.com/gradzikb/vim-abaqus

 

[JerryHe1]

Hi,Bartosz

Sorry for replying late.I double check the codes and simplify the model to see what is going on with that.This time,I replace UMATHT subroutine with HETVAL and make a very simple codes like this:

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
IF (KINC.EQ.1) THEN
STATEV(1)=1E-4
C
ELSE
END IF
C
RETURN
END
C
C Boundary condition
SUBROUTINE FILM(H,SINK,TEMP,KSTEP,KINC,TIME,NOEL,NPT,
1 COORDS,JLTYP,FIELD,NFIELD,SNAME,NODE,AREA)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION H(2),TIME(2),COORDS(3), FIELD(NFIELD)
CHARACTER*80 SNAME
C
C
IF(TIME(2).LE.30.) THEN
SINK=20.+TIME(2)*58/30.
ELSE IF(TIME(2).LE.90.) THEN
SINK=78.
ELSE IF(TIME(2).LE.150.) THEN
SINK=78+(TIME(2)-90.)*12./60.
ELSE IF(TIME(2).LE.210.) THEN
SINK=90
ELSE IF(TIME(2).LE.240.) THEN
SINK=90+(TIME(2)-210.)*36./30.
ELSE
SINK=126
ENDIF
C
RETURN
END
C
SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,
1 PREDEF,DPRED)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
C
DIMENSION TEMP(2),STATEV(*),PREDEF(*),TIME(2),FLUX(2),
1 DPRED(*)
DOUBLE PRECISION m,n,det_E,A,Hr,R,pr

m=0.524
n=1.476
det_E=1.674e5 !(J/mol^-1)
A=3.7e22 !(Kg*m^-3)
Hr=7.7e4 !(kJ*kg^-1)
R=8.314
pr=1417 !(Kg/m^3) mass density for the resin matrix
C
C --------user coding to define FLUX and update STATEV------------
C
STATEV(2)=A*exp(-det_E/R/(TEMP(1)+273))*(STATEV(1)**m)*
&(1-STATEV(1))**n !rate equation for degree of cure(K=273+C)
C
STATEV(1)=STATEV(1)+STATEV(2)*DTIME
C
FLUX(1)=pr*Hr*STATEV(2)
C
RETURN
END

However,the most weird thing is that all SDVs are zero during each increment including FLUX(1),but I made a hand calculation:
for the first increment,m=0.524,n=1.476,det_E=1.674e5,A=3.7e22,R=8.314,STATEV(1)=1e-4,TEMP(1)=20
STATEV(2)=A*exp(-det_E/R/(TEMP(1)+273))*(STATEV(1)**m)*(1-STATEV(1))**n=9.5121e-006
what's more,even if I set STATEV(1)=1E-4 at first increment in USDFLD,the SDV(1) still be zero.

so can you help me check the input and codes again,I appreciate it a lot.

 

[JerryHe1]

and here is the input file.Thanks!

 

[akabarten]

Hi Jerry,

USDLFD subroutine is never used in your example.
I added *FIELD, USER, VARIABLE=2 to fix it.

Variable KINC starts from 0 so I changed condition to set STATEV(1) to

IF (KINC.EQ.0) THEN
  STATEV(1)=1E-4
END IF

After these changes I got STATV(2)=4.2441E-010 and it fits for my hand calculations according to equation in the subroutine.

m=0.524,n=1.476,det_E=1.674e5,A=3.7e22,R=8.314,STATEV(1)=1e-4,TEMP(1)=20
STATEV(2)=A*exp(-det_E/R/(TEMP(1)+273))*(STATEV(1)**m)*(1-STATEV(1))**n=9.5121e-006

Please check your calculation.

Regards,
Bartosz

VIM filetype plugin for Abaqus

https://github.com/gradzikb/vim-abaqus

 

[JerryHe1]

Hi Bartosz,

Did the STATEV(1) and STATEV(2) have a value at each increment in your subroutine? I did as what you said and added "*FIELD, USER, VARIABLE=2" there:

** STEP: Step-1
**
*Step, name=Step-1, nlgeom=NO
*Heat Transfer, end=PERIOD, deltmx=5.
1., 250., 0.0025, 5.,
**
[highlight #EF2929]*FIELD,USER,VARIABLE=2[/highlight]
** INTERACTIONS

but nothing happened except these values still being zero at each increment.

Can you send me your revised input file and codes please?

Thanks!

 

[akabarten]

Hi Jerry,

I tested only for datatcheck and now I see SDV values are non zero only in first increment.
I check it once again and I found that STATEV values in USDFLD are reset to zero between KINC=0 and KINC=1, later on they are kept in memory as it should be.
I do not know why it works in that way.

I am attaching my latest model + subroutine. Looks it works for all increment now.

Regards,
Bartosz

VIM filetype plugin for Abaqus

https://github.com/gradzikb/vim-abaqus