How to add torsion to a beam problem?

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