Initial commit

#########################################################################

#

# To run in Abaqus:

# File -> Run Script

#

#########################################################################

#########################################################################

#

#

# SCRIPT FOR SMA BENDING TUBE SIMULATION (BY CASMART)

#

# This script generates a model of an SMA bending tube.

# For details in the model experimental setup see:

# Reedlunn et al, JMPS, 2013 http://dx.doi.org/10.1016/j.jmps.2012.12.012

#

#

# The model uses symmetry. 1/4 of the tube is modeled.

# Abaqus implemantation of the Auricchio/Taylor Model is utilized

#

#

# Presented at: SMASIS 2013

# P. Zhu, E. Peraza-Hernandez, D. Hartl, A. Stebner, C. Brinson,

# Comparison of three-dimensional shape memory alloy constitutive models:

# finite element analysis of actuation and superelastic responses of a

# shape memory alloy tube

#

#

# FEA model created by:

# Edwin Peraza-Hernandez, and Darren Hartl (darren.hartl@tamu.edu)

#

#

# Department of Aerospace Engineering

# Texas A&M University

# 3409 TAMU

# College Station, TX 77843-3409

#

#

#########################################################################

#########################################################################

# INPUT PARAMETERS

# Current units: m, kg, s

# Geometric parameters

inner_radius=0.00127 # tube inner radius

outer_radius=0.00159 # tube outer radius

tube_half_length=0.0405 # half length of the tube

bearing_radius=0.0024 # radius of the bearings

LE=0.00958 # examined length

bearing_dist=0.01577 # distance between the bearings

# Conditions

Room_Temperature=345.0 # room temperature

rotation_bearing=0.5 # bearing rotation [radians]

# Meshing

meshsize = 0.00075 # Approximate mesh size

# Material properties

SMA_Density = 6450.0

SMA_E_Austenite = 70.0e9 # Elastic modulus austenite

SMA_E_Martensite = 70.0e9 # Elastic modulus martensite

SMA_v_Austenite = 0.33 # Poisson's ratio austenite

SMA_v_Martensite = 0.33 # Poisson's ratio martensite

SMA_H = 0.047 # Transformation strain

SMA_C_Martensite = 6.31e6 # Stress influence coefficient martensite

SMA_Sigma_Ms = 448.0e6 # Initiation stress for transformation into martensite

SMA_Sigma_Mf = 562.0e6 # Completion stress for transformation into martensite

SMA_T0 = 350.0 # Reference temperature

SMA_C_Austenite = 9.21e6 # Stress influence coefficient austenite

SMA_Sigma_As = 257.0e6 # Initiation stress for transformation into austenite

SMA_Sigma_Af = 221.0e6 # Initiation stress for transformation into austenite

SMA_Sigma_Ms_comp = 448.0e6 # Initiation stress for transformation into martensite (Compression)

SMA_Hv = 0.047 # Volumetric transformation strain

#########################################################################

#########################################################################

# OUTPUTS

# In the Field Output section, select State Dependent Variables SDV to obtain:

# SDV 1 to 6 Linear elastic strains

# SDV 7 to 12 Transformation strains

# SDV19 Equivalent transformation strain

# SDV21 Martensite volume fraction

#########################################################################

#########################################################################

# -*- coding: mbcs -*-

# Do not delete the following import lines

from abaqus import *

from abaqusConstants import *

import __main__

import os

import section

import regionToolset

import displayGroupMdbToolset as dgm

import part

import material

import assembly

import step

import interaction

import load

import mesh

import optimization

import job

import sketch

import visualization

import xyPlot

import displayGroupOdbToolset as dgo

import connectorBehavior

a='/home'

b='apadhy'

c='tryout'

d='PhD'

e='ABAQUS'

f='UMAT_implementation'

g='Elastic/'

inp_dir= os.path.join(a,b,c,d,e,f,g)

os.chdir(inp_dir)

sub= os.path.join(a,b,c,d,e,f,g,'Elasticity.f')

# CREATING TUBE

s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=1.0)

g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints

s1.setPrimaryObject(option=STANDALONE)

s1.ConstructionLine(point1=(0.0, -0.5), point2=(0.0, 0.5))

s1.FixedConstraint(entity=g[2])

s1.rectangle(point1=(inner_radius, 0.0), point2=(outer_radius, tube_half_length))

p = mdb.models['Model-1'].Part(name='Tube', dimensionality=THREE_D,

    type=DEFORMABLE_BODY)

p = mdb.models['Model-1'].parts['Tube']

p.BaseSolidRevolve(sketch=s1, angle=180.0, flipRevolveDirection=OFF)

s1.unsetPrimaryObject()

p = mdb.models['Model-1'].parts['Tube']

del mdb.models['Model-1'].sketches['__profile__']

# CREATING BEARING

s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=1.0)

g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints

s1.setPrimaryObject(option=STANDALONE)

s1.ArcByCenterEnds(center=(0.0, 0.0), point1=(0.065, 0.0), point2=(-0.08,

    0.01), direction=COUNTERCLOCKWISE)

s1.FixedConstraint(entity=v[2])

s1.RadialDimension(curve=g[2], textPoint=(0.1,

    0.1), radius=bearing_radius)

p = mdb.models['Model-1'].Part(name='Bearing', dimensionality=THREE_D,

    type=ANALYTIC_RIGID_SURFACE)

p = mdb.models['Model-1'].parts['Bearing']

p.AnalyticRigidSurfExtrude(sketch=s1, depth=0.01)

s1.unsetPrimaryObject()

p = mdb.models['Model-1'].parts['Bearing']

del mdb.models['Model-1'].sketches['__profile__']

p = mdb.models['Model-1'].parts['Bearing']

p.ReferencePoint(point=(-bearing_dist/2.0, bearing_radius+outer_radius, 0.0))

# PARTITIONING TUBE

p = mdb.models['Model-1'].parts['Tube']

p.DatumPlaneByPrincipalPlane(principalPlane=XZPLANE, offset=LE/2.0)

p = mdb.models['Model-1'].parts['Tube']

c = p.cells

pickedCells = c.getSequenceFromMask(mask=('[#1 ]', ), )

d1 = p.datums

p.PartitionCellByDatumPlane(datumPlane=d1[2], cells=pickedCells)

# TUBE SETS

p = mdb.models['Model-1'].parts['Tube']

c = p.cells

cells = c.getSequenceFromMask(mask=('[#3 ]', ), )

p.Set(cells=cells, name='WholeTube')

p = mdb.models['Model-1'].parts['Tube']

s = p.faces

side1Faces = s.getSequenceFromMask(mask=('[#90 ]', ), )

p.Surface(side1Faces=side1Faces, name='OuterTubeSurface')

# CREATING ASSEMBLY

a = mdb.models['Model-1'].rootAssembly

a.DatumCsysByDefault(CARTESIAN)

p = mdb.models['Model-1'].parts['Tube']

a.Instance(name='Tube-1', part=p, dependent=ON)

a = mdb.models['Model-1'].rootAssembly

p = mdb.models['Model-1'].parts['Bearing']

a.Instance(name='Bearing-1', part=p, dependent=ON)

a = mdb.models['Model-1'].rootAssembly

p = mdb.models['Model-1'].parts['Bearing']

a.Instance(name='Bearing-2', part=p, dependent=ON)

# POSITIONING

a = mdb.models['Model-1'].rootAssembly

a.rotate(instanceList=('Bearing-2', ), axisPoint=(0.0, 0.0, 0.0),

    axisDirection=(0.0, 0.0, 1.0), angle=90.0)

a = mdb.models['Model-1'].rootAssembly

a.rotate(instanceList=('Bearing-1', ), axisPoint=(0.0, 0.0, 0.0),

    axisDirection=(0.0, 0.0, 1.0), angle=-90.0)

a = mdb.models['Model-1'].rootAssembly

a.translate(instanceList=('Bearing-2', ), vector=(outer_radius+bearing_radius, LE/2.0+bearing_dist, 0.0))

a = mdb.models['Model-1'].rootAssembly

a.translate(instanceList=('Bearing-1', ), vector=(-outer_radius-bearing_radius, LE/2.0, 0.0))

# BOUNDARY CONDITIONS

a = mdb.models['Model-1'].rootAssembly

region = a.instances['Tube-1'].sets['WholeTube']

mdb.models['Model-1'].Temperature(name='Initial_temperature',

    createStepName='Initial', region=region, distributionType=UNIFORM,

    crossSectionDistribution=CONSTANT_THROUGH_THICKNESS, magnitudes=(Room_Temperature,

    ))

a = mdb.models['Model-1'].rootAssembly

f1 = a.instances['Tube-1'].faces

faces1 = f1.getSequenceFromMask(mask=('[#100 ]', ), )

region = regionToolset.Region(faces=faces1)

mdb.models['Model-1'].YsymmBC(name='Symmetry_Y', createStepName='Initial',

    region=region, localCsys=None)

a = mdb.models['Model-1'].rootAssembly

f1 = a.instances['Tube-1'].faces

faces1 = f1.getSequenceFromMask(mask=('[#60a ]', ), )

region = regionToolset.Region(faces=faces1)

mdb.models['Model-1'].ZsymmBC(name='Symmetry_Z', createStepName='Initial',

    region=region, localCsys=None)

a = mdb.models['Model-1'].rootAssembly

r1 = a.instances['Bearing-2'].referencePoints

refPoints1=(r1[2], )

region = regionToolset.Region(referencePoints=refPoints1)

mdb.models['Model-1'].DisplacementBC(name='Bearing1', createStepName='Initial',

    region=region, u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET,

    amplitude=UNSET, distributionType=UNIFORM, fieldName='',

    localCsys=None)

a = mdb.models['Model-1'].rootAssembly

r1 = a.instances['Bearing-1'].referencePoints

refPoints1=(r1[2], )

region = regionToolset.Region(referencePoints=refPoints1)

mdb.models['Model-1'].DisplacementBC(name='Bearing2', createStepName='Initial',

    region=region, u1=SET, u2=SET, u3=SET, ur1=SET, ur2=SET, ur3=SET,

    amplitude=UNSET, distributionType=UNIFORM, fieldName='',

    localCsys=None)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Bearing-2'].faces

side2Faces1 = s1.getSequenceFromMask(mask=('[#1 ]', ), )

region5=regionToolset.Region(side2Faces=side2Faces1)

a = mdb.models['Model-1'].rootAssembly

r1 = a.instances['Bearing-2'].referencePoints

refPoints1=(r1[2], )

region1=regionToolset.Region(referencePoints=refPoints1)

mdb.models['Model-1'].RigidBody(name='Bearing1', refPointRegion=region1,

    surfaceRegion=region5)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Bearing-1'].faces

side2Faces1 = s1.getSequenceFromMask(mask=('[#1 ]', ), )

region5=regionToolset.Region(side2Faces=side2Faces1)

a = mdb.models['Model-1'].rootAssembly

r1 = a.instances['Bearing-1'].referencePoints

refPoints1=(r1[2], )

region1=regionToolset.Region(referencePoints=refPoints1)

mdb.models['Model-1'].RigidBody(name='Bearing2', refPointRegion=region1,

    surfaceRegion=region5)

# MATERIALS

#mdb.models['Model-1'].Material(name='Elastic_Test_Material')

#mdb.models['Model-1'].materials['Elastic_Test_Material'].Elastic(table=((

#    70000000000.0, 0.3), ))

#mdb.models['Model-1'].materials['Elastic_Test_Material'].Density(table=((

#    6450.0, ), ))

#

#mdb.models['Model-1'].HomogeneousSolidSection(name='Elastic_test_section',

#    material='Elastic_Test_Material', thickness=None)

# SMA:

mdb.models['Model-1'].Material(name='ELASTIC1')

mdb.models['Model-1'].materials['ELASTIC1'].Density(table=((SMA_Density, ), ))

mdb.models['Model-1'].materials['ELASTIC1'].Depvar(n=24)

mdb.models['Model-1'].materials['ELASTIC1'].UserMaterial(mechanicalConstants=(SMA_E_Austenite,SMA_v_Austenite))

    # SMA_E_Martensite, SMA_v_Martensite, SMA_H, SMA_C_Martensite,

    #    SMA_Sigma_Ms, SMA_Sigma_Mf, SMA_T0, SMA_C_Austenite,

    #    SMA_Sigma_As, SMA_Sigma_Af, SMA_Sigma_Ms_comp, SMA_Hv, 0.0))

mdb.models['Model-1'].HomogeneousSolidSection(name='SMA_AT_section',

    material='ELASTIC1', thickness=None)

p = mdb.models['Model-1'].parts['Tube']

c = p.cells

cells = c.getSequenceFromMask(mask=('[#3 ]', ), )

region = regionToolset.Region(cells=cells)

p = mdb.models['Model-1'].parts['Tube']

p.SectionAssignment(region=region, sectionName='SMA_AT_section',

    offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='',

    thicknessAssignment=FROM_SECTION)

p = mdb.models['Model-1'].parts['Bearing']

r = p.referencePoints

refPoints=(r[2], )

p.Set(referencePoints=refPoints, name='REF_POINT')

# CREATING INTERACTIONS

mdb.models['Model-1'].ContactProperty('IntProp-1')

mdb.models['Model-1'].interactionProperties['IntProp-1'].TangentialBehavior(

    formulation=FRICTIONLESS)

mdb.models['Model-1'].interactionProperties['IntProp-1'].NormalBehavior(

    pressureOverclosure=LINEAR, contactStiffness=10000000000000000.0,

    constraintEnforcementMethod=DEFAULT)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Bearing-2'].faces

side2Faces1 = s1.getSequenceFromMask(mask=('[#1 ]', ), )

region1=regionToolset.Region(side2Faces=side2Faces1)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Tube-1'].faces

side1Faces1 = s1.getSequenceFromMask(mask=('[#80 ]', ), )

region2=regionToolset.Region(side1Faces=side1Faces1)

mdb.models['Model-1'].SurfaceToSurfaceContactStd(name='Int-1',

    createStepName='Initial', master=region1, slave=region2,

    sliding=FINITE, thickness=ON, interactionProperty='IntProp-1',

    surfaceSmoothing=NONE, adjustMethod=TOLERANCE, initialClearance=OMIT,

    datumAxis=None, clearanceRegion=None, tied=OFF, adjustTolerance=1e-05)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Bearing-1'].faces

side2Faces1 = s1.getSequenceFromMask(mask=('[#1 ]', ), )

region1=regionToolset.Region(side2Faces=side2Faces1)

a = mdb.models['Model-1'].rootAssembly

s1 = a.instances['Tube-1'].faces

side1Faces1 = s1.getSequenceFromMask(mask=('[#90 ]', ), )

region2=regionToolset.Region(side1Faces=side1Faces1)

mdb.models['Model-1'].SurfaceToSurfaceContactStd(name='Int-2',

    createStepName='Initial', master=region1, slave=region2,

    sliding=FINITE, thickness=ON, interactionProperty='IntProp-1',

    adjustMethod=TOLERANCE, initialClearance=OMIT, datumAxis=None,

    clearanceRegion=None, tied=OFF, adjustTolerance=1e-05)

# CREATING STEP

mdb.models['Model-1'].ImplicitDynamicsStep(name='Loading', previous='Initial',

    maxNumInc=10000, initialInc=0.1, minInc=1e-06, nlgeom=ON)

mdb.models['Model-1'].ImplicitDynamicsStep(name='Unload', previous='Loading',

    maxNumInc=10000, initialInc=0.1, minInc=1e-06)

# AMPLITUDE

mdb.models['Model-1'].SmoothStepAmplitude(name='Amp-1', timeSpan=STEP, data=((

    0.0, 0.0), (1.0, 1.0)))

mdb.models['Model-1'].SmoothStepAmplitude(name='Amp-2', timeSpan=STEP, data=((

    0.0, 1.0), (1.0, 0.0)))

mdb.models['Model-1'].boundaryConditions['Bearing1'].setValuesInStep(

    stepName='Loading', ur3=rotation_bearing)

mdb.models['Model-1'].boundaryConditions['Bearing2'].setValuesInStep(

    stepName='Loading', ur3=rotation_bearing, amplitude='Amp-1')

mdb.models['Model-1'].boundaryConditions['Bearing1'].setValuesInStep(

    stepName='Loading', amplitude='Amp-1')

mdb.models['Model-1'].boundaryConditions['Bearing1'].setValuesInStep(

    stepName='Unload', amplitude='Amp-2')

mdb.models['Model-1'].boundaryConditions['Bearing2'].setValuesInStep(

    stepName='Unload', amplitude='Amp-2')

# MESHING

elemType1 = mesh.ElemType(elemCode=C3D20R, elemLibrary=STANDARD)

elemType2 = mesh.ElemType(elemCode=C3D15, elemLibrary=STANDARD)

elemType3 = mesh.ElemType(elemCode=C3D10, elemLibrary=STANDARD)

p = mdb.models['Model-1'].parts['Tube']

c = p.cells

cells = c.getSequenceFromMask(mask=('[#3 ]', ), )

pickedRegions =(cells, )

p.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2,

    elemType3))

p = mdb.models['Model-1'].parts['Tube']

p.seedPart(size=meshsize, deviationFactor=0.1, minSizeFactor=0.1)

p = mdb.models['Model-1'].parts['Tube']

p.generateMesh()

a = mdb.models['Model-1'].rootAssembly

a.regenerate()

# OUTPUT

mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(variables=(

    'S', 'PE', 'PEEQ', 'PEMAG', 'LE', 'U', 'V', 'A', 'RF', 'CF', 'CSTRESS',

    'CDISP', 'NT', 'COORD', 'SDV'))

# OUTPUT TIME FRAMES

mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(

    numIntervals=50, timeMarks=OFF)

# CREATING JOB

mdb.Job(name='Elastic', model='Model-1', description='bending tube job',

    type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None,

    memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True,

    explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF,

    modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine=sub,

    scratch='', multiprocessingMode=DEFAULT, numCpus=2, numDomains=2,

    numGPUs=0)

!DIR$ NOFREEFORM

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

     1 RPL,DDSDDT,DRPLDE,DRPLDT,

     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,

     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,

     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)

C

      INCLUDE 'ABA_PARAM.INC'

C

      CHARACTER*80 CMNAME

      DIMENSION STRESS(NTENS),STATEV(NSTATV),

     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),

     4 JSTEP(4)

C     ELASTIC USER SUBROUTINE

      PARAMETER (ONE=1.0D0, TWO=2.0D0)

	  E=PROPS(1)

	  ANU=PROPS(2)

      ALAMBDA=E/(ONE+ANU)/(ONE-TWO*ANU)

	  BLAMBDA=(ONE-ANU)

         CLAMBDA=(ONE-TWO*ANU)

	     DO I=1,NTENS

	      DO J=1,NTENS

	      DDSDDE(I,J)=0.0D0

	      ENDDO

	     ENDDO

		DDSDDE(1,1)=(ALAMBDA*BLAMBDA)

		DDSDDE(2,2)=(ALAMBDA*BLAMBDA)

		DDSDDE(3,3)=(ALAMBDA*BLAMBDA)

		DDSDDE(4,4)=(ALAMBDA*CLAMBDA)

		DDSDDE(5,5)=(ALAMBDA*CLAMBDA)

		DDSDDE(6,6)=(ALAMBDA*CLAMBDA)

		DDSDDE(1,2)=(ALAMBDA*ANU)

		DDSDDE(1,3)=(ALAMBDA*ANU)

		DDSDDE(2,3)=(ALAMBDA*ANU)

		DDSDDE(2,1)=(ALAMBDA*ANU)

		DDSDDE(3,1)=(ALAMBDA*ANU)

		DDSDDE(3,2)=(ALAMBDA*ANU)

         DO I=1,NTENS

	      DO J=1,NTENS

	      STRESS(I)=STRESS(I)+DDSDDE(I,J)*DSTRAN(J)

	      ENDDO

	     ENDDO

      RETURN

      END