Result visualization for uniaxial loading of a cube with two different maerials

1

HI,
I am trying to do uniaxial testing of a cube with two different nearly-incompressible NeoHookean materials for two sub-domains. My code is running fine. However, I do not know how to visualize the result. After the solver.solve() command, I calculate PK stress, but I can not figure out how to display both stresses in the mesh (As the mesh has two different materials, it should have two different stressed regions). Here is my code.

from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import math 
from mshr import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"


# Kinematics
def pk1Stress(u,c1,nu):
    G = 2*c1
    K = 2*G*(1+nu)/(3*(1-2*nu)) # Bulk modulus
    I = Identity(V.mesh().geometry().dim())  # Identity tensor
    F = I + grad(u)          # Deformation gradient ( I+ strain = L/L + dL/L = l/L = F)
    C = F.T*F                # Right Cauchy-Green tensor
    Ic, J = tr(C), det(F)    # Invariants of deformation tensors

    pk1 = 2*c1*F + (K*ln(J) - 2*c1) *inv(F).T # PK1 stress
    return pk1

def geometry_cube():
    mesh = UnitCubeMesh(4, 4, 4) 
    boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    x0 = AutoSubDomain(lambda x: near(x[0], 0))
    x1 = AutoSubDomain(lambda x: near(x[0], 1))
    y0 = AutoSubDomain(lambda x: near(x[1], 0))
    z0 = AutoSubDomain(lambda x: near(x[2], 0))
    x0.mark(boundary_parts, 1)
    y0.mark(boundary_parts, 2)
    z0.mark(boundary_parts, 3)
    x1.mark(boundary_parts, 4)
    
    # Define a MeshFunction over two subdomains
    subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())
    class Omega0(SubDomain):
        def inside(self, x, on_boundary):
            return True if x[1] <= 0.5 else False
    class Omega1(SubDomain):
        def inside(self, x, on_boundary):
            return True if x[1] >= 0.5 else False
        
    # Mark subdomains with numbers 0 and 1
    subdomain0 = Omega0()
    subdomain0.mark(subdomains, 0)
    subdomain1 = Omega1()
    subdomain1.mark(subdomains, 1)
    
    return boundary_parts,subdomains


# Create mesh and define function space ============================================
facet_function,subdomains = geometry_cube()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())

# Limit quadrature degree
dx = dx(degree=4)

# Create function space
V = VectorFunctionSpace(mesh, "CG", 1) 
du, _u = TrialFunction(V), TestFunction(V)
u = Function(V, name='u') 


# Create tensor function spaces for stress and stretch output 
W_DFnStress = TensorFunctionSpace(mesh, "CG", degree=1)
defGrad = Function(W_DFnStress, name='F')
PK1_stress = Function(W_DFnStress, name='PK1')


# Define Dirichlet boundary (symmetric Bc- no out  of plane motion)
bc0 = DirichletBC(V.sub(0), Constant(0.), facet_function, 1)# ux = 0 at x =0
bc1 = DirichletBC(V.sub(1), Constant(0.), facet_function, 2)# uy = 0 at y =0
bc2 = DirichletBC(V.sub(2), Constant(0.), facet_function, 3)# uz = 0 at z =0
x1disp = Constant(0.0)
bc3 = DirichletBC(V.sub(0), x1disp, facet_function, 4)# prescribed disp 
bcs = [bc0,bc1,bc2,bc3]


# material parameters for different subdomains, NeoHookean material
c0 = 1000  
c1 = 2000
nu =0.49


# # Total potential energy
pkstrs0 =  pk1Stress(u,c0,nu)  #stress in sub domain 0
pkstrs1 =  pk1Stress(u,c1,nu)  #stress in sub domain 1

F1_0 = inner(pkstrs0, grad(_u))*dx(0) 
F1_1 = inner(pkstrs1, grad(_u))*dx(1) 

FF = F1_0+F1_1
DJ = derivative(FF,u,du)  


# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(FF, u, bcs=bcs, J=DJ)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'


# Time stepping parameters
Tend = 1
Nsteps = 10
time = np.linspace(0, Tend, Nsteps+1)
dt = time[1]-time[0]

# Save solution in xdmf format for visualization using paraview
file_results = XDMFFile("./Results/TestCubeStretching/Uniaxial.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True


for i in range(len(time)):
    tt = time[i]
    print(i, 'time: ', tt)

    # Increase traction
    x1disp.assign(0.5*tt)

    # solve and save disp
    solver.solve()

    # Save DefGrad to file in xdmf format
    I = Identity(V.mesh().geometry().dim())  # Identity tensor
    DF = I + grad(u)
    defGrad.assign(project(DF, W_DFnStress))

    # Save Stress to file in VTK format (I have an issue here. for stress calculation what expression should I use)
    PK1t=  pk1Stress(u,c0,nu) # PK1t=  pk1Stress(u,c1,nu) which expression to use ????
    PK1_stress.assign(project(PK1t, W_DFnStress))

    # save xdmf file
    file_results.write(u, tt)
    file_results.write(defGrad, tt)
    file_results.write(PK1_stress, tt)

print('done')
2

You are projecting your stress into a CG 1 space, when the stressed are contained in a DG0 space. You should project them into an appropriate function space and use XDMFFile.write_checkpoint to write it to a suitable file format for visualization in paraview. See for instance: When using CG degree 1, the result corresponds to degree 2 - #10 by dokken

2 Likes