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=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.