Apply normal suction pressure on the top surface and deformed surface

dolfinx
1

I want to apply a hypobaric pressure on the top surface of a plate. The pressure is applied on the circular area. The mesh deforms and it grows. The pressure now is applied in the vertical direction, but not normal to the surface.

Please help.

import numpy as np
import ufl

from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
from dolfinx import log
from dolfinx.io import XDMFFile

L = 32.0e-3  # length of the plate in meters
W = 32.0e-3  # width of the plate in meters
H = 1.53e-3  # height of the plate in meters

# Creating mesh for the plate - increased resolution
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, W, H]], [20, 20, 4], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 2))

fdim = domain.topology.dim - 1

# The circular region of radius 8mm in the center of the plate
def circular_region(x):
    r = 8.0e-3
    return np.sqrt((x[0]-L/2)**2 + (x[1]-W/2)**2) <= r

# Boundary conditions
def boundary(x):
    return np.logical_and(np.isclose(x[2], H), np.logical_not(circular_region(x)))
# Boundary conditions
def circular_boundary(x):
    return np.logical_and(np.isclose(x[2], H), circular_region(x))
   
# Define all side faces excluding bottom and the top face
def side_faces(x):
    return np.logical_or.reduce([np.isclose(x[0], 0.0), np.isclose(x[0], L), np.isclose(x[1], 0.0), np.isclose(x[1], W)])

# Identify and mark side faces
side_facets = mesh.locate_entities_boundary(domain, fdim, side_faces)
# Identify and mark top face except circular region
top_facets = mesh.locate_entities_boundary(domain, fdim, boundary)
# Identify and mark circular region on top face
circular_facets = mesh.locate_entities_boundary(domain, fdim, circular_boundary)

# Concatenate and sort the arrays based on facet indices. Marked facets: side - 1, top - 2, circular - 3
marked_facets = np.hstack([side_facets, top_facets, circular_facets])
marked_values = np.hstack([np.full_like(side_facets, 1), np.full_like(top_facets, 2), np.full_like(circular_facets, 3)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

# Apply Dirichlet boundary condition
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)

# Apply boundary conditions for side faces and the top
side_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
top_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(2))
bcs = [fem.dirichletbc(u_bc, side_dofs, V), fem.dirichletbc(u_bc, top_dofs, V)]

# Define constants
B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

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

# Elasticity parameters for the Ogden model
mu_value = 0.03e6  # in Pascals
alpha_value = 13.57
mu = fem.Constant(domain, mu_value)
alpha = fem.Constant(domain, alpha_value)

# First Piola-Kirchhoff stress tensor
P = mu * (F * ufl.exp(alpha * (Ic - 3)) - F.T * ufl.exp(-0.5 * alpha * (Ic - 3)))

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(3)  # Pressure applied to circular region

problem = fem.petsc.NonlinearProblem(F, u, bcs)

from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-4
solver.rtol = 1e-4
solver.convergence_criterion = "incremental"

log.set_log_level(log.LogLevel.INFO)
# Write the mesh to a file once
with XDMFFile(MPI.COMM_WORLD, "deformed_mesh.xdmf", "w") as outfile:
    outfile.write_mesh(domain)

max_pressure = 10.0  # Maximum pressure (in PSI)
num_steps = 10  # Number of steps

for n in range(1, num_steps+1):
    pressure = n * max_pressure / num_steps
    T.value[2] = pressure
    num_its, converged = solver.solve(u)
    assert(converged)
    u.x.scatter_forward()

    print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
   
    # Open the XDMF file in append mode and write the current function u to the file
    with XDMFFile(MPI.COMM_WORLD, "deformed_mesh.xdmf", "a") as outfile:
        outfile.write_function(u, n)

outfile.close()
2

First of all, please format your code with 3x` encapsulation, ie.

```python
#add code here

```
1 Like
3

Done. Thank you. Looking for your support.

4

Just to confirm I understand your question: you wish for the traction \boldsymbol{T} to be normal to the (deformed) surface of the plate? I.e.

\boldsymbol{T} = p \boldsymbol{n}

where \boldsymbol{n} is the outward normal of the plate in the deformed configuration. If that’s the case, you can transform the integral using Nanson’s formula, \boldsymbol{n} ds = J \boldsymbol{F}^{-T} \boldsymbol{n}^{(0)} ds^{(0)}, i.e.

\int_{\partial \Omega_3} {\boldsymbol{T} \cdot \boldsymbol{v} ds} = \int_{\partial \Omega_3} {(p \boldsymbol{n} ds) \cdot \boldsymbol{v}} = \int_{\partial \Omega^{(0)}_3} {(p J \boldsymbol{F}^{-T} n^{(0)} ds^{(0)}) \cdot \boldsymbol{v}}

where all quantities with superscript “(0)” are expressed in the reference configuration. You can implement this in your weak form as

Finv = ufl.inv(F)
p = fem.Constant(domain, PETSc.ScalarType(0))
n = ufl.FacetNormal(domain)
F = (
    ufl.inner(ufl.grad(v), P)*dx
    - ufl.inner(v, B)*dx
    - ufl.inner(p * J * ufl.dot(n, Finv), v)*ds(3)
)

Note that p in this implementation is expressed relative to the deformed configuration.

1 Like