# Solving hertzian contact simple equation

Hello everyone, i am new using fenics, i have development a short script to model a hertzian contact and stresses in a cad geometry converted to xml file using dolfin, my question is what am i doing wrong in my code because Fenics says the next error: “o Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F’.
Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?”

i provide my code here, i appreciate your helping:

``````from fenics import *

# Crear malla a partir de un archivo XML
# Definir la malla y el espacio de funciones
filename = '/home/jalastrag/fenics_proj/mesh.xml'
mesh = Mesh(filename)

# Definir los espacios de funciones para la deformación (u) y la presión (p)
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)

# Definir las condiciones de contorno (aquí se asumen condiciones de Dirichlet homogéneas)
bc_u = DirichletBC(V, Constant((0, 0, 0)), "on_boundary")
bc_p = DirichletBC(Q, Constant(0), "on_boundary")

# Definir las funciones de prueba y de prueba de incremento
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Definir los parámetros elásticos y de Poisson
E = Constant(1e6)  # Módulo de Young
nu = Constant(0.3)  # Coeficiente de Poisson

# Definir la ecuación elástica lineal

# Definir la ecuación de la presión (puede variar dependiendo de la formulación)
F_pressure = -p*q*dx

# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure

# Resolver el problema variacional
u_solution = Function(V)
p_solution = Function(Q)
solve(F == 0, u_solution, bc_u)
solve(F_pressure == 0, p_solution, bc_p)

# Guardar resultados en archivos VTK
File("deformacion.pvd") << u_solution
File("presion.pvd") << p_solution
``````

To reproduce your error, I first replace the above quoted part by

``````mesh = UnitCubeMesh(16, 16, 16)
``````

as you did not provide your `mesh.xml` file.

Your problem comes from the following one:

and

as you are solving a coupled problem.

If you replace

by

``````solve(F_elastic == 0, u_solution, bc_u)
solve(F_pressure == 0, p_solution, bc_p)
``````

you will find it reports the nonlinear variational problem which should come with a `residual F` or `Jacobian J`.
Please consult to the FEniCS demo:
https://fenicsproject.org/olddocs/dolfin/2019.1.0/python/demos/hyperelasticity/demo_hyperelasticity.py.html

2 Likes

Hi.

First, since your are solving for `u` and `p`, is better to use a `MixedElement` approach to avoid the error

``````"Found different Arguments with same number and part.
Did you combine test or trial functions from different spaces?"
``````

Second, the error `Differentiating residual from F to obtain Jacobian J = F’.` comes from the fact that your are trying to use the automatic non-linear solver but using `TrialFunctions`.

The modified code below works and converges in 1 iteration. Note that I have used `sol` as a function of the mixed space `V`, followed by `u, p = split(sol)`, which is typical when solving Navier-Stokes equations in Legacy Fenics.

``````from dolfin import *

mesh = UnitCubeMesh(5, 5, 5)

# Definir los espacios de funciones para la deformación (u) y la presión (p)
P2 = VectorElement('CG', mesh.ufl_cell(), 2)
P1 = FiniteElement('CG', mesh.ufl_cell(), 1)

element = MixedElement([P2, P1])  # Con multiplicador
V = FunctionSpace(mesh, element)

# Definir las condiciones de contorno (aquí se asumen condiciones de Dirichlet homogéneas)
bc_u = DirichletBC(V.sub(0), Constant((0, 0, 0)), "on_boundary")
bc_p = DirichletBC(V.sub(1), Constant(0), "on_boundary")

# Definir las funciones de prueba y de prueba de incremento
sol = Function(V)
u, p = split(sol)
v, q = TestFunctions(V)

# Definir los parámetros elásticos y de Poisson
E = Constant(1e6)  # Módulo de Young
nu = Constant(0.3)  # Coeficiente de Poisson

# Definir la ecuación elástica lineal

# Definir la ecuación de la presión (puede variar dependiendo de la formulación)
F_pressure = -p*q*dx

# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure

# Resolver el problema variacional
# sol = Function(V)

solve(F == 0, sol, bc_u)

uh, ph = sol.split(True)

# Guardar resultados en archivos VTK
File("deformacion.pvd") << uh
File("presion.pvd") << uh

``````

Cheers.

2 Likes

Hello thank you for your answer sir, i apologize for don’t share my mesh because i was not be able to upload, soon i hope do it, how ever i appreciate your answer, this routine that i showed my intention is to make a nonlinear problem about wire helix spring and study behavior of them in hertzian contacts, do you have any docummentation about it where do i can find simulation using fenics of that kind of problem?.. again thank you very much for your answer

Hello Jesus thank oyu very much for your answer i used your script how ever i can’t see in paraview any answer (cube response in deformation or stresses) could you explain me please where can i change number or boundary condition to see any result? how i wrote before my intettion (i am a new user fenics) is to make a simulatin of wire helix spring studying their hertzian contact, if you have some documentation i will appreciate it, thank you for share this code, and a nother question this code repsonse means that my results are zero? : olving nonlinear variational problem.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.

Your system is not doing any work. You can add a body force:

``````fbody = Expression(('sin(x[1])', 'cos(x[0])', 'x[0]*x[1]'), degree=2)
F_body = inner(fbody, v)*dx
# Definir la ecuación de equilibrio variacional
F = F_elastic + F_pressure - F_body
``````

On the other hand, since you are new to FEniCS, I recommend you to use the new DOLFINx, which has a really nice tutorial (or this) and a very rich documentation. There are some results here about hertzian contact in DOLFINx, and here for the Legacy FEniCS.

Hope this helps.

2 Likes