Hello everyone,
I am trying to recreate the fenics example of plasticity (code and output attached). The NGSolve example uses the potential approach but I personally wanted to try the return mapping algorithm for trial stress shown by fenics. The fenics examples combines the elastic and plastic calculations in one step (such that there are no plastic corrections to the stress in the elastic region) and uses the concept of projection on quadrature spaces for computations of stresses and plastic strains with a custom integration measure for the quadrature space.
In the fenics example, they get a different initial norm for the right hand side compared to my version and after one successful solution of the system (enough for the initial elastic steps), the displacement increment that they get is totally different because they residuum norm is different and thus with the displacement increment they get, the external load and the internal strain energy are balanced. My initial external load gets a very different norm (as a Linear Form) compared to their external load and the displacement increment I get does not balance out the external work and internal strain energy thus my residual norm does not go to zero.
The weird thing is that my disp_inc is correct (albeit for the wrong external load) because substracting the external work and the internal energy (even though they are different) still gives me zero solution for the next iteration as it should because the internal energy is subtracted from the external work on the right hand side and thus the linear system gives a zero solution (as it should).
I have seen that there are some DOFs at which there is a discrepancy between external load.vec and internal energy.vec but this shouldn’t be so but then why is the solution zero?
I have attached python scripts for my example (Small_Strain_Plasticity_2D.py) and the fenics script (Fenics_plasticity.py) as well as outputs from both codes as well.
Please let me know whether there is an error in the computation of my external load vector and if so why?
Best regards
Khuldoon Usman
I want to add files but I cannot, gives me error " Unable to get properties from the image"
Hi, yes sorry we have problems with attachments in the forum right now due to kunena/joomla problems… Can you upload the files somewhere else (ie. a github gist) and give a link to it?
Hello Dr. Lehrenfeld,
I have uploaded all the files in the following repo:
The fenics’ files are the code (fenics_plasticity.py) and its corresponding output (fenics_plasticity_out.txt).
My own version in NG solve also has a code (Small_Strain_Plasticity_2D.py) and the corresponding output (Small_Strain_Plasticity_2D_out.txt)
One more thing to add is that the discrepancy between the internal strain energy and the external work purely arises from the DirichletBCs. If you define two linear forms, one for computing the strain energy with the known stress and the other for computing the external work with the known force. Then, the norms of these two forms is different due to the dirichlet BCs imposed when solving for the displacement increment from which stress is then calculated.
Moreover, further refining the mesh with finer elements not only gives a different solution norm and residual norm but even the neumann boundary load computed has a different norm
I’ll probably find some time to look at this issue today as well and keep you posted…
Best,
Matthias
So, after having a first look: in the NGSolve version, “source” does not have boundary conditions applied, i.e. you still have reaction forces corresponding to nodes with Dirichlet boundary conditions.
You can kick these out with
nRes = Norm(bc_projector * source.vec)
print("Norm res:", nRes)
where is defined upfront as
bc_projector = Projector(V.FreeDofs(), True)
This will give you vanishing residual in the elastic regime (i.e. Newton converges after first step).
In the inelastic regime, there seems to be an issue with “Jacobian” stress as Newton does not converge.
I did not spot the bug though.
Besides that, I found it unusual, that you have used (u…trial func, v… tst func)
InnerProduct(Strain(u), Jacobian_stress(Strain(v)))
This is consistent with your FEniCS implementation, but I would have guessed that you need
InnerProduct(Strain(v), Jacobian_stress(Strain(u)))
instead. In my test, it didn’t make a difference though, maybe shadowed by “the” bug?
(didn’t run FEniCS yet to see whether it makes a difference there)
Dear Matthias,
Thank you for the source Norm clarification.
The bug is clearly resulting in different solutions of the system at different mesh sizes even in the elastic regime.
The reason I swapped the test and trial function is because to my understanding, the “a_newton” is symetric even with the “jacobian_stress”.
I look forward to anymore info.
Best regards
Khuldoon