I have a system (here simplified) of a deformed hyperelastic 2D disk, currently under compression from two rigid plates on its top and bottom surface. I know the deformed geometry, distance between the plates, the applied force and pressure exerted. I now want to set this as the starting state and calculate the updated shape if I either decrease or increase the distance between the plates (or if it is easier, if I increase or decrease the force/preassure). I have found a workable starting code (below) from a great post and this tutorial but these consider the “disk/volume” to not be initially compressed. How would I go about implementing this initial compression condition and then deform the volume from there? Any links, code etc would be very helpful!
# %%
# https://fenicsproject.discourse.group/t/compression-of-elastic-sphere/7480
from dolfin import *
from mshr import *
from matplotlib import pyplot as plt
# Create mesh of quarter circle (to take
# advantage of symmetry).
N = 32
R = 1.0
c = Circle(Point(0,0),R,4*N)
r1 = Rectangle(Point(-2*R,-2*R),Point(0,2*R))
r2 = Rectangle(Point(-2*R,-2*R),Point(2*R,0))
mesh = generate_mesh(c-r1-r2,N)
X = SpatialCoordinate(mesh)
h = CellDiameter(mesh)
# Linear elasticity, for simplicity. Replace with a nonlinear
# model (e.g., neo-Hookean) if displacements are large.
V = VectorFunctionSpace(mesh,"CG",1)
u = Function(V)
K = Constant(1.0)
G = Constant(1.0)
I = Identity(2)
eps = sym(grad(u))
sigma = K*tr(eps)*I + 2*G*(eps-tr(eps)*I/3)
# Penalty force:
x = X + u
R = Constant(R)
E = 9*K*G/(3*K+G)
C_pen = Constant(1e3) # (Dimensionless constant)
k_pen = C_pen*E/h # (Typical choice of contact penalty)
chi_pen = conditional(gt(x[1],R),1,Constant(0))
t_pen = k_pen*chi_pen*as_vector([Constant(0),R-x[1]])
# Residual of variational problem:
v = TestFunction(V)
res = inner(sigma,grad(v))*dx - dot(t_pen,v)*ds
# Quarter-symmetry BCs, with a constant vertical displacement
# applied to the bottom.
bdry_disp = Constant(0)
bcs = [DirichletBC(V.sub(0),Constant(0),"x[0]<DOLFIN_EPS"),
DirichletBC(V.sub(1),bdry_disp,"x[1]<DOLFIN_EPS")]
# Increase displacement and output force:
N_STEP = 7
MAX_DISP = 0.25
disp_inc = MAX_DISP/N_STEP
for i in range(0,N_STEP):
# Increase the applied displacement and re-solve:
bdry_disp.assign(Constant(float(bdry_disp)+disp_inc))
solve(res==0,u,bcs=bcs)
# Assemble net penalty force to get resultant:
print("Displacement = "+str(float(bdry_disp))
+" , Force = "+str(assemble(-t_pen[1]*ds)))
# Plot results:
plot(u,mode="displacement")
plt.show()
# %%