Modeling stress dependency of stiffness with USDFLD

[pevoz23]

Hi,

I'm trying to model a stress dependency of the Young's modulus of an elasto-plastic soil (Mohr-Coulomb) in Abaqus using the subroutine USDFLD (in plane strain conditions with CPE4 mesh elements).

My first simple implementation is:

      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,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(*),COORD(*)
      DOUBLE PRECISION PE11, PE22, PE33, PE12, PE23, PE31

      call GETVRM('S',ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,
     1 LACCFLA)

	  FIELD(1) = 1.0D0/3.0D0*(ARRAY(1)+ARRAY(2)+ARRAY(3))
	  

      RETURN
      END

This implementation however leads to convergence problems as soon as it is used (Geostatic step if I use it on this step, on the first static increment otherwise).

Can anyone explain me why and how can I improve this first implementation?
In this old post, Bartosz said:

The subroutine works ok for reduced elements (with one integration points) if you want to use fully integrated elements then you need to modify the subroutine.
I believe the best solution is to calculate E value for all material points in an element and assign one average value.

Is this my problem? How can I achieve this?

Thanks in advance

EDIT: THe problem should be the fact that at the beginning of my geostatic step I have FIELD(1)=0 and the change in stiffness is too aggressive. But how can I initialize the field variable for t=0? The value of this variable namely changes for every node of my soil part.

 

[akabarten]

Hi,

I do not think the problem is related to full integrated element type you are using.
Each integration point is treated separately in you example so it should not be a problem.

But how can I initialize the field variable for t=0? The value of this variable namely changes for every node of my soil part.

Is there any rule for initial distribution which can be describe with an equation?

One more point, are you sure that the subroutine works well.
Did you check FV field output to see values of E in your model?
Please also check JRCD flag from GETVRM routine to be sure the results have been obtained with out any errors.

Regards,
Bartosz

VIM filetype plugin for Abaqus

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

 

[pevoz23]

Hi Bartosz!
Thank you for your reply!
Yes, the field variable is the volumetric stress and since I know the initial stress distribution I could calculate the values for each point using its coordinates. The problem is that I don't know how to initialize a field variable which is different for each node of my set (yes, I probably could initialize it in the input file with a line for each node, but it is not a very good implementation, also considering that I am generating and meshing the model using the abaqus python interface).
I checked the output for FV and I'm now pretty sure that the problem is that at Step 1, increment 0 I have FV=0. Indeed, if the Young's modulus of my material is defined over the range field=[-1,-50] and I set with USDFLD a constant FIELD(1)=-1.1 for all points and increments, the problem converges. If I set, e.g., a constant FIELD(1)<-1.6 I have no convergence. However, if I define the Young's modulus only over the range field=[-6,-50] and FIELD(1)=-6.1 I get a solution (because the change in stiffness due to the change of FIELD(1) from 0 to -6.1 is small - due to the fact that abaqus does not extrapolate values outside the definition range).

Now I tried to calculate the model once with only the geostatic step activated and without table definition for the Young's modulus (so that the stiffness does not change). The idea is then to save the FV in the .fil file (this is probably still not working; I think I have to use *NODE OUTPUT \n FV, but I still get errors) and to load these value in the true model with *INITIAL CONDITION, FILE=..., STEP=1, INC=1. What do you think? Is there a better solution?

Any idea is welcome. Thank you so much.

EDIT: the problem is that if I put in the first input file the following string:

*NODE FILE,FREQ=1
FV,

I then get the following error message in the .dat file:

 ***ERROR: UNIDENTIFIED OR INVALID OUTPUT REQUEST FV
 ***NOTE: DUE TO AN INPUT ERROR THE ANALYSIS PRE-PROCESSOR HAS BEEN UNABLE TO 
          INTERPRET SOME DATA.  SUBSEQUENT ERRORS MAY BE CAUSED BY THIS OMISSION

With *EL FILE it works, but then I get

***ERROR: INITIAL CONDITIONS FOR FIELD VARIABLE 1 WERE NOT FOUND FOR STEP 1 
           INCREMENT 1 OF FILE field_job
 LINE IMAGE: *initialconditions, type=FIELD, file=field_job.fil, step=1, inc=1
 ***NOTE: DUE TO AN INPUT ERROR THE ANALYSIS PRE-PROCESSOR HAS BEEN UNABLE TO 
          INTERPRET SOME DATA.  SUBSEQUENT ERRORS MAY BE CAUSED BY THIS OMISSION

Are the field from USDFLD not extrapolated to the nodes? How can I solve this problem?

 

[akabarten]

Hi,

I have no experience with FV initialization but in documentation is information that values are read from output database so I understand from odb file.
Try to save FV outputs to odb (*NODE OUTPUT) and than use *INITIAL CONDITION, FILE=...

I have two ideas:

1. USDLD and 2 steps
I would create first step to just initialize FV in the model with USDFLD and second step with my true calculation.
First step should have no loads and maybe even boundary condition for all nodes. In next step you can remove boundary conditions.
With such model we should not expect convergence problems in 1st step.
The advantage is that Abaqus in automatic way carry over FV values to from 1st to 2nd step.

! initialize FV in 1st step
if (KSTEP == 1) then

  ... initialize FV/E --> base on COORD array ...

! set FV for other steps
else

  ... set FV/E --> base on stresses from GETVRM routine ...

endif

2. USDLFD & UFIELD
Never used such combination but I found this:

Before user subroutine USDFLD is called, the values of the field variables at the material point are calculated by interpolation from the values defined at the nodes. Any changes to the field variables in the user subroutine are local to the material point: the nodal field variables retain the values defined as initial conditions, predefined field variables, or in user subroutine UFIELD.

Perhaps we can use UFIELD for initialization of FV in USDFLD in 1st increment.

Regards,
Bartosz

VIM filetype plugin for Abaqus

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

 

[pevoz23]

Thanks for the reply.
I discovered that *INITIAL CONDITION with odb is only possible for

nodal values of temperature (NT), normalized concentrations (NNC), and electric potential (EPOT)

.

I also had approximately your ideas. Today I tried to implement the first one: I added a step after my geostatic step and i controlled via subroutine the increment time, so that the change in field variable wasn't too aggressive. The problem is that I get deformations and changes in stresses. Therefore I should really try as you said before applying my forces. The problem is that as far as I know, the geostatic stress can only be initialized in the Initial step. After this step (I think) I therefore have calculate a geostatic step, and I can not insert a step between them (I also tried it, with BC for all nodes, but didn't lead to the desired result). Have you ideas for a better implementation?
On wednesday, I also quickly looked to the UFIELD subroutine, but I think I gave up too early because I hadn't time. Perhaps I should look at it again.

My last idea was to edit the keywords with python scripting (it should be possible, isn't it?) and add an initial condition for each horizontal row of nodes (by looping over all mesh nodes).

What do you think about?

Regards,
David

 

[akabarten]

Hi David,

The problem is that as far as I know, the geostatic stress can only be initialized in the Initial step

Initial step always initialize only first step, doesn't matter what type it is (static, dynamic, ...).
But all results are passed to another step. So if you have first step which doesn't change your initial values you will get them also in second step.

edit the keywords with python scripting

Sounds good, you can use python to save external file and just use *INCLUDE to add it to your model.

Regards,
Bartosz

VIM filetype plugin for Abaqus

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

 

[Moustafa Elsawy]

Hi, I am sorry to bring this post back to life. David if you see this can you post what finally solved this issue for you? I am stuck on the same issue.

Thank you in davance!

 

[pevoz23]

Hi Moustafa. Sorry I never posted my final solution. I intended to do it, but never had time in the last months...
In any case, I solved the problem by initializing the field value for each node in my pressure dependent material with python.
To do this, I had to add them to the keywordBlock object (check out in the documentation how to use it - there is a method called insert -) or try to register a short macro.
If you need further help, I'm here.

 

[Moustafa Elsawy]

Hi Pevoz,

So sorry I havent seen your response until now. I have tried to read the documentation section for the insert function but still dont understand how you are initialising the FV for each node? Are you using the UFIELD subroutine? If you dont mind if you can explain in some detail as I kind of lost .

Moustafa

 

[Moustafa Elsawy]

Or is it using the *Initial Conditions in the input file and applying an initial condition to all nodes? If so, is there an easy way to get all nodes on the same level using a function or something?

Moustafa

 

[pevoz23]

Hi Moustafa.

Here my code (it should be self explaining):

[...]

# Keywords without nodes/elements
m.keywordBlock.synchVersions(storeNodesAndElements=False)
# Get position of *INITIAL blocks
init_blocks = [block for block in m.keywordBlock.sieBlocks if ('*INITIAL' in block or '*Initial' in block)]
# Her is where I want to add my keyword
id_init = m.keywordBlock.sieBlocks.index(init_blocks[-1])
# Create string
strfield = '*INITIAL CONDITIONS, TYPE=FIELD, VARIABLE=1\n'
# Concatenate value for each node
for partobj in cl.Part:
	if partobj.Material == material_obj.Name:
		surfload = 0
		for load in cl.Load:
			if load.type == 'Pressure' and load.part.Name == partobj.Name:
				surfload += load.magnitude
		for node in a.instances[partobj.Name + df.inststr].nodes:
			field_value = -float(1) / float(3) * (df.g * material_obj.Density * (
			partobj.yDim - partobj.yStart - node.coordinates[1]) + surfload) * (1 + 2 * material_obj.latcoeff)
			strfield += partobj.Name + df.inststr + '.' + str(node.label) + ', ' + str(field_value) + '\n'

# Insert keyword
m.keywordBlock.insert(id_init, strfield)

[...]

Hope it can help!

Regards

 

[Moustafa Elsawy]

Hi Pevoz,

thanks your fast reply. I am not an expert yet when it comes to controlling ABAQUS using Python as I always do an example with the macro recorder and work on from there. Would you mind just explaining a few points for me:

1) You are running this script after creating the model right? Like what I mean this is not part of the script used to create the model right?

2) What is cl, df? I am still not very good with the Access lines that each function uses so I need the full code line or path

3) A general idea on how you implement this code? (Running it after you create your input file or when exactly?) I do understand the logic behind how it works though.

Sorry my question might sound a little bit on the beginner side, I am still learning .

Moustafa

 

[pevoz23]

Hi Moustafa

Sorry, I just pasted my code which actually contains things you can't know, like df and cl (which are aliases of other modules I wrote).

I first reply your questions in the right order:

1) No, this is part of the script I use to create the model. Of course, however, I firstly create the parts, materials etc. and in the end I initialize the field variables. If you want, you can also have a separate script for initializing the fields... it's the same thing.

2) Sorry, as I said df and cl are aliases of the modules "definition" and "classes". I will give you a better MWE in the following.

3) Ok, now I understand what you really meant in question 1. I create the entire model using python. If you want, however, you can create your model directly in CAE or with an input file and then initialize your field variables.

Here a MWE.

Let's assume you have a model containing your soil part (or whatever else) named MOUSTAFA_PART. In the assembly you have the instance MOUSTAFA_PART_instance. This model was already created in CAE (or imported via input file or created via python). Now you want to initialize a field variable for each node dependent on the coordinates of every single node.

Your python script will then be something like:

# Assign your model to the variable m
m = mdb.models['your model']

# Assign your assembly to the variable a
m = m.rootAssembly

# Keywords without nodes/elements (Optional, but processing of sieblocks will be faster if storeNodesAndElements=False)
m.keywordBlock.synchVersions(storeNodesAndElements=False)

# You have now to find the position where you want to add your initial value for the field.
# You can, for example, add them as the last thing before the steps definitions.
# You therefore want to get the index of the *STEP keyword by iterating over all blocks and checking if it contains *STEP
step_blocks = [block for block in m.keywordBlock.sieBlocks if ('*STEP' in block or '*Step' in block or '*step' in block)]

# Because you probably have more than one step, you now want to insert your keyword before the first step
# Because the insert method inserts keywords AFTER the given index, subtract 1
id_step = m.keywordBlock.sieBlocks.index(step_blocks[0])-1

# For initializing fields you need the following keyword (check abq keyword reference in the docs)
# Followed by the value for each node
strfield = '*INITIAL CONDITIONS, TYPE=FIELD, VARIABLE=1\n'

# Loop over every node in your instance and calculate the field variable according to your needs
for node in a.instances['MOUSTAFA_PART_instance'].nodes:
        field_value = -float(1) / float(3) * materialdensity * (y0 - node.coordinates[1]) * (1 + 2 * K)
	strfield += 'MOUSTAFA_PART_instance.' + str(node.label) + ', ' + str(field_value) + '\n'

# Where y0 is some reference value for the position, K is the lateral coefficient of the earth pressure and the field value 
# is assumed to be equal the volumetric stress in your soil part.

# Insert keywordblock
m.keywordBlock.insert(id_step, strfield)

Now it should be clearer.

Best regards.

 

[Moustafa Elsawy]

Thank you so much Pevoz!!

I have edited it to my needs and it looks like its working as I am seeing the change in the input file. Now whats left is to run the model and see if the results make sense. Will keep you updated if something goes wrong. If you dont hear from me means im good to go!

Thank you again