Hi, I’m comming from this tutorial, link, where the reaction forces and bending moments are computed using two methods, One by using the displacements fields and other using virtual work of the forces. I tried to apply this in my case which is a 3D time dependent problem, but seems that I’m missing something here because here both methods differs too much, and also the virtual work give the same values for all directions.
I prepare this working example, to clariffy this, based on another tutorial to keep it simple:
import numpy as np
from mpi4py import MPI
import ufl
from dolfinx import fem, mesh
import dolfinx.fem.petsc
from dolfinx.mesh import create_box, CellType
L = 8.0
H = 0.2
B = 0.1
domain = create_box(
MPI.COMM_WORLD,
[[0.0, -B / 2, -H / 2], [L, B / 2, H / 2]],
[8, 2, 2],
CellType.hexahedron,
)
dim = domain.topology.dim
dx = ufl.Measure("dx", domain=domain)
degree = 2
shape = (dim,)
V = fem.functionspace(domain, ("Q", degree, shape))
u = fem.Function(V, name="Displacement")
def left(x):
return np.isclose(x[0], 0.0)
fdim = domain.topology.dim - 1
marked_values = []
marked_facets = []
# Concatenate and sort the arrays based on facet indices
facets = mesh.locate_entities_boundary(domain, fdim, left)
marked_facets.append(facets)
marked_values.append(np.full_like(facets, 1))
marked_facets = np.hstack(marked_facets)
marked_values = np.hstack(marked_values)
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tag)
clamped_dofs = fem.locate_dofs_geometrical(V, left)
u_bc = fem.Function(V)
bcs = [fem.dirichletbc(u_bc, clamped_dofs)]
E = fem.Constant(domain, 210e3)
nu = fem.Constant(domain, 0.3)
rho = fem.Constant(domain, 7.8e-3)
f = fem.Constant(domain, (0.0,) * dim)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2 / (1 + nu)
def epsilon(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return lmbda * ufl.tr(epsilon(v)) * ufl.Identity(dim) + 2 * mu * epsilon(v)
u_old = fem.Function(V)
v_old = fem.Function(V)
a_old = fem.Function(V)
a_new = fem.Function(V)
v_new = fem.Function(V)
beta_ = 0.25
beta = fem.Constant(domain, beta_)
gamma_ = 0.5
gamma = fem.Constant(domain, gamma_)
dt = fem.Constant(domain, 0.0)
a = 1 / beta / dt**2 * (u - u_old - dt * v_old) + a_old * (1 - 1 / 2 / beta)
a_expr = fem.Expression(a, V.element.interpolation_points())
v = v_old + dt * ((1 - gamma) * a_old + gamma * a)
v_expr = fem.Expression(v, V.element.interpolation_points())
u_ = ufl.TestFunction(V)
du = ufl.TrialFunction(V)
def mass(u, u_):
return rho * ufl.dot(u, u_) * ufl.dx
def stiffness(u, u_):
return ufl.inner(sigma(u), epsilon(u_)) * ufl.dx
Residual = mass(a, u_) + stiffness(u, u_) - ufl.dot(f, u_) * ufl.dx
Residual_du = ufl.replace(Residual, {u: du})
a_form = ufl.lhs(Residual_du)
L_form = ufl.rhs(Residual_du)
problem = fem.petsc.LinearProblem(
a_form, L_form, u=u, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
# ---------------------------- #
# for virtual work computation #
# ---------------------------- #
residual_vw = ufl.action(a_form, u) - L_form
v_reac = fem.Function(V)
virtual_work_form = fem.form(ufl.action(residual_vw, v_reac))
x = ufl.SpatialCoordinate(domain)
def one(x):
values = np.zeros((1, x.shape[1]))
values[0] = 1.0
return values
def y(x):
values = np.zeros((1, x.shape[1]))
values[0] = x[1]
return values
# ---------- #
# solve loop #
# ---------- #
t = 0
Nsteps = 400
times = np.linspace(0, 2, Nsteps + 1)
for i, dti in enumerate(np.diff(times)):
dt.value = dti
t += dti
if t <= 0.2:
f.value = np.array([0.0, 1.0, 1.5]) * t / 0.2
else:
f.value *= 0.0
problem.solve()
u.x.scatter_forward()
v_new.interpolate(v_expr)
a_new.interpolate(a_expr)
u.vector.copy(u_old.vector)
v_new.vector.copy(v_old.vector)
a_new.vector.copy(a_old.vector)
Rx = fem.assemble_scalar(fem.form(-sigma(u)[0, 0] * ds))
Ry = fem.assemble_scalar(fem.form(-sigma(u)[0, 1] * ds))
Rz = fem.assemble_scalar(fem.form(-sigma(u)[0, 2] * ds))
Mx = fem.assemble_scalar(fem.form((x[1] * sigma(u)[0, 2] - x[2] * sigma(u)[0, 1]) * ds))
My = fem.assemble_scalar(fem.form((x[2] * sigma(u)[0, 0] - x[0] * sigma(u)[0, 2]) * ds))
Mz = fem.assemble_scalar(fem.form((x[0] * sigma(u)[0, 1] - x[1] * sigma(u)[0, 0]) * ds))
# virtual work
u_bc.sub(0).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rx_vw = fem.assemble_scalar(virtual_work_form)
u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(1).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Ry_vw = fem.assemble_scalar(virtual_work_form)
u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(2).interpolate(one)
fem.set_bc(v_reac.vector, bcs)
Rz_vw = fem.assemble_scalar(virtual_work_form)
u_bc.vector.set(0.0)
v_reac.vector.set(0.0)
u_bc.sub(0).interpolate(y)
fem.set_bc(v_reac.vector, bcs)
Mz_vw = fem.assemble_scalar(virtual_work_form)
print(f"Reaction Rx = {Rx:.2f} / {Rx_vw:.2f}")
print(f"Reaction Ry = {Ry:.2f} / {Ry_vw:.2f}")
print(f"Reaction Rz = {Rz:.2f} / {Rz_vw:.2f}")
print(f"Moment Mx = {Mx:.2f}")
print(f"Moment My = {My:.2f}")
print(f"Moment Mz = {Mz:.2f} / {Mz_vw:.2f}")
print("-" * 45)
which returns something like this:
Reaction Rx = 0.00 / 0.00
Reaction Ry = -0.01 / -0.00
Reaction Rz = -0.02 / -0.00
Moment Mx = -0.00
Moment My = -0.00
Moment Mz = 0.00 / -0.00
---------------------------------------------
Reaction Rx = -0.00 / 12.15
Reaction Ry = 108.17 / 12.15
Reaction Rz = 633.58 / 12.15
Moment Mx = 0.00
Moment My = 43.60
Moment Mz = 318.33 / 10.36
---------------------------------------------
Reaction Rx = -0.00 / -39.31
Reaction Ry = -426.01 / -39.31
Reaction Rz = -824.05 / -39.31
Moment Mx = 0.00
Moment My = -94.04
Moment Mz = 323.53 / -26.82
---------------------------------------------
Reaction Rx = -0.00 / 44.44
Reaction Ry = 475.62 / 44.44
Reaction Rz = 903.02 / 44.44
Moment Mx = 0.00
Moment My = -246.30
Moment Mz = 299.27 / 33.03
...
As you can see the results above are two diferents, and from my point of view, the wrong is the virtual work method here.