Occlude a cylinder with suction pressure

variational formulation
1

I am trying is to occlude a cylinder with internal suction pressure (the aim is to occlude it as much as possible). I am using purely incompressible hyperelastic material with the 2-field formulation in Fenics. However, I can not occlude the cylinder after a certain limit. The figures below are the initial configuration and the deformed configuration of the cylinder.

1
11454×922 31.9 KB
2
21454×922 81.1 KB

I can not occlude the cylinder any further. However, if I use other FE packages, e.g., FEBio, I can almost fully occlude the cylinder. I am wondering what could be the possible reason that I can not occlude it further.

There are 2-3 things that might be wrong. First, although I am applying pressure normal to the inner face using the FacetNormal function, after a certain stage of deformation the pressure is not acting in the normal direction. If that is the case somehow we need to calculate the FacetNormal function in every step of deformation (do not know how to do that). Second, I am using the 2-field formulation whereas FeBio uses the 3-field formulation. I really do not know whether that would be helpful or not. But, I think the 2-field formulation should be fine.

Could you please help me to find the issue?

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

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"


# Kinematics
def pk1Stress(u,c1,pressure):
    I = Identity(V.mesh().geometry().dim())  # Identity tensor
    F = I + grad(u)          # Deformation gradient
    C = F.T*F                # Right Cauchy-Green tensor
    Ic = tr(C)               # Invariants of deformation tensors
    J = det(F)
    E = 0.5*(C-I)
    pk1 = 2*c1*F-pressure*inv(F).T
    return pk1, (J-1)

## geometry
def geometry_3d():
    c_length = 2.0; # L
    c_inner_radius = 1.0; # A
    c_outer_radius = 2.0; # B 
    mesh = Mesh('./QuarterCylinder.xml')

    boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    x0 = AutoSubDomain(lambda x: np.absolute(x[0]) < 0.01)
    y0 = AutoSubDomain(lambda x: np.absolute(x[1]) < 0.01)
    z0 = AutoSubDomain(lambda x: near(x[2], 0))
    innerface = AutoSubDomain(lambda x: (np.sqrt(x[0]*x[0]+x[1]*x[1])- c_inner_radius) <= 0.1  )       #..
    x0 .mark(boundary_parts, 1)
    y0 .mark(boundary_parts, 2)
    z0 .mark(boundary_parts, 3)
    innerface.mark(boundary_parts, 4)
    
    return boundary_parts


# Create mesh and define function space
facet_function = geometry_3d()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure("dx")
ds = Measure("ds", subdomain_data=facet_function, subdomain_id=4)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())

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

# Create function space (in mixed space approach, have to define elements)
element_u = VectorElement("CG", mesh.ufl_cell(), 2) 
element_p = FiniteElement("CG", mesh.ufl_cell(), 1) 

mixed_element = MixedElement([element_u,element_p]) 
V = FunctionSpace(mesh, mixed_element) 


# Define test and trial functions
dup = TrialFunction(V)  
_u_p = Function(V)      
u, p = split(_u_p) 
_u, _p = TestFunctions(V)

# Define Dirichlet boundary
bc0 = DirichletBC(V.sub(0).sub(0), Constant(0.), facet_function, 1)
bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facet_function, 2)
bc2 = DirichletBC(V.sub(0).sub(2), Constant(0.), facet_function, 3)
bcs = [bc0,bc1,bc2]

# material parameters
c1 = 1000.0

# pressure on the inner face of the cylinder
facenormal = FacetNormal(mesh) 
innerpress = Constant(-0.0)
I = Identity(V.mesh().geometry().dim()) 
F = I + grad(u)
press_trac = det(F)*dot(inv(F).T, innerpress*facenormal) # pressure 


# # Total potential energy
I = Identity(V.mesh().geometry().dim())
pk1, p_bar =  pk1Stress(u,c1,p) 
dgF = I + grad(u) 
F1 = inner(pk1, grad(_u))*dx - dot(press_trac, _u)*ds
F2 = p_bar*_p*dx        
F = F1+F2
DJ = derivative(F, _u_p,dup)


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

# Extract solution components
u, p = _u_p.split() # seperate functions
u.rename("u", "displacement")
p.rename("p", "pressure")

# Save solution in VTK format
file_results = XDMFFile("./Results/TestTriaxialLoading/Triaxial_3field.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

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

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

    # Increase traction
    innerpress.assign(2900.0*time[i]) #Pa

    # solve and save disp
    solver.solve()

    # save xdmf file
    file_results.write(u, time[i])

print('done')
2

Equal-order interpolation of displacements and pressures is not a stable discretization of the two-field formulation. By analogy to the linearized problem (which is formally equivalent to Stokes flow), you might have better luck with the inf-sup-stable Taylor–Hood element, using quadratic elements for displacement (but still linear elements for the pressure):

# Create function space (in mixed space approach, have to define elements)
#element_u = VectorElement("CG", mesh.ufl_cell(), 1) 
element_u = VectorElement("CG", mesh.ufl_cell(), 2) 
element_p = FiniteElement("CG", mesh.ufl_cell(), 1)
3

I tried the same simulation with

element_u = VectorElement("CG", mesh.ufl_cell(), 2)

The solver failed again. Should I use hex elements instead of tet?

Here is a link to my cylinder mesh file. You may run the code using this mesh file

4

Hi, there are various things you can try such as taking smaller load steps, switching to PETScSNESSolver with line search.

5

I tried to change the solver with this

solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'petsc'  # 'mumps', 'lu'

This one is not helping. The solver fails after few steps. How can I use PETScSNESSolver with line search.

Also, I tried smaller loading steps. That is not helpful. Looks like its a volume locking problem.

6

Dear @dokken,

Could you please help me with it, as I am still stuck? I tried a smaller loading step, PETSC solver, and fine linear tetra mesh too. Also, I tried to import quadratic tetra and linear Hex mesh to FEnics but experiencing errors ( using Gmsh for mesh generation).

Your guidance would be appreciated.

Thank you

7

Unstructured Hex meshes aren’t supported in FEniCS, but as I understand the support for these has improved in dolfinx. Your mesh-to-be-imported in FEniCS has to be linear. You can then use a higher order function space over the linear mesh – since these two are decoupled in FEniCS.

As for using SNES, you can either use it through the NonlinearVariationalSolver interface via

snes_solver_parameters = {
    "nonlinear_solver": "snes",
    "snes_solver": {
        "linear_solver": "mumps",  
        "line_search": "bt", # stands for backtracking line search
        "maximum_iterations": 50,
        "report": True
    }
}
solver.parameters.update(snes_solver_parameters)

or for a finer control use PETSc.SNES directly as demonstrated here and in several other posts on the group.

2 Likes