bizarre results for function of young's modulus with principal stress defined by USDFLD 1

[smallpea14]

Hi all,

I got great help from reading previous posts relate to defining elasticity function by USDFLD. By imitating those examples, I defined young's modulus as a simple function of minor principal stress by following two methods. Young's modulus was defined as field variable by itself in method a while minor principal stress was defined as filed variable. I run a simple triaxial test by applying 10kpa confining stress and downward displacement of 0.3 m for a asymmetric block of 1m. I got bizarre stress after the displacement reached to 0.15m by method b and vary different stress at the end of calculation by two methods. Could anyone please give me some hints what is the reason for this? I don't really have any clue. Thanks so much in advance.
method a:
PARAMETER (PA=101.3D0,EK=300.0D0,EN=0.5D0,EMIN=50.0D0,EC=250.D0)
CALL GETVRM('SP',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,
& LACCFLA)
SP1=ARRAY(1)
SP2=ARRAY(2)
SP3=ARRAY(3)
WRITE (6,*) SP1,SP2,SP3
IF (SP3.EQ.0.0) EJ=EC
IF (SP3.GT.0.0) EJ=EJMIN
IF (SP3.LT.0.0) EJ=EK*PA*(-SP3/PA)**EN
FIELD(1)=EJ
STATEV(1)=EJ
*Elastic, dependencies=1
10., 0.3, , 10.
1e+06, 0.3, , 1e+06
method b:
PARAMETER (PA=101.3D0,EK=300.0D0,EN=0.5D0,EMIN=50.0D0,EC=250.D0)
CALL GETVRM('SP',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,
& LACCFLA)
SP1=ARRAY(1)
SP2=ARRAY(2)
SP3=ARRAY(3)
WRITE (6,*) SP1,SP2,SP3
IF (SP3.LT.0) THEN
FIELD(1)=-SP3
STATEV(1)=EK*PA*(-SP3/PA)**EN
ELSE
FIELD(1)=0
STATEV(1)=EC
END IF
*Elastic, dependencies=1
250., 0.3, , 0.
9548.3, 0.3, , 10.

 

[akabarten]

Hi,

Do you expect to have the same results for both method, because you will never get this.

To set E you need to set field variable value with FIELD(1) array.
In second approach you set it always to -SP3 or to 0.0 value.
But you set STATEV array to your material law which have no influence for field variable in your case at all.
What is a reason of using STATEV array in your case?
Perhaps you misunderstanding difference between field variable and state-dependent variable?

Regards,
Bartosz

 

[smallpea14]

Hi Bartosz,

Thanks so much for your reply, Bartosz. Actually, I got the idea from your previous posts answering other person's questions.
As for method b, yes, you are correct.I set field variable as -sp3 or zero but in the input file i calculated the E as the same function as defined in method A, ie., E=300*101.3*(-SP3/101.3)**0.5 (-sp3=10 in this simple triaxial test). in my case, E is only function of compression (it's positive for my case) while E is constant when tension. Why I set sdv was that I just wanna to see if sp3 was the same for both cases in odb file unfortunately they were different. sp3 for both cases differed after several increments and were no longer the -10 kpa supposedly. I also tried the simplest case without any subroutines, ie., i set E as constant as 9548.3 and got very different stress of mises with those from method a and b both. Could you please give me more comments? I really appreciate for your help. I attached the simple inp file which takes less than 2 mins.
Regards
Muller

 

[akabarten]

Hi,

I made two runs of your model without material dependencies.
One with constant E=250.0 and second one with E=9548.0.
There is one thing which I am not sure. For both runs S33(SP3)=-10.0.
Does it mean your E law which is a function of S33(SP3 variable in the subroutine) and it never changes?
Is it ok for you?

Since SP3 does not change when I changed E value I made two runs with sub which should give me the same results.
One where I always set constant field variable value (FIELD(1)=10.d). That runs is stable.
Next with FIELD(1)=-array(3) and here are very strange things. Because here SP3(S33) is not -10.0.
It various with many different values, strange values. I expected that SP3(S33) will be always -10.0 like in runs without material dependencies.

What do you think about it?

Regards,
Bartosz

 

[smallpea14]

Hi Bartosz,

Thanks so much for your efforts. Yes, I need to define a E law dependent on SP3(minor principal stress) with Mohr-Couloumb plasticity failure criterion eventually. In my case, sp3 should increase gradually. This very simple triaxial test was just meant to test if this simple usdfld should work as I expected, an exponential function of sp3. I had the same issue as yours, sp3 is not stable as it should be when without material dependencies. So, i am not sure if the problem is due to my inp file or just the problem of usdfld. Ideally, the results should be the same for both methods, either defining E as a field variable or defining E dependent sp3 as a field variable. Actually, I wrote a umat for this E law, which worked well. But I need this E law coupled with Mohr-Couloumb plasticity. In that case, umat is too complex when incorporating Mohr-coloumb. I think there should be an easier way that I just define elastic part by usdfld then I add Mohr-coloumb part available already in ABaqus. I imitated your previous postings and wrote the usdlf above. Do you think if this idea is applicable for my case?
I don't really understand how usdfld called by Abaqus. As show in this linke:

http://www.eng-tips.com/viewthread.cfm?qid=295995

,
You once prepared an example of E and commented that this usdfld only worked with reduced element rather than full-integration element. I can't really understand this comment. Does this issue apply to my case? Or this is due to time increment since usdfld is called at the beginning of current time increment. Could you please give me more comments?
I also tried to define permeability as a field variable based on the same idea. Again, it differs with the result for permeability dependent on Void ratio. I used the Abaqus example, Terzaghi 1D consolidation. Firstly, I defined permeability dependent on void ratio as below:
*Permeability, specific=1.
2e-05, 1.
2.24404e-05, 1.025
2.51785e-05, 1.05
2.82508e-05, 1.075
3.16979e-05, 1.1
3.55656e-05, 1.125
3.99052e-05, 1.15
4.47744e-05, 1.175
5.02377e-05, 1.2
5.63677e-05, 1.225
6.32456e-05, 1.25
7.09627e-05, 1.275
7.96214e-05, 1.3
8.93367e-05, 1.325
0.000100237, 1.35
0.000112468, 1.375
0.000126191, 1.4
0.000141589, 1.425
0.000158866, 1.45
0.00017825, 1.475
0.0002, 1.5
This is basically the linear approximation of the function defined as in usdfld
PARAMETER (PK0=0.0002D0,EV0=1.5D0,PCK=0.21714724D0)
CALL GETVRM('VOIDR',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,
& LACCFLA)
EV=ARRAY(1)
STATEV(1)=EV
PK=PK0*EXP((EV-EV0)/PCK)
FIELD(1)=PK
STATEV(2)=PK
*Permeability, specific=1, DEPENDENCIES=1
2e-05, 1, , 2e-5
0.0002, 1.5, , 2e-4
*User Defined Field
*Depvar
2,
Did I make the same mistake as I defined E law? I really don't have any clue. Could you please give me more comments, Bartosz? I really appreciate for your help.

 

[akabarten]

Hi,

Thanks for the explanation.

I don't really understand how usdfld called by Abaqus. As show in this linke:

http://www.eng-tips.com/viewthread.cfm?qid=295995

This is a different case as yours. USDDFLD subroutine were used to distribute material properties base on position.
Material properties are constant and do not change during calculation so they are set only once at first increment.

reduced element rather than full-integration element, Does this issue apply to my case?[quote

I do not think, I tested your model also with CAX4 and CAX4R elements and the model is also not stable at the end.
You can try do not use axisymmetric elements, I think.

time increment since usdfld is called at the beginning of current time increment

The point is that USDFLD has access for results from n-1 increment in n increment where usually abaqus has results
from n increment in n increment. In other way you have results from previous increment where in perfect case you should have results from current increment. When you compare two models (one with USDFLD and one without) results can different slightly. It's also a good idea to not compare the same increments but n increment in one model with n-1 increment in other model. The difference decrease when you cut down your time step since the difference between n and n-1 increment is small. I tested your model with 0.1 increment and it is still not stable.

So, i am not sure if the problem is due to my inp file or just the problem of usdfld.

I did not find anything strange in your model I would did it the same.
Maybe some Abaqus bug? I would sent the model to Simulia support if you can.

I am attaching nice document about user subroutines in Abaqus. One chapter is about USDFL. Hope it will be helpful for you.

Regards,
Bartosz

 

[smallpea14]

Hi Bartosz,

I really appreciate for your explanation and nice material.
Yes, please help forward this to Simulia help support if applicable. I have no access to that.
Also, thanks for your explanation about calling sequence of USDFLD. I realized this when I compared the output of state value variable (say sp3 defined by myself)and found out that the state value in n increment was the same as the field variable, say sp3 in n-1 increment calculated by abaqus itself.

Ba

Based on your suggestion, I tried the same test with 3D block element.
I am really frustrated by these results. I used state variable output and fv output to check but didn't find anything strange except the Young's modulus I really care about. I run 3 tests, 3D inp without any usdfld, simple E defined as 21350.6 corresponding to 50 sp3; 3DUSRA inp with method a (E defined as a field variable) and 3DUSRB wih method b (SP3 defined as a field variable). I calculated the Young's modulus for each case by "smises" divided by "E33" for elastic part. None of results by both usdflds are correct. Each inp won't take you less than 2 mins for each run. Could you please take a look and give me more help? I thought this simple subroutine should work easily. I can't really attribute this to Abaqus bug although I still have no clue what mistake I may make. Anyway, again, thanks so much for your efforts, my friend,Bartosz.

Best Regards,
Muller.

 

[akabarten]

Hi,

I run your model with solid element. I am wondering one thing.

I need to define a E law dependent on SP3(minor principal stress)

It means S.Min.Principal output in odb should be the same as your statev(3), but they are not.

In documentation 2.1.6 Obtaining material point information in an Abaqus/Standard analysis, Abaqus User Subroutines Reference Manual I found:
"Three values are always returned for principal value requests, the minimum value first and maximum value third, regardless of the dimensionality of the analysis."
It means you minor principal stress are first and you should use array(1) instead array(3) in your equation.

What do you think about it? I did not run any test with that modification.

Regards,
Bartosz

 

[smallpea14]

Hi Bartosz,

Thanks so much for your efforts.

Yes, sp3 is maximum pricipal stress in Abaqus which defines tension is positive by default. But in my case, compression is positive and I want to define E dependent on -sp3 (minimal absolute value). So I don't think there is problem at this matter. Also, I calculate E based on the -sp3 (50kPa) and assigned the corresponding value to E as 21350.6 without usdfld.

I thought this might be caused by the initial field value or time increment for cases with usdfld but I don't know why and how to fix it.

For 3dusrA (method a), it gets in plastic state right away in first geostatic step. It diverged when I first assigned E (field variable here)between 50 to 1E6. I increased it to 100 in the inp file I attached last posting and let the program finish running although I got the unexpected results. As for 3dusrB (method b), although it didn't go to plastic state rightaway but in elastic part, the yong's modulus is not the expected value corresponding to filed variable (50kpa). Can you please give me more comments?

Best Regards,
Muller

 

[smallpea14]

Hi Bartosz,

If you come back to the forum, could you please help me to check my files again? I am really in the need of your help.
To be simple, I deleted my Mohr-Coloumb part in input files. For the three cases as described before, case A(3Dela) without any subroutine; Case B(3DUSRAela) with method a (E as field variable); Case C(3DUSRBela) with method b(-SP3 as field variable), they are exactly the same except the definition of young's modulus. The Smises were produced very differently by three cases even though their young's moduli were defined the same in my opinion.

I noticed that there were large displacement and smises during the geostatic step when using subroutines. I am really puzzled why this could happen considering that all the conditions were the same.

I attached these simplified, purely elastic files. Hopefully, you would have time to take a look, Bartosz. Thanks so much for your help in advance.

Best Regards,

Muller