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)