Expression based on the volume of holes in a mesh

variational formulation
1

Hello,
I wish to simulate some deformable object with some air-chamber.
My input is a tetrahedral mesh. Starting from some surface mesh, it is first converted to volumetric mesh with tetgen, then to xdmf with meshio.
When loading the mesh, I can mark the boundaries of the air chamber, and use the measure ds_ac = Measure('ds', subdomain_data=boundaries)(1)
Till now, I use the pressure as parameter to inflate the system, I can express the force as pressure multiplied by normal.
Now, the problem is that the system is not behaving well when the pressure reaches some value, I suspect bistability is causing the solver to fail, so I want to switch to another parameter, namely the quantity of air in the air chamber, then use gas law pV=nRT.
My problem is that I had no idea how to express the volume V in the formulation. I know the function assemble could do an integration, however no domain is representing the hole as it is not meshed.

Let me give you an MWE of what I’m doing now (for simplicity I generate a simple 2D mesh here, and there is just one iteration, but in practice it will be 3D and I loop with increasing pressures)

from dolfin import *
import numpy as np 
from mshr import *

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, "eliminate_zeros": True, "precompute_basis_const": True, "precompute_ip_const": True}

# Create mesh and define function space
box = Rectangle(Point(0.0, 0.0), Point(1.0,0.1))
subboxes = [Rectangle(Point(0.1*i-0.01, 0.04), Point(0.1*i+0.01,0.06)) for i in range(1,10)]  
domain=box
for i in range(0,9):
  domain=domain-subboxes[i]
mesh=generate_mesh(domain, 100) 

V = VectorFunctionSpace(mesh, "Lagrange", 2)
N = FacetNormal(mesh)
left =  CompiledSubDomain("near(x[0], side)", side = 0.)
c = Expression(("0.0", "0.0"),degree=2)
bcl = DirichletBC(V, c, left)
bcs = [bcl]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
class myBubbles(SubDomain):
      def inside(self,x,on_boundary):
        return on_boundary and \
            not(near(x[0],0.0)) and not(near(x[0],1.0)) and \
                x[1]>0.01 and x[1]<0.09

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

# Elasticity parameters
E, nu = 100.0, 0.27
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
filex = XDMFFile('MWEresults.xdmf');
filex.parameters.update(
{
    "functions_share_mesh": True, "rewrite_function_mesh": False, "flush_output": True
})
u.rename("displacement","");
p=0.01 ##Here is the pressure, I would like to change for n (number of air molecules)
#computation of energy
F = I + grad(u)
C = F.T*F
Finv=inv(F)
Ic = tr(C)
J  = det(F)
n = J*Finv.T*N
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
pr=Constant((p))
## Here I'd like to change pr expression
## First compute volume of the bubbles 
## V_bubbles = ...
## Then express pr=n/V*Constant(R*T)
boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1, 0)
boundaries.set_all(0)
bubbles_bound=myBubbles()
bubbles_bound.mark(boundaries, 1)
ds_region = Measure('ds', subdomain_data=boundaries)
Pi = psi*dx - dot(u,-pr*n)*ds_region(1)
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
Fmat = derivative(psi, u, v)
      
# Compute Jacobian of F
J = derivative(F, u, du)
# Solve variational problem
solve(F == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options,
      solver_parameters={
        "newton_solver":{\
          "relaxation_parameter":1.0, "relative_tolerance":1e-6 ,
        "maximum_iterations":100,  'linear_solver': 'mumps'
          } })
                   
filex.write(u, p);
2

If you are solving in small increments, a linearized approximation of the volume change from increment n to increment n+1 would be

\Delta V^{n+1} ~\approx~ -\int_{S}(\mathbf{u}^{n+1} - \mathbf{u}^n)\cdot\mathbf{n}\,dS\text{ ,}

where \mathbf{u}^i is the displacement field at step i, S is the surface of the bubble, and \mathbf{n} is the outward unit normal to the solid (and thus points into the bubble, which is why there is a minus sign in front). However, you need to keep in mind that the normal vector and area element in the deformed configuration depend on \mathbf{u}, via Nanson’s formula, so

\mathbf{n}\,dS = J\mathbf{F}^{-T}\mathbf{N}\,dS_0\text{ ,}

where \mathbf{N} and dS_0 are the normal and area element in the original configuration, \mathbf{F} is the deformation gradient, and J = \operatorname{det}\mathbf{F}. (If you’re taking small enough steps to justify the linearized volume change approximation, then it’s probably reasonable to compute \mathbf{n}\,dS explicitly, using \mathbf{u}^n.)

3

Thanks, for your answer.
I have been investigating, a bit how to solve my problem.
I know I could in theory compute the volume with the divergence theorem

id=Expression(("x[0]", "x[1]", "x[2]"), degree=3)
vol = assemble(-Constant(1./3.)*inner(id,n)*ds_region(1,domain=mesh))

I suppose I can not use this expression like this, as the result is evaluation directly and is not an expression.
I found ufl has an Integral type, but this is not supposed to be used directly (TypeError: unsupported operand type(s) for /: 'Product' and 'Integral')
The solution could be to iterate like you suggest, but how to implement this exactly ?

pr=kGas*nn/vol
Pi = psi*dx  - dot(u,-pr*n)*ds_region(1)
F = derivative(Pi, u, v) +  derivative(Pi, vol, v)*dvoldu

What would be the dvoldu term (derivation of volume with respect to u ?)

4

I know I could in theory compute the volume with the divergence theorem

If you want the deformed volume, and you’re formulating the problem on the undeformed configuration, you’d want to define id with something like

X = SpatialCoordinate(mesh)
id = X + u

instead of using an Expression, where the "x[i]" components would refer to coordinates in the original configuration.

I suppose I can not use this expression like this, as the result is evaluation directly and is not an expression.

One way to include a global scalar equal to an integral would be to use a mixed space with a "Real"-type element, i.e., a single global scalar degree of freedom (DoF) in the space \mathbb{R}. This is demonstrated for applying an integral constraint with a global Lagrange multiplier here. To instead introduce a real DoF I that is enforced to equal \int_S f\,dS, you could add

+\int_S\left(\frac{I}{\vert S\vert}-f\right)\delta I\,dS

to the variational problem’s residual, where \delta I is a test function in \mathbb{R} and \vert S\vert is the measure of S, which can be computed via assemble. If this holds \forall\, \delta I\in\mathbb{R}, then I = \int_S f\,dS.

What would be the dvoldu term (derivation of volume with respect to u ?)

Keep in mind that the approach of adding an integral of \mathbf{f}\cdot\mathbf{u} to the energy functional \Pi(\mathbf{u}) only works if \mathbf{f} is a conservative force field. In particular, a follower load in the deformed normal direction is not conservative, and would need to be added directly to the residual. This may actually be why your original constant-pressure approach was causing problems. The correct way to include a follower load of the form -p\mathbf{n} would be to add

-\int_S(-p\mathbf{n})\cdot\mathbf{v}\,dS

directly to the residual F, since the derivative of \int_Sp\mathbf{n}\cdot\mathbf{u}\,dS with respect to \mathbf{u} would include extra terms from the \mathbf{u}-dependence of \mathbf{n}\,dS.

5

Thank you very much !