Solving contact
The resolution of contact problems is done with classes that inherit from tamaas::ContactSolver
. These usually take as argument a tamaas::Model
object, a surface described by a tamaas::Grid
or a 2D numpy array, and a tolerance. We will see the specificities of the different solver objects below.
Elastic contact
The most common case is normal elastic contact, and is most efficiently solved with tamaas::PolonskyKeerRey
. The advantage of this class is that it combines two algorithms into one. By default, it considers that the contact pressure field is the unknown, and tries to minimize the complementary energy of the system under the constraint that the mean pressure should be equal to the value supplied by the user, for example:
# ...
solver = tm.PolonskyKeerRey(model, surface, 1e-12)
solver.solve(1e-2)
Here the average pressure is 1e-2
. The solver can also be instanciated by specifying the the constraint should be on the mean gap instead of mean pressure:
solver = tm.PolonskyKeerRey(model, surface, 1e-12, constraint_type=tm.PolonskyKeerRey.gap)
solver.solve(1e-2)
The contact problem is now solved for a mean gap of 1e-2
. Note that the choice of constraint affects the performance of the algorithm.
Non-periodic contact
When the primal type is pressure, it is possible to solve a non-periodic normal contact problem by first registering a non-periodic integral operator (based on the DC-FFT method) then telling the contact solver to use it in its energy functional:
tm.ModelFactory.registerNonPeriodic(model, 'dcfft')
solver.setIntegralOperator('dcfft')
solver.solve(1e-2)
Note however that the solver still shifts the displacement field so that the gap is everywhere positive and zero in the contact zones. This loses the information of the absolute displacement relative to the initial position of the surface. This can be recovered post-solve with:
model.operators['dcfft'](model.traction, model.displacement)
The example nonperiodic.py
shows an Hertzian contact example.
Computing solutions for loading sequence
The module tamaas.utils
defines a convenience function
load_path
, which generates solution models
for a sequence of loads. This allows lazy evalution and reduces boiler-plate:
from tamaas.utils import load_path
loads = np.linspace(0.01, 0.1, 10)
for model in load_path(solver, loads):
... # do some computation on model, e.g. compute contact clusters
The function load_path
accepts the following
optional arguments:
verbose
Prints solver output (i.e. iteration, cost function and error)
callback
A function to call after each solve, before the next load step. Useful for dumping model in generator expressions.
Contact with adhesion
The second algorithm hidden in tamaas::PolonskyKeerRey
considers the gap to be the unknown. The constaint on the mean value can be chosen as above:
solver = tm.PolonskyKeerRey(model, surface, 1e-12,
primal_type=tm.PolonskyKeerRey.gap,
constraint_type=tm.PolonskyKeerRey.gap)
solver.solve(1e-2)
The advantage of this formulation is to be able to solve adhesion problems (Rey et al.). This is done by adding a term to the potential energy functional that the solver tries to minimize:
adhesion_params = {
"rho": 2e-3,
"surface_energy": 2e-5
}
adhesion = tm.ExponentialAdhesionFunctional(surface)
adhesion.setParameters(adhesion)
solver.addFunctionalterm(adhesion)
solver.solve(1e-2)
Custom classes can be used in place of the example term here. One has to inherit from tamaas::Functional
:
import numpy
class AdhesionPython(tm.Functional):
"""
Functional class that extends a C++ class and implements the virtual
methods
"""
def __init__(self, rho, gamma):
super().__init__(self)
self.rho = rho
self.gamma = gamma
def computeF(self, gap, pressure):
return -self.gamma * numpy.sum(np.exp(-gap / self.rho))
def computeGradF(self, gap, gradient):
gradient += self.gamma * numpy.exp(-gap / self.rho) / self.rho
This example is actually equivalent to tamaas::functional::ExponentialAdhesionFunctional
.
Viscoelastic contact
The solver class tamaas::MaxwellViscoelastic
solves the
quasi-static problem of a viscoelastic material in contact with a rigid rough
surface. The viscoelastic behaviour is described with a generalized Maxwell
model, i.e. a
parallel assembly of Maxwell elements, with one purely elastic chain, defined by
the elastic properties of the tamaas::Model
object:
# Defining the elastic branch (i.e. the behavior at t = ∞)
model.E = 3
model.nu = 0.5
# Characteristic times of the relaxation function
times = [0.1, 1, 10]
# Shear moduli for each branch of the model
shear_moduli = [3, 3/2, 1/3]
# Time step
Δt = 1e-3
# Solver instanciation
solver = tm.MaxwellViscoelastic(model, surface, 1e-10,
time_step=Δt,
shear_moduli=shear_moduli,
characteristic_times=times)
# Solve one timestep with given load
solver.solve(load)
Warning
The current viscoelastic formulation is limited to incompressible
(model.nu = 0.5
) materials.
Tip
The viscoelastic solver keeps track of loading and deformation history.
Use solver.reset()
to erase this data and start loading from a rest
configuration.
Tangential contact
For frictional contact, several solver classes are available. Among them, tamaas::Condat
is able to solve a coupled normal/tangential contact problem regardless of the material properties. It however solves an associated version of the Coulomb friction law. In general, the Coulomb friction used in contact makes the problem ill-posed.
Solving a tangential contact problem is not much different from normal contact:
mu = 0.3 # Friction coefficient
solver = tm.Condat(model, surface, 1e-12, mu)
solver.max_iter = 5000 # The default of 1000 may be too little
solver.solve([1e-2, 0, 1e-2]) # 3D components of applied mean pressure
Elasto-plastic contact
For elastic-plastic contact, one needs three different solvers: an elastic contact solver like the ones described above, a non-linear solver and a coupling solver. The non-linear solvers available in Tamaas are implemented in python and inherit from the C++ class tamaas::EPSolver
. They make use of the non-linear solvers available in scipy:
DFSANESolver
Implements an interface to
scipy.optimize.root()
wiht the DFSANE method.NewtonKrylovSolver
Implements an interface to
scipy.optimize.newton_krylov()
.
These solvers require a residual vector to cancel, the creation of which is handled with tamaas::ModelFactory
. Then, an instance of tamaas::EPICSolver
is needed to couple the contact and non-linear solvers for the elastic-plastic contact problem:
import tamaas as tm
from tamaas.nonlinear_solvers import DFSANESolver
# Definition of modeled domain
model_type = tm.model_type.volume_2d
discretization = [32, 51, 51] # Order: [z, x, y]
flat_domain = [1, 1]
system_size = [0.5] + flat_domain
# Setup for plasticity
material = tm.materials.IsotropicHardening(model,
sigma_y=0.1 * model.E,
hardening=0.01 * model.E)
residual = tm.Residual(model, material)
epsolver = DFSANESolver(residual)
# ...
csolver = tm.PolonskyKeerRey(model, surface, 1e-12)
epic = tm.EPICSolver(csolver, epsolver, 1e-7, relaxation=0.3)
epic.solve(1e-3)
# or use an accelerated scheme (which can sometimes not converge)
epic.acceleratedSolve(1e-3)
By default, tamaas::EPICSolver::solve()
uses a relaxed fixed point. It can be tricky to make it converge: you need to decrease the relaxation parameter passed as argument of the constructor, but this also hinders the convergence rate. The function tamaas::EPICSolver::acceleratedSolve()
does not require the tweaking of a relaxation parameter, so it can be faster if the latter does not have an optimal value. However, it is not guaranteed to converge.
Finally, during the first iterations, the fixed point error will be large compared to the error of the non-linear solver. It can improve performance to have the tolerance of the non-linear solver (which is the most expensive step of the fixed point solve) decrease over the iterations. This can be achieved with the decorator class ToleranceManager
:
from tamaas.nonlinear_solvers import ToleranceManager, DFSANESolver
@ToleranceManager(start=1e-2,
end=1e-9,
rate=1 / 3)
class Solver(DFSANESolver):
pass
# ...
epsolver = Solver(residual)
# or
epsolver = ToleranceManager(1e-2, 1e-9, 1/3)(DFSANESolver)(residual)
Using PETSc Solvers
Two C++ classes in Tamaas interface with PETSc solvers, provided
use_petsc=True
was used at compilation:
tamaas::petsc::OptimizationSolver
uses PETSc’s Toolkit for Advanced Optimization (TAO) to solve normal contact problems.tamaas::petsc::NonlinearSolver
uses PETSc’s SNES interface to solve non-linear problems (e.g. plasticity).
Both classes accept PETSc’s command-line parameter lists, so that everything about these solvers can be adjusted at run-time by the user, e.g. the linear solver, the tolerance, the line search strategy, etc.
Solving contact
Due to the definition of a standard quadratic constrained program in TAO, the contact
problem solved by tamaas::petsc::OptimizationSolver
is constrained
only by the condition that the gap be non-negative. This means that the average
pressure cannot be prescribed. Instead, the contact surface and be uniformly
shifted by \(\delta\) to control the contact (positive values of the mean
height generally induce larger contact), with contact occuring when:
Custom Contact Solvers
The tamaas::ContactSolver
class can be derived in Python so that
users can interface with Scipy’s scipy.optimize.minimize()
function or
PETSc’s solvers accessible in the Python interface. See
examples/scipy_penalty.py
for an example on how to proceed.