# How to add torsion to a beam problem?

**URL:** <https://fenicsproject.discourse.group/t/how-to-add-torsion-to-a-beam-problem/1681>\
**Category:** Uncategorized\
**Created:** [October 7, 2019, 12:20pm UTC](https://fenicsproject.discourse.group/t/how-to-add-torsion-to-a-beam-problem/1681 "2019-10-07T12:20:29Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![jkrokowski](https://avatars.discourse-cdn.com/v4/letter/j/0ea827/32.png) [@jkrokowski](https://fenicsproject.discourse.group/u/jkrokowski)\
**Post date:** [April 3, 2023, 9:32pm UTC](https://fenicsproject.discourse.group/t/how-to-add-torsion-to-a-beam-problem/1681/5 "2023-04-03T21:32:14Z")

</div>

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](https://en.wikipedia.org/wiki/Torsion_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](https://fenicsproject.discourse.group/t/mixed-function-space-with-dolfinx-how-to-implement-an-integral-constraint/9672/2) ) or we can [remove a constant nullspace](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468).

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

```auto
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)

```

---

_[View the full topic](https://fenicsproject.discourse.group/t/how-to-add-torsion-to-a-beam-problem/1681)._