Hi,
I’m working on a dynamic fracture simulation in legacy DOLFINx, where I’ve implemented the HHT-α time‐integration scheme (see the elastodynamics demo, implemented before here: https://comet-fenics.readthedocs.io/en/latest/demo/elastodynamics/demo_elastodynamics.py.html ) together with a phase‐field fracture model. To extract reaction forces, I’m applying both the method 1 (using the post-processed stress) and method 2 (virtual work) techniques outlined in the https://bleyerj.github.io/comet-fenicsx/tips/computing_reactions/computing_reactions.html.
Despite this, my dynamic code is giving erroneous reaction values for both methods, and I am not sure where is the issue. Below is the snippet for elastodynamics I’m using for the Single Edge Notch problem with Mode-I displacement boundary conditions. Any insights would be greatly appreciated.
# Impact parameters
t0_imp = ScalarType(20e-6)
v0_imp = ScalarType(75)
def update_disp_imp(t: float):
t_val = ScalarType(t)
if t_val < t0_imp:
disp_imp.value = ScalarType(v0_imp * t_val)
else:
disp_imp.value = ScalarType(v0_imp * t0_imp)
disp_imp_bot.value = -disp_imp.value
dofs_rtop_pt = fem.locate_dofs_topological(V.sub(0), 0, rtop_pt)
dofs_bot = fem.locate_dofs_topological(V.sub(1), fdim - 1, bottom_facets)
dofs_top = fem.locate_dofs_topological(V.sub(1), fdim - 1, top_facets)
bclb = fem.dirichletbc(ScalarType(0), dofs_rtop_pt, V.sub(0))
bcl = fem.dirichletbc(disp_imp_bot, dofs_bot, V.sub(1))
bct = fem.dirichletbc(disp_imp, dofs_top, V.sub(1))
bcs_du = [bclb, bcl, bct]
#setting initial conditions
u.x.array[:] = 0.
u.x.scatter_forward()
def m_form(u, v):
return rho_s*ufl.inner(u, v)*dx
def k_form(u, v, z):
return ufl.inner(sigma_d(u, z), ufl.sym(ufl.grad(v)))*dx
def update_a(u, u_prev, v_prev, a_prev, dt_):
b = float(beta)
return (u - u_prev - dt_ * v_prev) / (b * dt_**2) - ((1 - 2 * b) / (2 * b)) * a_prev
def update_v(a, u_prev, v_prev, a_prev, dt_):
g = float(gamma)
return v_prev + dt_ * ((1 - g) * a_prev + g * a)
# Define function for the weighted displacement
def avg(x_old, x_new, alpha):
return ufl.as_vector([alpha * x_old[i] + (1 - alpha) * x_new[i] for i in range(d)])
def update_fields(u, u_prev, v_prev, a_prev):
"""Update fields at the end of each time step."""
u_arr = u.x.array
u_prev_arr = u_prev.x.array
v_prev_arr = v_prev.x.array
a_prev_arr = a_prev.x.array
dt = float(stepsize)
# Compute new acceleration and velocity arrays
a_arr = update_a(u_arr, u_prev_arr, v_prev_arr, a_prev_arr, dt)
v_arr = update_v(a_arr, u_prev_arr, v_prev_arr, a_prev_arr, dt)
# Overwrite the previous‐step arrays in place
a_prev.x.array[:] = a_arr
v_prev.x.array[:] = v_arr
u_prev.x.array[:] = u_arr
# Scatter the updated arrays to synchronize ghost values
a_prev.x.scatter_forward()
v_prev.x.scatter_forward()
u_prev.x.scatter_forward()
# These lines can remain as before
a_new = update_a(du, u_prev, v_prev, a_prev, stepsize)
v_new = update_v(a_new, u_prev, v_prev, a_prev, stepsize)
# ------------------- PDE1: Disp Equation --------------------
Res_u = k_form(avg(u_prev, du, alpha_f), v, z) + m_form(avg(a_prev, a_new, alpha_m), v)
a_form = ufl.lhs(Res_u)
L_form = ufl.rhs(Res_u)
problem_u = fem.petsc.LinearProblem(a_form, L_form, bcs=bcs_du, u=u, petsc_options=petsc_options_lin)
while t-stepsize < T:
update_disp_imp(t)
u_trial.x.array[:] = u.x.array
out = problem_u.solve() #Solving for u
u.x.scatter_forward()
update_fields(u, u_prev, v_prev, a_prev)
u.x.scatter_forward()
# Method 1
# # calculating reaction force by integration of stress over a subdomain
F_sig_top = comm.allreduce(fem.assemble_scalar(fem.form(sigma(u)[1, 1] * ds_m(TOP))), op=MPI.SUM)
F_sig_bot = comm.allreduce(fem.assemble_scalar(fem.form(sigma(u)[1, 1] * ds_m(BOT))), op=MPI.SUM)
# ---------------------------------------------------------------------------------------------------------------
# Method 2
# # For virtual work computation
aform = ufl.inner(sigma(du), epsilon(v)) * dx
residual_vw = ufl.action(aform, u)
# Create test function for virtual work
v_reac = fem.Function(V)
# Top reaction
v_reac.x.array[:] = 0.0
dofs_top = fem.locate_dofs_topological(V.sub(1), fdim - 1, top_facets)
bcRy_top = fem.dirichletbc(ScalarType(1), dofs_top, V.sub(1))
fem.set_bc(v_reac.x.array, [bcRy_top])
vw_form = fem.form(ufl.action(residual_vw, v_reac))
F_action_top = comm.allreduce(fem.assemble_scalar(vw_form), op=MPI.SUM)
# Bottom reaction
v_reac.x.array[:] = 0.0
dofs_bot = fem.locate_dofs_topological(V.sub(1), fdim - 1, bottom_facets)
bcRy_bot = fem.dirichletbc(ScalarType(1.0), dofs_bot, V.sub(1))
fem.set_bc(v_reac.x.array, [bcRy_bot])
vw_form = fem.form(ufl.action(residual_vw, v_reac))
F_action_bot = comm.allreduce(fem.assemble_scalar(vw_form), op=MPI.SUM)