Nonlinear Timoshenko beam implementation
First of all, import os
import os
followed by dolfin
from dolfin import *
Set the form compiler
parameters.
The representation
uflacs
is needed because
it correctly handles the derivation of conditional UFL expressions
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["representation"] = "tsfc"
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
Limit the quadrature order
parameters["form_compiler"]["quadrature_degree"] = 3
Import the rotation handling library
import RotVector
and the customized assembly routines
from rotation_assembling import *
Build the mesh with 200 elements along a straight line of lenght 10
beam_len = 10.
mesh = IntervalMesh(200, 0., beam_len)
Define two Lagrange (CG
) vector element of size
equal to 3.
The first one has a degree
equal to 2, and will be used to approximate
the beam displacement. The second one has a degree
equal to 1,
and will be used to approximate rotations. This choice greatly alleviates, if
not eliminates completely, the problem of parasitic shear locking,
at the expense of an higher number on unknowns.
U_ELEMENT = VectorElement("CG", mesh.ufl_cell(), 2, 3)
PHI_ELEMENT = VectorElement("CG", mesh.ufl_cell(), 1, 3)
Define a Mixed element
comprising the displacement and rotation elements
UNKS_ELEMENT = MixedElement(U_ELEMENT, PHI_ELEMENT)
and the approximating function space of the problem
UNKS_SPACE = FunctionSpace(mesh, UNKS_ELEMENT)
Define the Functions
u
to store the total displacement and rotation,uinc
to store the increment of the unknowns during the nonlinear iterative solution procedure
u = Function(UNKS_SPACE)
uinc = Function(UNKS_SPACE)
Define also the TrialFunction
and the TestFunction
du = TrialFunction(UNKS_SPACE)
v = TestFunction(UNKS_SPACE)
Take a reference to the subspace basis function values using the UFL split
method
u_u, u_phi = split(u)
du_u, du_phi = split(du)
v_u, v_phi = split(v)
The helper class SwitchNodes
is used to ease the computation
of the complement rotation.
The class stores a reference to
the unknown function u
, and a local function utmp
, that should be
defined on the same space of u.
It also stores a reference to the Expression
RotVector.cppcode_unwrap_phi
,
and assigns to its public member phi_function
the function utmp
.
In this way, whenever the eval
method of RotVector.cppcode_unwrap_phi
,
it will be applied to utmp
.
When the Switch
metod is called the values of u
are copied into utmp
. After that, the values of u
are overwritten by
evaluating the expression cppcode_unwrap_phi
.
Note that applying cppcode_unwrap_phi
to u
and not to its copy utmp
would not work because it would at the same time overwrite the content of u
while looping in interpolate
.
class SwitchNodes:
def __init__(self, u):
self.u = u
self.utmp = Function(u.function_space())
self.unwrap_nodes = CompiledExpression(compile_cpp_code(RotVector.cppcode_unwrap_phi).unwrap_phi(),
element=u.function_space().ufl_element())
self.unwrap_nodes.phi_function = self.utmp.cpp_object()
def Switch(self):
self.utmp.vector()[:] = self.u.vector()
self.u.interpolate(self.unwrap_nodes)
swe = SwitchNodes(u)
The Dirichelet boundary condition is a clamp at the origin: both the displacement and the rotation are constrained to be zero.
left = DirichletBC(UNKS_SPACE, (0,0,0,0,0,0), "fabs(x[0]) <= DOLFIN_EPS")
bc = [left]
A concentrated load is applied at the other end of the beam.
To this end, the Right
class is defined; its inside
method should return True
only fo a point on the boundary
and with the \(x_0\) coordinate equal to the beam lenght.
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], beam_len)
right = Right()
Since facets are actually points for an IntervalMesh
a FacetFunction
is used to mark the end of the beam and define the Measure
ds
:
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
right.mark(boundaries, 1)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh)
Define the rotation tensor alpha
and the rotation differential tensor Gamma
as a function of the the rotation vector u_phi
:
alpha = RotVector.RotTensor(u_phi)
Gamma = RotVector.RotDiffTensor(u_phi)
Define the linear and angular strains epsilon
and beta
, i.e.
\(\boldsymbol{\varepsilon} =
\boldsymbol{\alpha}'^T\boldsymbol{x}'_{/s}-
\boldsymbol{\alpha}^T\boldsymbol{x}'_{/s}\)
and
\(\boldsymbol{\beta} =
\boldsymbol{\varPhi}^T\mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^T)\):
epsilon = dot(alpha.T, grad(u_u)[:, 0] + as_vector([1., 0., 0.])) - as_vector([1., 0., 0.])
beta = as_vector(dot(Gamma.T, grad(u_phi))[:,0])
Define the finale values of the ditrubuted force and coule per unit of beam lenght and of the concentrated force and couple applied at the beam end
distributed_final_force = Constant([0., 0., 0.])
distributed_final_moment = Constant([0., 0., 0.])
concentrated_final_force = Constant([0., 0., 50.])
concentrated_final_moment = Constant([0., 0., 200. * pi])
the current value of the loads
distributed_force = Constant([0., 0., 0.])
distributed_moment = Constant([0., 0., 0.])
concentrated_force = Constant((0., 0., 0.))
concentrated_moment = Constant([0., 0., 0.])
and the four second order constant tensors \(\overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}\), \(\overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\beta}}\), \(\overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\varepsilon}}\) and \(\overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\beta}}\) defining the beam section linear constitutive law
k_ee = Constant([[1.E4, 0., 0.], [0., 1.E4, 0.],[0., 0., 1.E4]])
k_eb = Constant([[0., 0., 0.], [0., 0., 0.],[0., 0., 0.]])
k_bb = Constant([[1.E2, 0., 0.], [0., 1.E2, 0.],[0., 0., 1.E2]])
k_be = k_eb.T
The internal actions \(\boldsymbol{T}\) and \(\boldsymbol{M}\) are readily computed as
traction = dot(k_ee, epsilon) + dot(k_eb, beta)
moment = dot(k_be, epsilon) + dot(k_bb, beta)
The variation \(\delta\boldsymbol{\varepsilon}\) and \(\delta\boldsymbol{\beta}\) of the strain vectors \(\boldsymbol{\varepsilon}\) and \(\boldsymbol{\beta}\) can be computed as
delta_epsilon = derivative(epsilon, u, v)
delta_beta = derivative(beta, u, v)
Finally, the expresison of the functional \(\int_L(\delta\boldsymbol{\varepsilon}\boldsymbol{T}+\delta\boldsymbol{\beta}\boldsymbol{M})\mathrm{d}s- \int_L (\delta\boldsymbol{x}'\boldsymbol{t} + \boldsymbol{\varphi}_\delta\boldsymbol{m}) \mathrm{d}s - \sum_i\delta\boldsymbol{x}'_i\boldsymbol{F}_i - \boldsymbol{\varphi}_{\delta i}\boldsymbol{C}_i\) is easily written
functional = inner(delta_epsilon, traction) * dx + \
inner(delta_beta, moment) * dx - \
inner(v_u, distributed_force) * dx - \
inner(v_phi, dot(Gamma.T, distributed_moment)) * dx - \
inner(v_u, concentrated_force) * ds(1) - \
inner(v_phi, dot(Gamma.T, concentrated_moment)) * ds(1)
and the Jacobian matrix is defined by deriving it with respect to the unknown Function
u
J = derivative(functional, u, du)
The optput will be written to two files, the first for the displacement field (u.sub(0)
)
and the second for the rotation field (u.sub(1)
). The undeformed configuration, for time equal to 0,
is written immediately
file_res = XDMFFile('beam_result.xdmf')
file_res.parameters['functions_share_mesh'] = True
file_res.parameters['rewrite_function_mesh'] = False
file_res.parameters["flush_output"] = True
file_res.write(u.sub(0), t=0.)
file_res.write(u.sub(1), t=0.)
As customary in nonlinear structural analysis the load is divided into num_incs
load
increments.
num_incs = 1200
for i in range(1, num_incs + 1, 1):
print('=====================')
print('Load fraction ', 1. / num_incs * i)
print('=====================')
For each of them the corresponding loads are first computed
distributed_force.assign(Constant(distributed_final_force.values() / num_incs * i))
distributed_moment.assign(Constant(distributed_final_moment.values() / num_incs * i))
concentrated_force.assign(Constant(concentrated_final_force.values() / num_incs * i))
concentrated_moment.assign(Constant(concentrated_final_moment.values() / num_incs * i))
Then, the residual -functional
and the Jacobian matrix J
are assembled together
using the newly defined r_assemble_system
method. The third argument, u.sub(1)
specifies which subspace defines the rotation field
A, b = r_assemble_system(J, -functional, u.sub(1))
After assembly, the boudary conditions can be applied. The nonlinear problem can be solved by iterating until the residual norm is small enough:
[bcs.apply(A, b) for bcs in bc]
iteration = 0
resnorm = b.norm('l2')
print('Iteration = ', iteration, '\tResNorm = ', resnorm)
while resnorm > 1.E-6:
Within each iteration, solve the linear system and put the result into uinc
; after that, increment the solution u
solve(A, uinc.vector(), b)
u.vector()[:] += uinc.vector()
Immediately after computing a new configuration call swe.Switch()
, so that the
rotation vector magnitude is kept under \(\pi\)
swe.Switch()
Assemble the Jacobian matrix and residual vector in the new configuration,
apply the boundary conditions, and compute the residual norm that will
be checked in the while
loop test for convergence
A, b = r_assemble_system(J, -functional, u.sub(1))
[bcs.apply(A, b) for bcs in bc]
resnorm = b.norm('l2')
iteration += 1
print('Iteration = ', iteration, '\tResNorm = ', resnorm)
After convergence of the load step write the solution at time float(i) / num_incs
file_res.write(u.sub(0), t=float(i) / num_incs)
file_res.write(u.sub(1), t=float(i) / num_incs)