Hi all,
I was wondering if there is any examples on modelling with axisymmetric formulation that has an out-of-plane twisting (for example, torsion of a shaft). I found this tutorial very helpful but it doesn’t have any out of plane terms for the displacement.
Thank you in advance!
You can use a 1D torsional model if you compute a “torsional constant” following this link. The 1D torsional model is quite simple (only a function of rotation of the cross section).
Additionally, a 3D implementation exists here at this link
Thank you! While that’s not really what I wanted to implement as in the end I would like to model the torsion of a composite with subdomains, which I don’t know the torsional constant a priori.
To make my question clearer. Here is a mwe that models the problem (but does not gives the correct value).
In this example, I am trying to model the simple torsion of a cylinder with axisymmetric formulation. The height of the cylinder is 1 and the radius is also 1. I apply the torsion angle per unit length is \gamma = 0.001. All boundaries remains flat during the deformaiton. The material is neo-Hookean with \lambda = \mu = 1.
from dolfinx import fem, mesh, nls, log, io
from mpi4py import MPI
import numpy as np
import ufl
from petsc4py import PETSc
from dolfinx.nls.petsc import NewtonSolver
# Material parameters
mu = 1.0
lambda_ = 1.0
# Create the domain
msh = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
# Create a function space
u_el = ufl.VectorElement("CG", msh.ufl_cell(), 2)
u_theta_el = ufl.FiniteElement("CG", msh.ufl_cell(), 2)
state_space = fem.FunctionSpace(msh, u_el * u_theta_el)
u_space, _ = state_space.sub(0).collapse()
theta_space, _ = state_space.sub(1).collapse()
ux_space, _ = u_space.sub(0).collapse()
uy_space, _ = u_space.sub(1).collapse()
def left(x):
return np.isclose(x[0], 0.0)
def right(x):
return np.isclose(x[0], 1.0)
def bottom(x):
return np.isclose(x[1], 0.0)
def top(x):
return np.isclose(x[1], 1.0)
# Define Dirichlet boundary conditions
fdim = msh.topology.dim - 1
left_facets = mesh.locate_entities_boundary(msh, fdim, left)
right_facets = mesh.locate_entities_boundary(msh, fdim, right)
bottom_facets = mesh.locate_entities_boundary(msh, fdim, bottom)
top_facets = mesh.locate_entities_boundary(msh, fdim, top)
dofs_left_ux = fem.locate_dofs_topological(state_space.sub(0).sub(0), fdim, left_facets)
dofs_bottom_uy = fem.locate_dofs_topological(state_space.sub(0).sub(1), fdim, bottom_facets)
dofs_right_ux = fem.locate_dofs_topological(state_space.sub(0).sub(0), fdim, right_facets)
dofs_top_uy = fem.locate_dofs_topological(state_space.sub(0).sub(1), fdim, top_facets)
dofs_left_utheta = fem.locate_dofs_topological(state_space.sub(1), fdim, left_facets)
dofs_bottom_utheta = fem.locate_dofs_topological(state_space.sub(1), fdim, bottom_facets)
dofs_right_utheta = fem.locate_dofs_topological((state_space.sub(1), theta_space), fdim, right_facets)
bc_left_ux = fem.dirichletbc(PETSc.ScalarType(0), dofs_left_ux, state_space.sub(0).sub(0))
bc_bottom_uy = fem.dirichletbc(PETSc.ScalarType(0), dofs_bottom_uy, state_space.sub(0).sub(1))
bc_right_ux = fem.dirichletbc(PETSc.ScalarType(0), dofs_right_ux, state_space.sub(0).sub(0))
bc_top_uy = fem.dirichletbc(PETSc.ScalarType(0), dofs_top_uy, state_space.sub(0).sub(1))
bc_left_utheta = fem.dirichletbc(PETSc.ScalarType(0), dofs_left_utheta, state_space.sub(1))
bc_bottom_utheta = fem.dirichletbc(PETSc.ScalarType(0), dofs_bottom_utheta, state_space.sub(1))
# Applied load
gamma = 0.001
u_theta_right_expr = lambda x: gamma * x[1] / x[0]
u_theta_right = fem.Function(theta_space)
u_theta_right.interpolate(u_theta_right_expr)
bc_right_utheta = fem.dirichletbc(u_theta_right, dofs_right_utheta, state_space.sub(1))
bcs = [bc_left_ux, bc_bottom_uy, bc_right_ux, bc_top_uy, bc_left_utheta, bc_bottom_utheta, bc_right_utheta]
# Define the gradient in cylindrical coordinate
x = ufl.SpatialCoordinate(msh)
def cyl_Grad(u, utheta):
Gradu = ufl.as_tensor([[u[0].dx(0), -utheta / x[0], u[0].dx(1)],
[utheta.dx(0), u[0]/x[0], utheta.dx(1)],
[u[1].dx(0), 0, u[1].dx(1)]])
return Gradu
# Define the volume and boundary measures
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=msh, metadata=metadata)
ds = ufl.Measure("ds", domain=msh, metadata=metadata)
state = fem.Function(state_space)
u, utheta = ufl.split(state)
# Kinematics
F = ufl.variable(cyl_Grad(u, utheta) + ufl.Identity(3))
C = F.T * F
Ic = ufl.tr(C)
J = ufl.det(F)
psi_strain = 0.5 * mu * (Ic - 3 - 2 * ufl.ln(J)) + 0.5 * lambda_ * (J - 1) ** 2
U_strain = psi_strain * x[0] * dx
# Define the weak form
R = ufl.derivative(U_strain, state)
problem = fem.petsc.NonlinearProblem(R, state, bcs=bcs)
solver = NewtonSolver(msh.comm, problem)
solver.atol = 1e-10
solver.rtol = 1e-10
solver.convergence_criterion = "incremental"
log.set_log_level(log.LogLevel.INFO)
num_its, converged = solver.solve(state)
assert(converged)
u_sol, utheta_sol = state.split()
u_sol.name = "u"
utheta_sol.name = "u_theta"
V_u_P1 = fem.functionspace(msh, ("CG", 1, (2,)))
u_P1 = fem.Function(V_u_P1)
u_P1.interpolate(u_sol)
u_P1.name = "u_P1"
V_utheta_P1 = fem.functionspace(msh, ("CG", 1))
utheta_P1 = fem.Function(V_utheta_P1)
utheta_P1.interpolate(utheta_sol)
utheta_P1.name = "u_theta_P1"
with io.XDMFFile(MPI.COMM_WORLD, "results.xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_function(u_P1)
xdmf.write_function(utheta_P1)
Here are the comparison between Abaqus and FEniCSx. Does anyone have any idea about where this is wrong?
One simple question because it is a bit hard for me to tell from the above pictures… Are these figures made with mesh/mesh sizes? I can’t exactly tell, but it looks like the ABAQUS one may be a slight more refined mesh?
