### Project different functions of the solution in different subdomains (Von Mises Stress and Strain Energy in a 3-Materials Body)

80
views
1
4 weeks ago by
Hello everyone,

I'm modelling the deformation of a 3-layer Golf-ball using "demo_contact-vi-tao.py". With the displacement $u(x)$ (Solution) I want to calculate the Von-Mises Stress in each layer separately
$T_1 = f_1 (u) , \quad \forall x \in \Omega_1$
$T_2 = f_2 (u) , \quad \forall x \in \Omega_2$
$T_3 = f_3 (u) , \quad \forall x \in \Omega_3$
and somehow combine the solutions so I can view the Von-Mises Stress in all layers simultaneously in ParaView
$T = f_1 (u) + f_2 (u) + f_3 (u) , \quad \forall x \in \Omega$
$\Omega = \Omega_1 \cup \Omega_2 \cup \Omega_3$
I'm pretty new with FEniCS so I'm not really sure how to specify that. I tried the solution of a similar question (https://www.allanswered.com/post/ggmxb/how-to-evaluate-an-expression-only-on-a-subdomain/) but it didn't work.

from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt

if not has_petsc():
print("DOLFIN must be compiled with PETSc to run this demo.")
exit(0)
tic()
mesh=Mesh("./Geometrie/GB2DH.xml")
# Create function space
V=VectorFunctionSpace(mesh, "Lagrange", 1)
cells=MeshFunction("size_t", mesh, "./Geometrie/GB2DH_physical_region.xml")
facets=MeshFunction("size_t", mesh, "./Geometrie/GB2DH_facet_region.xml")
dA=Measure("ds", domain=mesh, subdomain_data=facets,\
dV=Measure("dx", domain=mesh, subdomain_data=cells,\

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, -1.0))     # Body force per unit volume

# Kinematics
I = Identity(len(u))  # Identity tensor
F = variable(F)
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
# Material constants
#Core
E0_core=50                     #Initial E-Module (MPa)
nu_core=0.45                  #Poisson ratio
rho0_core=1               #Density (ton/mm^3)
lmbda_core=Constant(E0_core*nu_core/((1. + nu_core)*(1. - 2.*nu_core)))
mu_core=Constant(E0_core/(2.*(1. + nu_core)))
#Mid
E0_mid=25                     #Initial E-Module (MPa)
nu_mid=0.45                    #Poisson ratio
rho0_mid=1               #Density (ton/mm^3)
lmbda_mid=Constant(E0_mid*nu_mid/((1. + nu_mid)*(1. - 2.*nu_mid)))
mu_mid=Constant(E0_mid/(2.*(1. + nu_mid)))
#Cover
E0_cover=400                   #Initial E-Module (MPa)
nu_cover=0.45                  #Poisson ratio
rho0_cover=1             #Density (ton/mm^3)
lmbda_cover=Constant(E0_cover*nu_cover/((1. + nu_cover)*(1. - 2.*nu_cover)))
mu_cover=Constant(E0_cover/(2.*(1. + nu_cover)))

# Stored strain energy density (compressible neo-Hookean model)
def psi(mu,lmbda):
return (mu/2)*(Ic-2) - mu*ln(J) + (lmbda/2)*(ln(J))**2

#Von Mises (First layer)
T_1PK = diff(psi(mu_cover,lmbda_cover), F)
T = 1.0/det(F)*dot(T_1PK, F.T)
T_dev = T-1.0/3.0*tr(T)*I
VM = sqrt(3.0/2.0*inner(T_dev,T_dev))

def P(mu, lmbda, rho0, G):
return psi(mu,lmbda)*dV(G)-rho0*dot(B,u)*dV(G)
# Total potential energy
Pi=P(mu_core,lmbda_core,rho0_core,1)+P(mu_mid,lmbda_mid,rho0_mid,6)\
+P(mu_cover,lmbda_cover,rho0_cover,2)
# Compute first variation of Pi (directional derivative about u in the
# direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Define the minimisation problem by using OptimisationProblem class
class ContactProblem(OptimisationProblem):

def __init__(self):
OptimisationProblem.__init__(self)

# Objective function
def f(self, x):
u.vector()[:] = x
return assemble(Pi)

# Gradient of the objective function
def F(self, b, x):
u.vector()[:] = x
assemble(F, tensor=b)

# Hessian of the objective function
def J(self, A, x):
u.vector()[:] = x
assemble(J, tensor=A)

# The displacement u must be such that the current configuration
# doesn't escape the box [xmin, xmax] x [ymin, ymax]
constraint_u = Expression(("xmax-x[0]", "ymax-x[1]"), xmax=30.0+DOLFIN_EPS, ymax=21.4+DOLFIN_EPS, degree=1)
constraint_l = Expression(("xmin-x[0]", "ymin-x[1]"), xmin=0.0-DOLFIN_EPS, ymin=-21.4-DOLFIN_EPS, degree=1)
u_min = interpolate(constraint_l, V)
u_max = interpolate(constraint_u, V)

# Symmetry condition (to block rigid body rotations)
#def symmetry_line(x, on_boundary):
#    return abs(x[0])<DOLFIN_EPS
bc = DirichletBC(V.sub(0), Constant(0.0), facets, 4)
bc.apply(u_min.vector())  # BC will be incorporated into the lower
bc.apply(u_max.vector())  # and upper bounds
# Create the PETScTAOSolver
solver = PETScTAOSolver()

# Set some parameters
solver.parameters["method"] = "tron"  # when using gpcg make sure that you have a constant Hessian
solver.parameters["linear_solver"] = "nash"
solver.parameters["line_search"] = "gpcg"
solver.parameters["preconditioner"] = "ml_amg"
solver.parameters["monitor_convergence"] = True
solver.parameters["report"] = True

# Uncomment this line to see the available parameters
# info(parameters, True)

# Parse (PETSc) parameters
parameters.parse()

# Solve the problem
solver.solve(ContactProblem(), u.vector(), u_min.vector(), u_max.vector())

# Save solution in XDMF format if available
out = XDMFFile(mesh.mpi_comm(), "./Ergebnis/GB2DHtao_defDISP.xdmf")
if has_hdf5():
out.write(u)
elif MPI.size(mesh.mpi_comm()) == 1:
encoding = XDMFFile.Encoding_ASCII
out.write(u, encoding)
else:
# Save solution in vtk format
out = File("./Ergebnis/DispGB2DHtao.pvd")
out << u

Core = SubMesh(mesh, cells, 1)
FS_core=FunctionSpace(Core, "Lagrange", 1)

VM_core = project(VM, FS_core)

File_VM = File("./Ergebnis/GB2DHtao_defVM3.pvd")
File_VM << VM_core
print ("Simulation time", toc())​
Community: FEniCS Project