How to add torsion to a beam problem?

1

Hi,

So, I’ve managed to script a simple 2D cantilever beam under several point loads and now wish to add torque (twisting moments) at several points also. From this I’d like to output the maximum rotation at the end of the beam.

How should I go about doing this? Any guidance/help will be much appreciated.

Cheers.

2

Hi,

I am also doing the same thing as you, I am wondering if you already got a solution to your question? I’d appreciate it if you can share it. Thank you!

T.Z.

3

Are you looking to do this in legacy DOLFIN or DOLFINx? Do you have a beam model that already captures torsion? If not, I’d suggest looking at this demo. If so, what does you formulation look like? (Mixed? Hermetian?)

Applying point moments is typically handled mathematically through a dirac delta applied a specific location. There are two sorts of cases that are relevant for us:

  1. the point moment is applied at a node. In this case, we can assemble our system and directly modify the RHS PETSC vector.
  2. the point moment is not applied at a node. In this case, we need to use some sort of PointSource functionality. In legacy DOLFIN, we can directly use PointSource(). In DOLFINx, PointSource() has not yet been reimplemented, so in this case, I’m not sure how to implement that yet ¯\_(ツ)_/¯
4

I’m using dolfinx, and I am trying to bring a single torsion to my mesh generated by using GMSH(a 3D beam) or even just to a minimal 3D beam mesh with rectangular section. The demo you posted here maybe not fits my situation. To 2: Yes it’s a pity we don’t have PointSource in dolfinx. To 1: maybe could you please explain it a little bit more with a minimal example? I’d appreciate it very much.

5

If the cross-section is symmetric, isotropic and homogenous, then the most complicated part about adding twisting into your problem is computing a torsional constant (often denoted “J”). Outside of a circle, there are no exact closed form solutions for a torsional constant. Otherwise, it is much more complicated, as there are various couplings between the different deformations (bending-twisting, bending-bending, etc).

Computing the torsional constant J involves solving a pure neumann Laplace’s equation. To solve pure neumann problems like this, we have to either use Lagrange multipliers (which are currently problematic in dolfinx ) or we can remove a constant nullspace.

An alternative (and more approachable) formulation involves using a Prandtl stress function. A MWE for this approach is:

from dolfinx import mesh
from dolfinx.fem import (FunctionSpace,Constant,Function,form,dirichletbc,
                         locate_dofs_geometrical,assemble_scalar)
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from ufl import TrialFunction,TestFunction,inner,grad,dx,ds,SpatialCoordinate
from petsc4py.PETSc import ScalarType
import numpy as np
'''
Computing the torsional constant of a square using a Prandtl Stress function 

'''

# Create mesh
square = mesh.create_unit_square(MPI.COMM_WORLD, 20,20, mesh.CellType.quadrilateral)
domain = square
V = FunctionSpace(domain, ("CG", 1))

phi = TrialFunction(V)
v = TestFunction(V)

x = SpatialCoordinate(domain)
#G and Theta are virutal loads that can be anything within numerical reason
# their purpose is to perturb the solution away from a "boring" LaPlace equation
# with the boundary equal to 0 (which would just be 0 everywhere)
G = 10
Theta = 0.1
f = Constant(domain, ScalarType(2*G*Theta))
a = (inner(grad(phi), grad(v)) )*dx
L = f*v*dx 

def left(x):
    return np.isclose(x[0], 0)
def right(x):
    return np.isclose(x[0],1)
def top(x):
    return np.isclose(x[1],0)
def bot(x):
    return np.isclose(x[1],1)

bc1 = dirichletbc(ScalarType(0), locate_dofs_geometrical(V, left),V)
bc2 = dirichletbc(ScalarType(0), locate_dofs_geometrical(V, right),V)
bc3 = dirichletbc(ScalarType(0), locate_dofs_geometrical(V, top),V)
bc4 = dirichletbc(ScalarType(0), locate_dofs_geometrical(V, bot),V)

# Compute solution
phih = Function(V)
problem = LinearProblem(a,L,bcs=[bc1,bc2,bc3,bc4])
phih = problem.solve()

It = (2/(G*Theta))*assemble_scalar(form(phih*dx))
print("torsional constant:")
print(It)
1 Like