Hi
I was trying out modifying hypereslaticity demo for dynamic simulation. What I did:
- I don’t apply volumetric force and boundary traction, i.e.
E_strain = psi * dx - I used mixed element for simulate displacement (u) and velocity (v).
Physical problem: with a Neo-hook cube, I fixed one side (zero displacement) and pull the other side along x-axis under a fixed speed. With that said, I applied DirichletBC for uy, uz, vy, vz = 0, while ux = 0.1*t and vx = 0.1
Observation (result at t=1 attached, original and zoomed ): if we look at one surface parallel to x-axis, there are surfaces that displacement is not symmetric about the middle of that surface.
Here’s a piece of code how I modified the demo
U, V = VectorElement('CG', mesh.ufl_cell(), 1), VectorElement('CG', mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, MixedElement([U, V]))
test_func = TestFunction(ME)
tu, tv = split(test_func) # test function components
solution, solution_old = Function(ME), Function(ME)
u, v = split(solution) # solution components
u0, v0 = split(solution_old) # solution components from previous time step
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 # psi from the hyperelastic demo
Pi = psi*dx
residual = inv_dt * rho * inner(tv, (v - v0)) * dx +\ #inv_dt = 1/dt
derivative(Pi, u, tv) +\
inv_dt * inner(tu, (u - u0)) * dx - inner(tu, v) * dx
u_right = Expression(('0.1 * t', '0', '0'), t=0, degree = 2)
bcul = DirichletBC(u_space, Constant((0., 0., 0.)), left)
bcur = DirichletBC(u_space, u_right, right)
bcvl = DirichletBC(v_space, Constant((0., 0., 0.)), left)
bcvr = DirichletBC(v_space, Constant((0.1, 0., 0.)), right)
bcs = [bcul, bcur, bcvl, bcvr]
# dt, t setup ...
while t < t_total:
t += dt
u_right = t
solver.solve ...
solution_old.assign(solution)
I don’t know what I missed. Supposedly as I uniformly apply Dirichlet BC on side x[0]=1, on each boundary surface that is parallel to x-axis, the displacement should be symmetric about the middle of that surface, but in my result, look at the dark line in the red region, apparently it’s not what I expected.
I thought it was because of I used derivative, so I hard-coded pk1 stress of the Neohook but the results are the same. Another observation is that, I don’t see this weird result in static solid simulation. I am not sure if it’s because I mis-used mixed element
Any help is appreciated.
Thanks
Victor