1D Linear elasticity with multi-vector force application

variational formulation
1

Hi,

I’m currently attempting to model something that looks similar to this simplified Hyperworks representation.

Zoomedinversion
Zoomedinversion1059×578 22.8 KB

With forces applied orthogonal to the beam and the moments parallel, the outputs are displacement and twist. My current model allows the application of forces in the vertical direction (Z-axis) only, where as to create the desired scenario would require the application of forces in the horizontal axis (X-axis) also in order to represent the resultant, orthogonal force.

I’ll provide the working code here.

import meshio
from matplotlib import pyplot as plt
import time
from dolfin import *

t0 = time.time()

# Read in mesh
msh = meshio.read("Linemesh.msh")
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"line": msh.cells["line"]}))
mesh = Mesh()
filename = "mesh.xdmf"
f = XDMFFile(filename)
f.read(mesh)
k = 1

# Define function space
UE = FiniteElement("CG",mesh.ufl_cell(),k) # Displacement u
VE = FiniteElement("CG",mesh.ufl_cell(),k) # EI*u''
Phi = FunctionSpace(mesh,"CG", k) # Phi
W = FunctionSpace(mesh,MixedElement([UE,VE]))

# This is the most versatile way to use spatial coordinates in a formulation:
x = SpatialCoordinate(mesh)

# Characteristic functions for each end of the beam, for use in boundary
# integrals for weak enforcement of boundary conditions:
leftChar = 10000-x[0]
rightChar = x[0]

# Defining properties
# Outer Diameter
D0 = 30 #mm
D1 = 30 #mm
D2 = 30 #mm
Do = Expression('x[0] <= 250 ? D0 : ((250 < x[0] && (x[0] <= 500)) ? D1 : D2)'\
               ,degree = 0, D0 = D0, D1 = D1, D2 = D2) #mm
               
# Inner Diameter
D0 = 25 #mm
D1 = 25 #mm
D2 = 25 #mm
Di = Expression('x[0] <= 250 ? D0 : ((250 < x[0] && (x[0] <= 500)) ? D1 : D2)'\
               ,degree = 0, D0 = D0, D1 = D1, D2 = D2) #mm

# Second moment of area
I = (3.141592654/64)*((Do**4)-(Di**4)) #mm^4

# Young's Modulus
E0 = 150*10**3 
E1 = 150*10**3 
E2 = 150*10**3 
E = Expression('x[0] <= 250 ? E0 : ((250 < x[0] && (x[0] <= 500)) ? E1 : E2)'\
               ,degree = 0, E0 = E0, E1 = E1, E2 = E2) #Pa   
EI = E*I

# Polar moment of intertia
J = (3.141592654/32)*((Do**4)-(Di**4))#mm^4

# Shear modulus
G0 = 75*10**3
G1 = 75*10**3 
G2 = 75*10**3 
G = Expression('x[0] <= 250 ? G0 : ((250 < x[0] && (x[0] <= 500)) ? G1 : G2)'\
               ,degree = 0, G0 = G0, G1 = G1, G2 = G2) #Pa   
GJ = G*J
qx = Constant(0) #Nm/m Uniformly distributed torque
q = Constant(0) #N/m Uniformly distributed lift

n = FacetNormal(mesh)
h = CellDiameter(mesh)

# PDE residual Euler-Bernoulli bending
def pdeRes(u,v,du,dv):
    return inner(grad(u),grad(du))*dx + inner(v/EI,du)*dx \
        + inner(grad(v),grad(dv))*dx + inner(q,dv)*dx \
        + leftChar*u*inner(n,grad(dv))*ds \
        + rightChar*v*inner(n,grad(du))*ds \
        - leftChar*dot(grad(v),n)*dv*ds \
        - rightChar*dot(grad(u),n)*du*ds

# Set up functions and assemble
u,v= split(TrialFunction(W))
du,dv = split(TestFunction(W))
residual = pdeRes(u,v,du,dv)
a, L = lhs(residual), rhs(residual)
w = Function(W)
A = assemble(a)
b = assemble(L)

# Torsion forumulation
phi = TrialFunction(Phi)
dphi = TestFunction(Phi)
LT = qx*dphi*dx # Linear part
aT = inner(GJ*grad(phi),grad(dphi))*dx #Bilinear part

# Set BCs and solve:
bc = DirichletBC(Phi,Constant(0.0),"x[0]<DOLFIN_EPS")
AT, bT = assemble_system(aT, LT, bc)

# Apply PointSources for force at the endpoint of beam
x = 10000
y = 0
z = 1000
sub1 = W.sub(1)
ptSrc = PointSource(sub1, Point(x, y ,z), 10)
ptSrc.apply(b)

# Apply PointSources for conentrated torque (Nmm)
ptSrcT = PointSource(Phi, Point(6660, 0, 0), 10)
ptSrcT.apply(bT)

# Solve
u,v = w.split()
phi_h = Function(Phi)
solve(A, u.vector() ,b)
solve(AT ,phi_h.vector(),bT)

# Compute magnitude of displacement                    
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('max deflection:',
u_magnitude.vector().get_local().max())

# Compute magnitude of twist
V = FunctionSpace(mesh, 'CG', 2)
u_magnitude = sqrt(dot(phi_h, phi_h))
u_magnitude = project(u_magnitude, V)
print('max twist:',
u_magnitude.vector().get_local().max())

# Write results to file
fileResults = XDMFFile('Results.xdmf')
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True

U = FunctionSpace(mesh, UE)
u = project(u, U)
disp = Function(U, name='Displacement')
disp.assign(u)
fileResults.write(disp,0.)

Twist = Function(V, name='Twist')
Twist.assign(phi_h)
fileResults.write(Twist,0.)

fileResults.close()
t1 = time.time()

total = t1-t0
print(total)

I believe that in order to apply forces in the X direction I need to change
VE = FiniteElement("CG",mesh.ufl_cell(),k)
to
VE = VectorElement('CG', mesh.ufl_cell(), k)

From here I adapted the residual to

# PDE residual - Euler Bernoulli
def pdeRes(u,v,du,dv):
    return inner(grad(u),grad(du))*dx + inner(v[0]/EI,du)*dx \
        + inner(grad(v[0]),grad(dv[0]))*dx + inner(q,dv[0])*dx \
        + inner(grad(u),grad(du))*dx + inner(v[1]/EI,du)*dx \
        + inner(grad(v[1]),grad(dv[1]))*dx + inner(q,dv[1])*dx \
        + inner(grad(u),grad(du))*dx + inner(v[2]/EI,du)*dx \
        + inner(grad(v[2]),grad(dv[2]))*dx + inner(q,dv[2])*dx \

# Define root boundary condition
def root_boundary(x, on_boundary):
    tol = 1
    return on_boundary and x[0] < tol

bc = DirichletBC(W, Constant((0., 0., 0., 0.)), root_boundary)

However, I don’t believe that it is correct. I’m not too familiar with FEA in general and so unsure about the correct formulation. I do however get the same answer for deflection. But when applying the Pointsource to a subspace of VE, which I believe would be the correct way to apply the force at along a different axis, it simply gives me a third of the deflection.

sub1 = W.sub(1)
ptSrc = PointSource(sub1, Point(x, y ,z), 10)

Gives an answer of 1751.17mm deflection.

sub1 = W.sub(1)
ptSrc = PointSource(sub1.sub(0), Point(x, y ,z), 10)

Gives an answer of 583.72mm deflection.

This is incorrect. I’m assuming that my formulation must be incorrect and that simply copying the previous residual for the FiniteElements isn’t the right way to go about it for conversion to VectorElements.

Any help or guidance would be extremely appreciated!

Edit:
I also attempted something like this

# PDE residual - Euler Bernoulli
def pdeRes(u,v,du,dv):
    return inner(grad(u),grad(du))*dx(0) + inner(v[0]/EI,du)*dx(0) \
        + inner(grad(v[0]),grad(dv[0]))*dx(0) + inner(q,dv[0])*dx(0) \
        + inner(grad(u),grad(du))*dx(1) + inner(v[1]/EI,du)*dx(1) \
        + inner(grad(v[1]),grad(dv[1]))*dx(1) + inner(q,dv[1])*dx(1) \
        + inner(grad(u),grad(du))*dx(2) + inner(v[2]/EI,du)*dx(2) \
        + inner(grad(v[2]),grad(dv[2]))*dx(2) + inner(q,dv[2])*dx(2) \

But this is a bit of a guess now and gives me the deflection as ‘nan’.