Hello’ I am trying to model a cantilever under distributed load considering large deformation assumption and voigt notation, but when I run the code, I received the following error message:
ValueError: too many values to unpack (expected 2)
I think it would be related to nonlinear terms that are existed in strain tensor!
I would be appreciated if anyone can help me about this issue.
here is the entire code:
from dolfin import *
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.patches import Rectangle
import numpy as np
from ufl.core.multiindex import indices
Problem definition
Len = 10;
W = 1;
t=0.1;
Load=1150000.0
Mesh for cantilever beam
mesh = BoxMesh(Point(0, 0, 0), Point(Len, W, t), 60, 4, 4)
mesh.init()
V = VectorFunctionSpace(mesh, ‘P’, 2) # displacements
Define boundary condition
Left boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
Right boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], Len)
Top boundary
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], t)
#Boundary segments
left = Left()
top_boundary = TopBoundary()
boundaries = MeshFunction(“size_t”, mesh,0)
boundaries.set_all(0)
left.mark(boundaries, 1)
top_boundary.mark(boundaries, 2)
ds= Measure(“ds”, domain=mesh, subdomain_data=boundaries)
define variational variables
du=TrialFunction(V) # Incremental displacement
del_u=TestFunction(V) # Test function
u=Function(V) # Displacement
Geometric dimension
d = u.geometric_dimension()
d x d identiy matrix
Kronecker delta
I = Identity(d)
#index notation
i,j,k,l= indices(4)
----------------------Parameters------------
c_11 = 2.695652e9
c_12 = 1.65217e9
c_13 = 1.65217e9
c_33 = 2.695652e9
c_22= 2.695652e9
c_23= 1.65217e9
c_55= 1.04347e9
c_44 = 1.04347e9
c_66= 1.04347e9
#---------------------------------------------
#------------------Elasticity tensor c^E
C = as_tensor([[c_11, c_12, c_13, 0., 0., 0.],
[c_12, c_22, c_23, 0., 0., 0.],
[c_13, c_23, c_33, 0., 0., 0.],
[0., 0., 0., c_44, 0., 0],
[0., 0., 0., 0., c_55, 0.],
[0., 0., 0., 0., 0., c_66]])
#-----------------rewrite the tensor into a vector using Voigt notation
def strain_voigt(ten):
return as_vector([ten[0,0],ten[1,1],ten[2,2],2ten[1,2],2ten[0,2],2*ten[0,1]])
#------------------rewrite a vector into a tensor using Voigt notation
def stress_voigt(vec):
return as_tensor([[vec[0], vec[5], vec[4]],
[vec[5], vec[1], vec[3]],
[vec[4], vec[3], vec[2]]])
#----------------define the function for the strain tensor in fenics notation
def epsilon(u):
return as_tensor((1./2.(u[i].dx(j)+u[j].dx(i))+(u[k].dx(i))(u[k].dx(j))),[i,j])
return as_tensor((1./2.*(u[i].dx(j)+u[j].dx(i))),[i,j])
def sigma(u):
return stress_voigt((dot(C, strain_voigt(epsilon(u)))))
loading in the surface normal direction , the traction vector in Pa
tr=Constant((0.0,0.0,Load))
Body Force , the gravity acceleration in m/s^2
bf=Constant((0,0,0))
left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left)]
variational formulation
Form=sigma(u)[i,j]*epsilon(del_u)[i,j]*dx-bf[j]*del_u[j]*dx + tr[j]*del_u[j]*ds(2)
Gain=derivative(Form,u,du)
solve(Form == 0, u, bc, J=Gain,solver_parameters={“newton_solver”:{“linear_solver”: “mumps”, “absolute_tolerance”: 1E-9, “relative_tolerance”: 1E-3, “maximum_iterations”: 20, “error_on_nonconvergence”: True } },form_compiler_parameters={“cpp_optimize”: True, “representation”: “quadrature”, “quadrature_degree”: 2} )
print(‘min/max u:’,u.vector().get_local().min(),u.vector().get_local().max())