Hello
I am trying to implement a dynamic hyperelasticity problem (Neohookean Model). A static demo is available here
My question is all about implementation of the variational formulation. I looked at this document (page.9) and here is what I did for a time dependent hyperelasticity problem:
discretization in time:
Lets say \rho=1 and Stress is defined as:
Where J = det(F) and b is same as C in the FEniCS demo (Right Cauchy-Green tensor)
Finally after multiplying a test function (v) and integration by part, here is the variational form :
Here is a minimal code:
from dolfin import *
import numpy as np
A1 = [0,2,4,6,8]
A2 = [0,3,6,9,12]
mesh1 = Mesh("mesh.xml")
facets1 = MeshFunction("size_t", mesh1, "mesh_facet_region.xml")
domains1 = MeshFunction("size_t", mesh1, "mesh_physical_region.xml")
#function space
V = VectorFunctionSpace(mesh1, "Lagrange", 1)
W_const = Constant((0.0, 0.0))
u1= interpolate(W_const, V)
u0= interpolate(W_const, V)
dt = 0.5
t = 0
T = 2
# Define Dirichlet boundary conditions
U_const = Constant((0.0, 0.0))
bc_top = DirichletBC(V, U_const, facets1,35)
bcs = [bc_top]
# Define functions
du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Elasticity parameters
E, nu = 2E6, 0.37
mu = Constant(E/(2*(1 + nu)))
lmbdaa= Constant(E*nu/((1 + nu)*(1 - 2*nu)))
sigma = inv(J)*(lmbdaa*ln(J)*I+mu*(C_1-I))
# Define vectorial body force
delta_A2 = Expression(' alpha*t',degree=1, alpha=1, t=0)
delta_A1 = Expression(' beta*t',degree=1, beta=1, t=0)
force = (delta_A2 - delta_A1)
def BF():
return as_vector([force] + [0.0])
# Variational Form. Is that right?!
F = inner(u,v)*dx + dt*dt*inner(sigma, grad(v)) * dx -dt*dt*inner(v, BF()) * dx - 2*inner(u1,v)*dx+inner(u0,v)*dx
counter = 0
while t <= T:
delta_A1.t = A1[counter]
delta_A2.t = A2[counter]
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
# Solve variational problem
solver.solve()
u0.assign(u1)
u1.assign(u)
t += dt
counter = counter + 1
As I said this code works without error but I am not sure about the variational formulation. Is this variational form correct or I should take a new approach?