Twist of two connected bars with circular cross-sections

[adtzlr]

Twist of two connected bars with circular cross-sections

• ✅ Mixed-Field Formulation (Displacement, Pressure, Volume Ratio)
• ✅ Nearly-Incompressible Neo-Hookean Material Formulation
• ✅ Q2-Q1 Hexahedron Element Formulation (Continuous tri-quadratic displacement and (dis-) continuous linear pressure/volume-ratio interpolations)
  

import felupe as fem
import numpy as np

from pypardiso import spsolve

# create a cube
mesh = fem.mesh.Cube(a=(-1, -1, 0), b=(1, 1, 20), n=(3, 3, 11))
mesh = fem.mesh.convert(mesh, 2, True, True, True)

# https://www.xarg.org/2017/07/how-to-map-a-square-to-a-circle/
mesh.points[:, :2] *= np.sqrt(1 - mesh.points[:, [1, 0]] **2 / 2)

msh2 = mesh.copy()
msh2.points[:, 1] -= 2

mesh = fem.mesh.concatenate([mesh, msh2])
mesh = fem.mesh.sweep(mesh)
mesh.save()

import felupe as fem
import numpy as np

# create a hexahedron-region on a cube
region = fem.RegionTriQuadraticHexahedron(mesh)

# add a mixed field container (with displacement, pressure and volume ratio)
mesh_dual = mesh.copy()
mesh_dual.update(mesh.cells[:, :8])
mesh_dual.cell_type = "hexahedron"
mesh_dual.disconnect()
field = fem.FieldsMixed(region, n=3, mesh=mesh_dual)

# apply a uniaxial elongation on the cube
boundaries = dict(
    bottom=fem.Boundary(field[0], fz=mesh.points[:, 2].min()),
    top=fem.Boundary(field[0], fz=mesh.points[:, 2].max()),
)

# create the values for the rotation
move_top = []
for angle_deg in fem.math.linsteps([0, 360, 360], num=30):
    
    # rotate point coordinates on top boundary
    top = mesh.points[boundaries["top"].points]
    top_rotated = fem.mesh.rotate(
        points=top, 
        cells=None, 
        cell_type=None,
        angle_deg=angle_deg, 
        axis=2, 
        center=[0, 0, 0]
    )[0]

    move_top.append((top_rotated - top).ravel())

# create the values for the rotation
move_bottom = []
for angle_deg in fem.math.linsteps([0, 0, 360], num=30):
    
    # rotate point coordinates on top boundary
    bottom = mesh.points[boundaries["bottom"].points]
    bottom_rotated = fem.mesh.rotate(
        points=bottom, 
        cells=None, 
        cell_type=None,
        angle_deg=angle_deg, 
        axis=2, 
        center=[0, 0, 0]
    )[0]

    move_bottom.append((bottom_rotated - bottom).ravel())

# define the constitutive material behaviour and create a solid body
umat = fem.ThreeFieldVariation(
    fem.NeoHooke(mu=1, bulk=5000, parallel=True),
    parallel=True
)
solid = fem.SolidBody(umat, field)

# prepare a step with substeps
step = fem.Step(
    items=[solid], 
    ramp={
        boundaries["bottom"]: move_bottom,
        boundaries["top"]: move_top,
    }, 
    boundaries=boundaries
)

# add the step to a job, evaluate all substeps and create a plot
job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["top"])
job.evaluate(
    filename="result.xdmf", 
    solver=spsolve, 
    tol=1e-4, 
    kwargs={"parallel": True},
)

The code was adopted from here:
Originally posted by @adtzlr in #309 (comment)

[adtzlr]

In Paraview it is necessary to combine two Warp by Vector Filters - one with surface and another one with wireframe representation. This avoids the scattered outline representation in combination with subdivision (nonlinear surface representation) turned on.