Storing traction and displacement in each iteration

1

Hello Everyone,

In my displacement controlled code below, I tried to store the calculated traction and displacement in each iteration to create the Force-Displacement diagram, but I could not find a way to do this.
If you have any suggestions, it would be very helpful.
Thanks for your help in advance!

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# 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 = Mesh("FinerTriangularMesh05.xml")
V = VectorFunctionSpace(mesh, "Lagrange", 2)
def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)

bc_bot = DirichletBC(V, [0,0,0], boundary_bot)
def boundary_top(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[1], 9, tol)

bc_top = DirichletBC(V, [0.0,-1.8,0], boundary_top)
bcs = [bc_bot, bc_top]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

T  = Constant((0.0, 0.0, 0.0))     # Traction force on the boundary

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)


dss = ds(subdomain_data=boundary_subdomains)


# Kinematics
d = u.geometric_dimension()
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)

# 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


TS = TensorFunctionSpace(mesh, "Lagrange", 2)

# Total potential energy
Pi = psi*dx 

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

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

# Solve variational problem
solve(dPi == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)
storedu=[]
storedu.append(u)
print(storedu)


#Piola-Kirchhoff Stress
def S(u):
    pi = psi
    stress = diff(psi, F)
    return stress

#Traction
Stress = project(S(u), TS)

traction = Function(V)
e2 = Constant((0, 1, 0))
traction = project(dot(Stress,e2), V)
storedtraction = assemble(traction[1]*dss) #top surface traction
Traclist.append(storedtraction)
print(Traclist)
print(storedtraction)



# Save solution in VTK format
file = File("displacement.pvd");
file << u;

file = File("stress.pvd");
file << Stress;

file = File("traction.pvd");
file << traction;

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

print("Done")

# plot(storedtraction,storedu)
# plt.show()

[quote="HaticeGokcenGuner, post:1, topic:7351, full:true"]
Hello Everyone,

In my displacement controlled code below, I tried to store the calculated traction and displacement in each iteration to create the Force-Displacement diagram, but I could not find a way to do this. 
If you have any suggestions, it would be very helpful.
Thanks for your help in advance!

---
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# 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 = Mesh("FinerTriangularMesh05.xml")
V = VectorFunctionSpace(mesh, "Lagrange", 2)
def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)

bc_bot = DirichletBC(V, [0,0,0], boundary_bot)
def boundary_top(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[1], 9, tol)

bc_top = DirichletBC(V, [0.0,-1.8,0], boundary_top)
bcs = [bc_bot, bc_top]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

T  = Constant((0.0, 0.0, 0.0))     # Traction force on the boundary

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)


dss = ds(subdomain_data=boundary_subdomains)


# Kinematics
d = u.geometric_dimension()
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)

# 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


TS = TensorFunctionSpace(mesh, "Lagrange", 2)

# Total potential energy
Pi = psi*dx 

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

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

# Solve variational problem
solve(dPi == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)
storedu=[]
storedu.append(u)
print(storedu)


#Piola-Kirchhoff Stress
def S(u):
    pi = psi
    stress = diff(psi, F)
    return stress

#Traction
Stress = project(S(u), TS)
traction = Function(V)
e2 = Constant((0, 1, 0))
traction = project(dot(Stress,e2), V)
storedtraction = assemble(traction[1]*dss) #top surface traction
Traclist.append(storedtraction)




# Save solution in VTK format
file = File("displacement.pvd");
file << u;

file = File("stress.pvd");
file << Stress;

file = File("traction.pvd");
file << traction;

print("Done")

# plot(storedtraction,storedu)
# plt.show()
2

As I understand the question, your goal is to produce a simulated force-displacement curve. In your code, you only solve for a single displacement BC, which will only produce one data point, not an entire curve. I think what you want to do is loop through several different values of the applied displacement, call solve repeatedly for each step of the loop, and store the result for the total force on the displacement boundary. Consider the code below.

This code computes the total force on the +y surface of a box. Note the following:

  1. I have removed the projection of the stress into a CG-2 tensor function space and the subsequent projection of the traction into a CG-2 vector function space. This is not necessary to compute the total force on the top surface.
  2. The form passed to assemble should have dss(1), not dss, since you only want to integrate over the portion of the boundary marked 1.
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# 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 = BoxMesh(Point(0,0,0), Point(1,9,1), 1, 9, 1)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
def boundary_bot(x, on_boundary):
    tol = 1E-7
    return on_boundary and near(x[1], 0, tol)

bc_bot = DirichletBC(V, [0,0,0], boundary_bot)
def boundary_top(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[1], 9, tol)

max_disp = 1.8
disp = Constant([0.0,0.0,0.0])
bc_top = DirichletBC(V, disp, boundary_top)
bcs = [bc_bot, bc_top]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

T  = Constant((0.0, 0.0, 0.0))     # Traction force on the boundary

boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(boundary_top).mark(boundary_subdomains, 1)


dss = ds(subdomain_data=boundary_subdomains)


# Kinematics
d = u.geometric_dimension()
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)

# 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 

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

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

#Piola-Kirchhoff Stress
stress = diff(psi, F)

force_FS = VectorFunctionSpace(mesh, "Real", 0)
df = TestFunction(force_FS)

# Solve variational problem
storedu=np.linspace(0, max_disp, 10)
Traclist = []
normal = Constant((0,1,0))
file = File("displacement.pvd")
for d in storedu:
    print(d)
    disp.assign(Constant([0.0, d, 0.0]))
    solve(dPi == 0, u, bcs, J=Jac,
          form_compiler_parameters=ffc_options)
    file << u
    storedtraction = assemble(inner(df, dot(stress,normal))*dss(1)) #top surface traction
    Traclist.append(storedtraction[1])

from matplotlib import pyplot as plt
plt.plot(storedu, E*storedu/9, '-', storedu, Traclist, '-o')
plt.xlabel("Displacement")
plt.ylabel("Force")
plt.legend(["Linear","Neo-Hookean"])

Producing:
force-displ

3 Likes
3

Thank you very much for your great help!

2 Likes