Nonlinear problem

Question

Hi all, I want to solve the hyperelasticity problem from the demo. The cube is fixed on one face and a Traction force is applied to the opposite face. But I want to apply this traction force by increments Traction(t). For this I created a time loop.

Please could you take a look at my code. I am not sure that is correct. Do I need to update my solution after each timestep?

Thanks!

from dolfin 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
mesh = UnitCubeMesh(16, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

# Initialize sub-domain instances
left = Left()
right = Right()

#New measures
domains = CellFunction("size_t", mesh)
boundaries = FacetFunction("size_t", mesh)

left.mark(boundaries, 1)
right.mark(boundaries, 3)

dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define Dirichlet boundary (x = 0)
c = Expression(("0.0", "0.0", "0.0"))
bcl = DirichletBC(V, c, boundaries, 1)
bcs = [bcl]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
Traction  = Expression(("t", "0.0", "0.0"),t=0.0)  # Traction force on the boundary

# Kinematics
I = Identity(len(v))        # 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 = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# 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 - dot(Traction, u)*ds(3)

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)



# Timestepping
t = 0.0; dt = 0.01; T = 0.05
inc = 0
file_u = File("displacement.pvd");

while t < T:
        inc += 1
        print 'Increment =', inc
        print 'time =', t
        t += dt
        Traction.t=t

        # Solve variational problem
        solve(F == 0, u, bcs, J=J,
          form_compiler_parameters=ffc_options)

        # Save solution in VTK format
        file_u << u;

# Plot and hold solution
plot(u, mode = "displacement", interactive = True)

Nonlinear problem

1 Answer