Homec4science
Diffusion tamaas (ci-fix-static-checks) Browse

tamaas/doc/sphinx/source/contact.rst
bd65ceb3903fci-fix-static-checks

Last Change
Raw File

contact.rst
View Options

Solving contact
===============
The resolution of contact problems is done with classes that inherit from :cpp:class:`tamaas::ContactSolver`. These usually take as argument a :cpp:class:`tamaas::Model` object, a surface described by a :cpp:class:`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 :cpp:class:`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.
Contact with adhesion
---------------------
The second algorithm hidden in :cpp:class:`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. <https://doi.org/10.1007/s00466-017-1392-5>`_). 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 :cpp:class:`tamaas::Functional`::
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 * np.sum(np.exp(-gap / self.rho))
def computeGradF(self, gap, gradient):
gradient += self.gamma * np.exp(-gap / self.rho) / self.rho
This example is actually equivalent to :cpp:class:`tamaas::functional::ExponentialAdhesionFunctional`.
Tangential contact
------------------
For frictional contact, several solver classes are available. Among them, :cpp:class:`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 :cpp:class:`tamaas::EPSolver`. They make use of the non-linear solvers available in scipy:
:py:class:`DFSANESolver <tamaas.nonlinear_solvers.DFSANESolver>`
Implements an interface to :py:func:`scipy.optimize.root` wiht the DFSANE method.
:py:class:`NewtonKrylovSolver <tamaas.nonlinear_solvers.NewtonKrylovSolver>`
Implements an interface to :py:func:`scipy.optimize.newton_krylov`.
These solvers require a residual vector to cancel, the creation of which is handled with :cpp:class:`tamaas::ModelFactory`. Then, an instance of :cpp:class:`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
residual = tm.ModelFactory.createResidual(model,
sigma_y=0.1 * model.E,
hardening=0.01 * model.E)
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, :cpp:func:`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 :cpp:func:`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 :py:class:`ToleranceManager <tamaas.nonlinear_solvers.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)
c4science · Help