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

4 weeks ago by
Hello everyone,

I'm modelling the deformation of a 3-layer Golf-ball using "". 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 ( but it didn't work.

Thanks in advance!
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.")
# Read mesh
# 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,\
           metadata={"quadrature_degree": 2} )
dV=Measure("dx", domain=mesh, subdomain_data=cells,\
           metadata={"quadrature_degree": 2} )

# 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 = I + grad(u)             # Deformation gradient
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
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)))
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)))
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
# 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):

    # 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

# 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():
elif MPI.size(mesh.mpi_comm()) == 1:
    encoding = XDMFFile.Encoding_ASCII
    out.write(u, encoding)
    # 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

1 Answer

4 weeks ago by
You are trying to extract the part of the subdomain ‘cells’ which is marked with number 1, but you did not mark any part of the subdomain. Try with:
Core = SubMesh(mesh, cells)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »