I am having a problem updating the pressure normals when applying pressure to the face of a beam. After the beam is deformed, the pressure continues to be applied in the same direction it was being applied before the beam was deformed. I have tried using Nanson’s formula n_deformed = J * inv(F).T * n_initial to update these but it does not seem to work. There’s no errors, just incorrect output. Any guesses on how to make this right?
Here the pressure is defined:
// Pressure
auto V_Pr = std::make_shared<dolfinx::fem::FunctionSpace>(
fem::create_functionspace(Functionspace, "Pr", mesh));
auto Pr = std::make_shared<dolfinx::fem::Function<PetscScalar>>(V_Pr);
Pr->name = "Pr";
Here its applied:
if(bnd_type=="Pressure")
{
std::vector<std::int32_t> facetst_values(bcst_facets.size(), itype32);
vst_facets.insert(vst_facets.end(), bcst_facets.begin(), bcst_facets.end());
vst_values.insert(vst_values.end(), facetst_values.begin(), facetst_values.end());
auto entity_to_cell = mesh->topology().connectivity(2, 3);
std::int32_t num_local_cells = mesh->topology().index_map(3)->size_local();
for (std::size_t i = 0; i < bcst_facets.size(); ++i)
{
auto cell = entity_to_cell->links(bcst_facets[i]);
if (cell[0] < num_local_cells) Pr_cells.push_back(cell[0]);
This is a modified version of the hyper elasticity Fenics demo, with inclusion of Pressure surface integrals, based on deformed configuration.
# UFL input for hyperleasticity
# =============================
#
# The first step is to define the variational problem at hand. We define
# the variational problem in UFL terms in a separate form file
# :download:`hyperElasticity.py`.
#
# We are interested in solving for a discrete vector field in three
# dimensions, so first we need the appropriate finite element space and
# trial and test functions on this space::
from ufl import (Coefficient, Identity, TestFunction, TrialFunction,
VectorElement, derivative, det, diff, dx, grad, ln,
tetrahedron, tr, variable)
# Function spaces
element = VectorElement("Lagrange", tetrahedron, 1)
mesh = Mesh(element)
n = FacetNormal(mesh)
# Trial and test functions
du = TrialFunction(element) # Incremental displacement
v = TestFunction(element) # Test function
# Note that ``VectorElement`` creates a finite element space of vector
# fields. The dimension of the vector field (the number of components)
# is assumed to be the same as the spatial dimension (in this case 3),
# unless otherwise specified.
#
# Next, we will be needing functions for the boundary source ``B``, the
# traction ``T`` and the displacement solution itself ``u``::
# Functions
u = Coefficient(element) # Displacement from previous iteration
# B = Coefficient(element) # Body force per unit volume
# T = Coefficient(element) # Traction force on the boundary
# Now, we can define the kinematic quantities involved in the model::
# Kinematics
d = len(u)
I = Identity(d) # Identity tensor
F = variable(I + grad(u)) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
nup=n*J*inv(F.T)
nup_norm = nup/inner(nup,nup)
# Before defining the energy density and thus the total potential
# energy, it only remains to specify constants for the elasticity
# parameters::
# Elasticity parameters
E = 10.0
nu = 0.3
mu = E/(2*(1 + nu))
lmbda = E*nu/((1 + nu)*(1 - 2*nu))
# Both the first variation of the potential energy, and the Jacobian of
# the variation, can be automatically computed by a call to
# ``derivative``::
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx - inner(P*nup_norm, u)*ds
# First variation of Pi (directional derivative about u in the direction of v)
F_form = derivative(Pi, u, v)
# Compute Jacobian of F
So, a pressure follower load by design is not conservative, and thus does not have a potential like you’ve written it out here. You can show that there exists a potential p* \int_{Omega_0} J dV if the pressure acts on a boundary surface enclosing a region \Omega_0.
So, either change your potential term, or just add a term \int_{Gamma_0} p J F^{-T} * v dA to your F_form. That should work.
We are trying to apply a pressure on one face of the beam (in the deformed configuration). So we want to use the normals of the deformed configuration inside a Non-Linear Newton-Rapson iteration.
I saw two errors in the .py code I sent (normalization of n, and the sign in one of the components of Pi). Either way, the results are still inaccurate.
(...)
n0=FacetNormal(mesh)
#Updated Normals
nup = J*(inv(F).T)*n
nup_norm = nup/sqrt(inner(nup,nup))
# Elasticity parameters
E = 10.0
nu = 0.3
mu = E/(2*(1 + nu))
lmbda = E*nu/((1 + nu)*(1 - 2*nu))
# Both the first variation of the potential energy, and the Jacobian of
# the variation, can be automatically computed by a call to
# ``derivative``::
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy
Pi = psi*dx + inner(P*nup_norm, u)*ds
# First variation of Pi (directional derivative about u in the direction of v)
F_form = derivative(Pi, u, v)
F_form += P*J*dot(inv(F).T*n0, v)*ds
# Compute Jacobian of F
J_form = derivative(F_form, u, du)
As I said, this kind of load does not have a potential like this. So the term inner(P*nup_norm, u)*ds is wrong - especially since you’re adding the virtual work of the load directly to your F_form, it does not make sense to include it in the potential.
So modify the potential to just Pi = psi*dx.
Also I think that there is no need to normalize the facet normal field; it should be normalized already.