Hi all,
I have a mixed element space as follows:
V_el = ufl.VectorElement("CG", cell, 1)
S_el = ufl.FiniteElement("DG", cell, 0)
state_space = fem.FunctionSpace(msh, ufl.MixedElement([V_el, S_el]))
state = fem.Function(state_space, name="state")
u, dp = ufl.split(state)
And an internal variable of plastic strain
Ts = fem.TensorFunctionSpace(msh, ('DG', 0))
plastic_strain_old = fem.Function(Ts)
My script continues as follows, bearing in mind these are the only references to plastic_strain_old until solving:
#set up boundary conditions and material parameters etc
#====================================================
def eps(u):
return ufl.sym(ufl.grad(u))
def sigma(eps_el):
return lam*ufl.tr(eps_el)*ufl.Identity(gdim) + 2*mu*eps_el
theta = np.pi/2 # theta is a known constant for now, although in future it will be a variable added to the state
N = ufl.as_tensor([[ufl.sin(theta),ufl.cos(theta)],
[ufl.cos(theta),-ufl,sin(theta)]])
def strain_energy(total_strain, dp):
eps_el = total_strain-plastic_strain_old - dp*N
return 0.5*(ufl.inner(eps_el, sigma(eps_el)))
def dissipation_potential(dp):
return dp/dt*sy
J_p = strain_energy(eps(u),dp)*dx + dt*dissipation_potential(dp)*dx # Energy functional
R = ufl.derivative(J_p, state)
dR = ufl.derivative(R, state)
#set up a solver etc
I then enter a for loop increasing the displacement applied via a dirichlet bc and solve for the state which minimises the functional, as is standard. Until this point I was under the impression that plastic_strain_old was unaffected by the solver, it was an accumulation variable which should not change until I manually update it later.
However, if I print the plastic_strain_old as follows:
print(f'plastic_strain_old = {l2_projection(plastic_strain_old, Ts).x.array}')
solver_tao.solve(state.vector) # solve the problem
print(f'plastic_strain_old = {l2_projection(plastic_strain_old, Ts).x.array}')
#store solution and update plastic_strain_old
u_sol,dp_sol = state.split()
print(f'N*dp = {l2_projection(N*dp_sol, Ts).x.array}')
plastic_strain_old += N*dp_sol
The sample should yield between steps 7 and 8, but the output looks like this:
=================================================
Step 6 of 10
plastic_strain_old = [0. 0. 0. 0.]
iter = 0, Function value: 0.05, Residual: 4.66224 Constraint: 0.
iter = 1, Function value: 0.0730769, Residual: 10.8786 Constraint: 0.
iter = 2, Function value: 0.072, Residual: 5.07343e-15 Constraint: 0.
plastic_strain_old = [0. 0. 0. 0.]
N*dp = [0. 0. 0. 0.]
=================================================
Step 7 of 10
plastic_strain_old = [0. 0. 0. 0.]
iter = 0, Function value: 0.072, Residual: 4.66224 Constraint: 0.
iter = 1, Function value: 0.098, Residual: 5.72554e-11 Constraint: 0.
plastic_strain_old = [ 4.93159389e-15 0.00000000e+00 0.00000000e+00 -4.93159389e-15]
N*dp = [ 7.04513413e-16 0.00000000e+00 0.00000000e+00 -7.04513413e-16]
=================================================
Step 8 of 10
plastic_strain_old = [ 5.63610731e-15 0.00000000e+00 0.00000000e+00 -5.63610731e-15]
iter = 0, Function value: 0.098, Residual: 16.0755 Constraint: 0.
iter = 1, Function value: 0.126, Residual: 4.11485e-10 Constraint: 0.
plastic_strain_old = [ 0.0016 0. 0. -0.0016]
N*dp = [ 0.0002 0. 0. -0.0002]
Step 9 of 10
plastic_strain_old = [ 0.0018 0. 0. -0.0018]
iter = 0, Function value: 0.126, Residual: 16.0755 Constraint: 0.
iter = 1, Function value: 0.154, Residual: 4.05623e-10 Constraint: 0.
plastic_strain_old = [ 0.0036 0. 0. -0.0036]
N*dp = [ 0.0004 0. 0. -0.0004]
What I expect to see before and after solving at step 8 is plastic_strain_old being zero in both, as before it has not yielded and is zero, and directly after it is zero, as it has not been updated yet.
Somehow during step 8 plastic_strain_old is being updated to the exact value of the total_strain, leading to zero elastic_strain. I don’t understand how plastic_strain_old is being updated during solving, as it should not be affected by the state, should be purely an internal variable that is only updated after the incremental solution is found and a constant whilst the problem is being solved.
Can anyone help please.
Kind regards,
Magnus