Incremental Stress Projection

dolfinx
1

Fenics Community,

I am trying to replicate an Abaqus User Material in dolfinx (docker v0.5.1-r1) that relates the increment of stress to the increment in strain, i.e.

\sigma{} = \sigma{}_{n} + \frac{\Delta{}\sigma{}}{\Delta{}\varepsilon{}}\delta{}\varepsilon{}

as the solution to

\nabla{}\cdot{}\sigma{}=0
\sigma{}\cdot{}n=0 on \Gamma{}
u=\bar{u} on \Gamma{}^{1}

The problem I am trying to solve involves many time dependent terms but so far I am unable to solve the above set of equations even using only elastic plane strain. The mechanical scenario is a three point flexural bending test subject to displacement controlled loading. I’ve simplified the stiffness tensor \frac{\Delta{}\sigma{}}{\Delta{}\varepsilon{}} to the plane strain Hooke’s law tensor for 2D. Additionally I define \delta{}\varepsilon{} as:

\delta{}\varepsilon{}(\delta{}u)=\frac{1}{2}(\nabla{}\delta{}u+(\nabla{}\delta{}u)^{T}) where \delta{}u=u-u_{n}, u_{n} is the previous displacement and u is the unknown displacement

I believe my issue is that I am not properly assigning \sigma{}_{n} at the end of each timestep as it does not match the solution given when \sigma{} = \frac{\Delta{}\sigma{}}{\Delta{}\varepsilon{}}\varepsilon{}(u) (perhaps something related to the change in function spaces during project()). In the application scenario \sigma{}_{n} evolves nonlinearly due anelastic/viscous material behavior so I can’t fully define \sigma{}_{n} only in terms of the previous displacement u_{n} (it would effectively require knowledge of \sigma{}_{n-1}). Any help on how to better implement the above constraints or corrections to mistakes I’ve overlooked is much appreciated. My attempt at an MWE is below:

from mpi4py import MPI
from dolfinx import mesh, fem, nls, io, cpp, log
import numpy as np
import ufl
import petsc4py.PETSc as pets

# MPI #
commWorld = MPI.COMM_WORLD
thisCore = commWorld.Get_rank()

# Saving #
fullDir = "Outputs/"

# Mesh #
domain = mesh.create_rectangle(commWorld, [[-1,0],[1,1]], [100, 5], ghost_mode=mesh.GhostMode.shared_facet)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
f_to_c = domain.topology.connectivity(fdim, tdim)
num_cells = domain.topology.index_map(tdim).size_local + domain.topology.index_map(tdim).num_ghosts
num_facets = domain.topology.index_map(fdim).size_local + domain.topology.index_map(fdim).num_ghosts

# Files #
loadFileName = fullDir + "_LOADDISP.txt"

# Files #
if (thisCore == 0):
     loadFile = open(loadFileName,'w')

# Constants #
scalar1 = fem.Constant(domain, pets.ScalarType(1.0))

# Inputs #
_E = 9e9
_nu = 0.0
_mu = _E/(2*(1+_nu))
_lmbda = (_nu*_E)/((1+_nu)*(1-2*_nu))
_coeff = _E/((1+_nu)*(1-2*_nu))
_c11 = _coeff * (1 - _nu)
_c12 = _nu * _coeff
_c33 = 0.5 * (1 - 2*_nu) * _coeff

_Compliance = ufl.as_tensor([[_c11, _c12, 0],
                             [_c12, _c11, 0],
                             [0,   0, _c33]])

# Boundary Parameters
tol = 1e-8 # Tolerance for mesh tagging (not for computations) [m]
comprTest = 0
LoadingTol = 0.05 + 1e3 * comprTest
relSupportDistance = 0.1 + 1e3 * comprTest

# Sim Conditions #
deltaLoad = 1.67e-5 # Displacement Rate [m/s]
dt1 = 1 # Starting timestep [s]
totalTime = 7 # Total Sim Time

# FunctionSpaces #
Vs_CG = fem.VectorFunctionSpace(domain, ("CG", 1))
Vs_CG_3 = fem.VectorFunctionSpace(domain, ("CG", 1), dim = 3)

# Functions #
u_k = fem.Function(Vs_CG)
u_k_old = fem.Function(Vs_CG, name="Displacement")
v_k = ufl.TestFunction(Vs_CG)
sigma_old, sigma_old_old = fem.Function(Vs_CG_3), fem.Function(Vs_CG_3)

# Boundary Definitions #

## SubDomain
class SubDomain:
    ## SubDomain.__init__()
    # Default class constructor without args
    def __init__(self):
        return

    ## Unused
    def on_boundary(self):
        return False

    ## SubDomain.lor()
    # A shorthand for np.logical_or()
    # - a,b: logical np arrays
    def lor(self, a, b): # Numpy Logical Or
        return np.logical_or(a, b)

    ## SubDomain.land()
    # A shorthand for np.logical_and()
    # - a,b: logical np arrays
    def land(self, a, b): # Numpy logical And
        return np.logical_and(a, b)

    ## SubDomain.inside()
    # A placeholder function to be overridden by the class instance.
    # Use a combination of SubDomain.lor(), SubDomain.land(), and SubDomain.near() to determine marked facets
    # - x: The point coordinates, provided when called by other methods such as SubDomain.mark()
    def inside(self, x): # Place Holder
        return True

    ## Unused
    def insideBoundary(self, x):
        print(self.inside(x))
        print(self._boundary)
        print(f"Inside {self.inside(x).shape}")
        print(self._boundary.shape)
        return self.land(self.inside(x), self._boundary)

    ## SubDomain.near()
    # Same functionality as dolfin.near()
    # Return whether a given value is absolutely close to another point
    # - xIn: The point to check
    # - val: The reference point
    # - atol: The absolute tolerance to be considered "near"
    def near(self, xIn, val, atol): # Same functionality as dolfin.near()
        return np.isclose(xIn, val, rtol=0, atol=atol)

    def assign(self, matProp, value): # Return Numpy array of matching indices
        markedArray = mesh.locate_entities(matProp._domain, matProp._dim, self.inside)
        matProp.x.array[markedArray] = np.full(len(markedArray), value)
        return

    ## SubDomain.mark()
    # Given a facetFunction, mark appropriate facets with value as determined by SubDomain.inside()
    # - facetFunction: The cdfx.FacetFunction to mark
    # - value: The value to assign to all facets
    def mark(self, facetFunction, value):
        if (not(self.on_boundary())):
            markedArray = mesh.locate_entities(facetFunction._domain, facetFunction._fdim, self.inside)
        else:
            self._boundary = np.array(mesh.compute_boundary_facets(facetFunction._domain.topology))
            markedArray = mesh.locate_entities(facetFunction._domain, facetFunction._fdim, self.insideBoundary)

        facetFunction._indices.append(markedArray)
        facetFunction._markers[markedArray] = value
        return

def DirichletBCs(argList, fs, facetTags):
    fdim = fs.mesh.topology.dim - 1
    bcs = []
    for i in argList:
        if 'Dirichlet' in argList[i]:
            if isinstance(argList[i]['Dirichlet'], fem.Constant):
                value = argList[i]['Dirichlet']
            else:
                value = fem.Constant(fs.mesh, pets.ScalarType(argList[i]['Dirichlet']))

            if 'Direction' in argList[i]:
                facets = np.array(facetTags.indices[facetTags.values==argList[i]['Surface']])
                left_dofs = fem.locate_dofs_topological(fs.sub(argList[i]['Direction']), fdim, facets)
                bcs.append(fem.dirichletbc(value, left_dofs, fs.sub(argList[i]['Direction'])))
            else:
                facets = np.array(facetTags.indices[facetTags.values==argList[i]['Surface']])
                left_dofs = fem.locate_dofs_topological(fs, fdim, facets)
                bcs.append(fem.dirichletbc(value, left_dofs, fs))

    return bcs

## FacetFunction
class FacetFunction:
    ## FacetFunction.__init__()
    # Default constructor with args
    # - domain: The dolfinx.mesh entity
    # - fdim: The facet dimension for this function [int]
    def __init__(self, domain, fdim, numEntities):
        self._domain = domain
        self._fdim = fdim
        self._indices = []
        self._numEntities = numEntities
        # self._markers = []
        self._markers = np.zeros(self._numEntities, dtype=np.int32)

    ## FacetFunction.CreateMeshTag()
    # Call this after using SubDomain.mark() to generate a dolfinx.mesh.meshtags object
    def CreateMeshTag(self):
        return mesh.meshtags(self._domain, self._fdim, np.arange(self._numEntities, dtype=np.int32), self._markers)

class RigidConstraint(SubDomain): # Bending Supports
    def inside(self, x):
        arg1 = self.near(x[1], 0, tol)
        arg2 = self.near(x[0], -0.8, relSupportDistance)
        arg3  = self.near(x[0], 0.8, relSupportDistance)
        return  self.land(arg1, self.lor(arg2, arg3))

class LoadingArea(SubDomain): # Loading Area
    def inside(self, x):
        arg1 = self.near(x[1], 1, tol)
        arg2 = self.near(x[0], 0, LoadingTol)
        return  self.land(arg1, arg2)

# Facets #
## FacetFunction instance
# Hold the marked facets of topological dimension 2
facetFunction = FacetFunction(domain, fdim, num_facets)
RigidConstraint().mark(facetFunction, 1)
LoadingArea().mark(facetFunction, 2)
facetTags = facetFunction.CreateMeshTag()

# Setup dz and ds #

dz = ufl.Measure('dx', domain = domain, metadata={"quadrature_degree": 5})
ds = ufl.Measure("ds", domain = domain, subdomain_data = facetTags)

# Load Expression #
loadArea = commWorld.allreduce(fem.assemble_scalar(fem.form(scalar1*ds(2))), op = MPI.SUM)
displacementDriven = fem.Constant(domain, pets.ScalarType(0.0))

# Boundary Conditions #

boundary_conditions_k = {1: {'Dirichlet': 0.0, 'Surface' : 1, 'Direction' : 1},
                         2: {'Dirichlet': displacementDriven, 'Surface' : 2, 'Direction' : 1},
                         3: {'Dirichlet': 0.0, 'Surface' : 2, 'Direction' : 0},
                         }

bcs_k = DirichletBCs(boundary_conditions_k, Vs_CG, facetTags)

integrals_N_k = []

def epsilon(u):
    return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def voigt2stress(vec):
    return ufl.as_tensor([[vec[0], vec[2]],
                      [vec[2], vec[1]]])

def StrainVector(e):
    return ufl.as_vector([e[0,0], e[1,1], 2*e[0,1]])

def Sigma(u_k, u_k_old, _oldsigma):
    dU = u_k - u_k_old
    dE = StrainVector(epsilon(dU))
    return _oldsigma + _Compliance * dE

def SigmaVec(sigmaIn):
    return ufl.as_vector([sigmaIn[0,0],sigmaIn[1,1],sigmaIn[0,1]])

def project(function, space):
    p = ufl.TrialFunction(space)
    q = ufl.TestFunction(space)
    a = ufl.inner(p, q) * dz
    L = ufl.inner(function, q) * dz
    problem = fem.petsc.LinearProblem(a, L, bcs = [])
    return problem.solve()

# Problem Declaration #

# Displacement #
F_k = ufl.inner(voigt2stress(Sigma(u_k, u_k_old, sigma_old_old)), epsilon(v_k))*dz - sum(integrals_N_k)

# Problem Creation #
Prob_k = fem.petsc.NonlinearProblem(F_k, u_k, bcs=bcs_k)

## Displacement #
Solver_k = nls.petsc.NewtonSolver(commWorld, Prob_k)
Solver_k.convergence_criterion = "residual"
Solver_k.atol = 1e-8
Solver_k.max_it = 25
Solver_k.report = True


# Write Load File #
if (thisCore == 0):
    loadFile.write("Time [s]\tDisplacement [m]\tActualLoad [N]\tStress [Pa]\tDelta Stress [Pa]\tError Value\n")
    loadFile.close()

## Loop ##
t = 0
dt = dt1
deltaLoad0 = 0
log.set_log_level(log.LogLevel.INFO)
while (t < totalTime):
    displacementDriven.value = deltaLoad0
    # Solve for k
    Solver_k.solve(u_k)

    # Solve for the current sigma using the current and previous displacement
    sigma_old.x.array[:] = project(Sigma(u_k, u_k_old, sigma_old_old), Vs_CG_3).x.array
    sigma_old.x.scatter_forward()

    # Assignments
    sigmaExact = _Compliance * StrainVector(epsilon(u_k))
    deltaValue = commWorld.allreduce(fem.assemble_scalar(fem.form((sigma_old - sigma_old_old)[1]*ds(2)))/loadArea, op = MPI.SUM)
    errorValue = commWorld.allreduce(fem.assemble_scalar(fem.form(ufl.inner(sigma_old - sigmaExact, sigma_old - sigmaExact)*ds(2)))/loadArea, op = MPI.SUM)
    sigma_old_old.x.array[:] = sigma_old.x.array
    u_k_old.x.array[:] = u_k.x.array
    # Saving
    dispValue = commWorld.allreduce(fem.assemble_scalar(fem.form(u_k[1]*ds(2)))/loadArea, op = MPI.SUM)
    loadValue = commWorld.allreduce(fem.assemble_scalar(fem.form(sigma_old_old[1]*ds(2))), op = MPI.SUM)
    stressValue = commWorld.allreduce(fem.assemble_scalar(fem.form(sigma_old_old[1]*ds(2)))/loadArea, op = MPI.SUM)
    if (thisCore == 0):
        loadFile = open(loadFileName,'a')
        loadFile.write("{:.4f}".format(t))
        loadFile.write("\t{:1.3E}".format(dispValue))
        loadFile.write("\t{:1.3E}".format(loadValue))
        loadFile.write("\t{:1.3E}".format(stressValue))
        loadFile.write("\t{:1.3E}".format(deltaValue))
        loadFile.write("\t{:1.3E}".format(errorValue))
        loadFile.write("\n")
        loadFile.close()

    deltaLoad0 -= deltaLoad*dt
    t += dt

The comprTest variable simply turns the loading scenario into a compression test if = 1 and flexural loading if = 0.

2

Your problem and code is too long for me to parse, but perhaps @bleyerj 's excellent tutorial will help Welcome to Numerical tours of Computational Mechanics using FEniCS — Numerical tours of continuum mechanics using FEniCS master documentation.

2 Likes
3

@nate ,

Thank you for the response. I believe that I’ve overlooked the element family I chose for the previous stress state \sigma{}_{n}. Thanks to @bleyerj 's page, specifically Elasto-plastic analysis of a 2D von Mises material and the post Quadrature function space - I realized that for my problem declaring:

is the issue and instead:

Vs_CG_3 = fem.VectorFunctionSpace(domain, (“DG”, 0), dim = 3)

is sufficient given my choice of function space order for the displacement u. I also am using fem.function.interpolate() for going between two different function spaces rather than project() such as from Dolfinx Discontinuous Expressions SIGSEGV.

Thanks,

Sven

1 Like