Obtaining strain tensor and Newton-Raphson solver

[HoseongJeong]

Hello everyone!

I am currently working on solving a nonlinear FEM problem that includes material nonlinearity.

Which solver should I use in scikit-fem for this type of problem? (like Newton-raphson solver).

Additionally, I need some guidance on deriving the strain tensor.

• I tried checking the strain tensor by including "print(sym_grad(u))" in the Bilinear form. However, I noticed that the tensor didn't change when I altered the boundary condition from u^1=-0.5 to u^2=-0.5....
  For instance, I expected the shape of the "sym_grad(u)" matrix to be (3,3,number of elements,number of quadrature points) for a 3D element. Thus, I assumed the data at (:,:,0,0) would represent the strain tensor at element 1 and quadrature point 1.
  When i apply u^1=-0.5, and get data using "sym_grad(u)(:,:,0,0)", the result was as follows
  [[0. 0. 0.31100423]
  [0. 0. 0.31100423]
  [0.31100423 0.31100423 0.62200847]]
  and this result didn't change when I apply u^2=-0.5.

• Furthermore, when defining the strain tensor as:
  [e11, e12. e13]
  [e21, e22, e23]
  [e31, e32, e33].
  When I applied axial diaplacement at the boundary,I expected e11 to be significantly higher than the other components when applying axial displacement at the boundary, but this was not the case.

• Additionally, when I defined the boundary condition with a very small displacement, the values of "(sym_grad(u)" were much higher than I anticipated.

Could anyone explain how to correctly derive the strain tensor?

Thank you to everyone in the Scikit-FEM community!
--------------------------------------Printed result-------------------------------
Hey
[[[[0. 0. 0. 0. 0. 0. 0. 0. ]]
[[0. 0. 0. 0. 0. 0. 0. 0. ]]
[[0.31100423 0.08333333 0.31100423 0.08333333 0.08333333 0.0223291 0.08333333 0.0223291 ]]]

[[[0. 0. 0. 0. 0. 0. 0. 0. ]]
[[0. 0. 0. 0. 0. 0. 0. 0. ]]
[[0.31100423 0.31100423 0.08333333 0.08333333 0.08333333 0.08333333 0.0223291 0.0223291 ]]]

[[[0.31100423 0.08333333 0.31100423 0.08333333 0.08333333 0.0223291 0.08333333 0.0223291 ]]
[[0.31100423 0.31100423 0.08333333 0.08333333 0.08333333 0.08333333 0.0223291 0.0223291 ]]
[[0.62200847 0.16666667 0.16666667 0.0446582 0.62200847 0.16666667 0.16666667 0.0446582 ]]]]

--------------------------------------My code-------------------------------

-- coding: utf-8 --

"""
Created on Mon Aug 12 13:50:00 2024

@author: Administrator
"""

r"""Linear elasticity.

This example solves the linear elasticity problem using trilinear elements.

"""
import numpy as np
from skfem import *
from skfem.helpers import ddot, sym_grad, eye, trace
from skfem.models.elasticity import lame_parameters

m = MeshHex().refined(0).with_defaults()
e = ElementVector(ElementHex1())
basis = Basis(m, e, intorder=3)

calculate Lamé parameters from Young's modulus and Poisson ratio

lam, mu = lame_parameters(1e3, 0.3)

def C(T):
return 2. * mu * T + lam * eye(trace(T), T.shape[0])

@BilinearForm
def stiffness(u, v, w):
print("Hey")
### print(sym_grad(u))
return ddot(C(sym_grad(u)), sym_grad(v))

K = stiffness.assemble(basis)

u = basis.zeros()

u[basis.get_dofs('right').nodal['u^3']] = -0.5

u = solve(*condense(K, x=u, D=basis.get_dofs({'left', 'right'})))

sf = 1.0
m = m.translated(sf * u[basis.nodal_dofs])

if name == "main":
from os.path import splitext
from sys import argv

# save to VTK for visualization in Paraview
m.save(splitext(argv[0])[0] + '.vtk')

[kinnala]

Which solver should I use in scikit-fem for this type of problem?

There is no Newton solver built in scikit-fem, so if you need to perform nonlinear iterations, such solvers must be implemented by the user, and the same applies to time discretization. In fact, there are no linear solvers either and solve (linear solver) is simply a lightweight wrapper to the linear solvers in scipy.sparse package for convenience. I would characterize scikit-fem as a FEM matrix assembly package.

However, Newton-Raphson is probably the most common solver for the purpose you are describing. You just need to create a for loop and assemble the tangent linear system and solve it during each iteration.

[kinnala]

I tried checking the strain tensor by including "print(sym_grad(u))" in the Bilinear form. However, I noticed that the tensor didn't change when I altered the boundary condition from u^1=-0.5 to u^2=-0.5.

Boundary conditions do not affect the matrix unless they are special ones such as Robin-type boundary conditions. The types of boundary conditions you are describing are Dirichlet-type conditions and they are implemented at the linear algebra level, after the matrices have been created, by modifying the contents of the matrix.

Edit: In practice, this modification of the matrix is done by one of the following functions: condense, enforce, penalize, or you can also implement your own quite easily.

[kinnala]

Could anyone explain how to correctly derive the strain tensor?

sym_grad(u) is the infinitesimal strain tensor.

[HoseongJeong]

Thank you so much for your kind response!

I am using a nonlinear plastic material model (class) to reflect the characteristics of steel, where stiffness changes according to the strain history of the material.
This class has internal variables that record the experienced strain.
Currently, I am trying to apply "sym_grad(u)" as an argument to this class (Material_model(sym_grad(u))). However, "sym_grad(u)" provides values related to the derivatives of shape functions at the quadrature points, and even without applying displacement boundary conditions, values greater than 0.002 (the strain at which steel begins plastic behavior) are input into the material model during the assembly process. As a result, the material model records it as if the material has undergone deformation.
In such cases, what would be the appropriate way to input the strain into the bilinear form?

In addition, thank you so much for kindly answering even my trivial and insignificant questions!! Compared to other FEM tools that are difficult to install, scikit-fem is easy to install and very easy to use. I will make sure to widely promote it to my fellow researchers, so that citations will be made actively.

[kinnala]

I believe you are doing an iteration and using the solution from the previous iteration to calculate strain which is used in the bilinear form.

You can do it like this:

x = solve(…) # previous solution
xip = basis.interpolate(x) # interpolate at quadrature points

@BilinearFoem
def bilinf(u, v, w):
    strain = sym_grad(w['xip'])
    return … # write your form here

K = bilinf.assemble(basis, xip=xip) # provide previous solution to the form

[HoseongJeong]

Thanks!!! It is a very clear answer!!!