Abaqus scripting - Lode Angle 1

[Manuel Costa]

Hello, guys!

I've been trying to write a script to express the Lode angle (Polar deviatoric angle), as defined here: Ductile damage criterion wit lode angle dependance (Abaqus Analysis User's Guide).

Shortly I'd like to express the variable defined in the link, ξ (in range of -1<ξ<1) in a viewport .odb, for every step, for every frame.

Do you have any idea on how to do that?

This is what a did so far:

_____________________________
from abaqusConstants import *
from odbAccess import *

odbPath = "/path/to/.odb" # path to output database

odb = session.openOdb(name=odbPath, readOnly=False)
allSteps = session.odbData[odbPath].steps.keys('')

for i in range(len(allSteps)):
step = odb.steps[allSteps]
allFrames = session.odbData[odbPath].steps[allSteps].frames.keys('')

for j in range(len(allFrames)):
frame = step.frames[j+1]
mises = frame.fieldOutputs['S'].getScalarField(MISES)
inv3 = frame.fieldOutputs['S'].getScalarField(PRESS)
ksi = somefunctionoflodeangle
newField = frame.FieldOutput(name='Lodetheta', description='Lode angle', field=somefunctionoflodeangle)
print 'stepName = ', allSteps, ' frameNumber = ', allFrames[j]

odb.save()
odb.close()

 

[FEA way]

Maybe this could be done without scripting. Try using Tools --> Create Field Output --> From Fields. There you can create new field output based on available variables. You will need MISES and INV3.

 

[Shaun E.]

Let me know if this works for you. I had to use a variable "eps" to correct for divide by 0 and acos() range errors.

import math
import numpy as np
import odbAccess
import os
import shutil
import time
import datetime
from abaqus import *
from abaqusConstants import *
from numpy.lib.recfunctions import append_fields
from operator import itemgetter

################################################################################
# Functions for ABAQUS material modeling - Shaun E. 2020
################################################################################

def getODB(ODBNum=None):
    '''
    Gets ODB object by number in current session. Does not need a number if only odb is open.
    '''
    if ODBNum is None:
        odbs = session.odbs.items()
        if len(odbs) == 1:
            odb = session.odbs.items()[0][1]
        else:
            print "Only 1 ODB file can be open for getODB() to work without an integer input"
            odb = None
    else:
        odb = session.odbs.items()[ODBNum][1]
    return odb

def addLodeField(myODB, elementType=None, instanceName=None):
	'''
	Calculate normalized third invariant and normalized lode angle
	e.g. addLodeField(myODB,elementType='C3D8R',region='SOLID-PATCH-INNER-1')
	'''
	steps = myODB.steps
	ii = steps.keys()
	for i in ii:
	    myStep = steps[i]
	    # Iterate over frames for given step
	    frames = myStep.frames
	    n_f = len(frames)
	    for j in xrange(n_f):
	        frame=frames[j]
	        t_f_start = time.time()
	        if 'NormInv3' in frame.fieldOutputs.keys():
	            continue
	        if 'Lode' in frame.fieldOutputs.keys():
	            continue
	        S_all = frame.fieldOutputs['S'] # Stress Field Output Object 
	        # Get subsets of field output object if specified
	        if instanceName is not None:
	        	instance = myODB.rootAssembly.instances[instanceName]
	        	S_all = S_all.getSubset(INTEGRATION_POINT, region=instance)
	        if elementType is not None:
	        	S_all = S_all.getSubset(INTEGRATION_POINT, elementType=elementType)
	        mises = S_all.getScalarField(invariant=MISES)
	        inv3 = S_all.getScalarField(invariant=INV3)
	        # eps ensures (mises + eps) > abs(inv3) for acos calculation 
	        # eps prevents division by 0
	        # eps determined empirically with units of psi 
	        eps = 0.1 
	        normInv3 = power(inv3/(mises + eps), 3) 
	        lode = 1.-2./math.pi*acos(normInv3)
	        normInv3Field = frame.FieldOutput(name='NormInv3', 
	            description='Normalized Third Invariant of Deviatoric Stress Tensor', 
	            field=normInv3)
	        lodeField = frame.FieldOutput(name='Lode', 
	            description='Normalized Lode Angle', 
	            field=lode)
	        print('Processed frame '+ str(j) + ' / ' + str(n_f - 1) + ' in ' + 
	            myStep.name + ' step.')


################################################################################
# End of definitions
################################################################################

# Backup ODB, close, and open backup with write permissions
myODB = getODB()
path1 = myODB.path
path2 = path1[:-4] + '_lode.odb'
shutil.copy(path1, path2)
myODB.close()
myODB = odbAccess.openOdb(path=path2)

# Add New Field Outputs for normalized Lode angle and normalized third invariant
addLodeField(myODB, elementType=None, instanceName=None)

# Save, close, and reopen ODB to get fields to display properly
myODB.save()
myODB.close()
myODB = odbAccess.openOdb(path=path2, readOnly=True)

# Change viewport to show normalized Lode angle
session.viewports['Viewport: 1'].setValues(displayedObject=myODB)
session.viewports['Viewport: 1'].odbDisplay.display.setValues(plotState=(CONTOURS_ON_DEF, ))
session.viewports['Viewport: 1'].odbDisplay.setFrame(step=0, frame=0 )
session.viewports['Viewport: 1'].odbDisplay.contourOptions.setValues(
    maxAutoCompute=OFF, minAutoCompute=OFF)
session.viewports['Viewport: 1'].odbDisplay.contourOptions.setValues(maxValue=1.0001, 
    minValue=-1.0001)
session.viewports['Viewport: 1'].view.setValues(session.views['Iso'])
session.viewports['Viewport: 1'].odbDisplay.setPrimaryVariable(variableLabel='Lode',
    outputPosition=INTEGRATION_POINT)
session.viewports['Viewport: 1'].view.setProjection(projection=PARALLEL)