Nonsmooth problems formulation and solve

When dynamical systems and their interactions have been properly defined inside a model and its nonsmooth dynamical system, a proper formulation for the latter must be chosen, associated to a nonsmooth solver.

Linear nonsmooth problems

\[w = q + M z, M \in R^{n \times n }, q \in R^{n}\]

where \(w \in R^{n}, z \in R^{n}\) are the unknowns.

  • Linear Complementarity Problems (LCP)
\[\begin{split}w = q + M z, M \in R^{n \times n }, q \in R^{n} \\ w \geq 0, z \geq 0, z^{T} w =0\end{split}\]
  • Mixed Linear Complementarity Problems (MLCP)
\[\begin{split}0 = Au + Cv + a\\ z = Du + Bv + b\\ v \geq 0, z \geq 0, z^{T} v =0\end{split}\]

where

\[\begin{split}u \in R^{n}, v \in R^{m}, z \in R^{m} \ the \ unknowns, \\ a \in R^{n}, b \in R^{m} \\ A \in R^{n \times n }, B \in R^{m \times m }\\ C \in R^{n \times m }, D \in R^{m \times n }\end{split}\]
  • Relay
  • Equality
  • AVI
  • 2D or 3D friction contact problem FrictionContact
\[\begin{split}velocity = q + M reaction \\ velocity \geq 0, reaction \geq 0, reaction^{T} velocity =0\end{split}\]

and a Coulomb friction law. With \(velocity \in R^{n}, reaction \in R^{n}\) the unknowns, and \(M \in R^{n \times n }, q \in R^{n}\)

  • multiple-impact problem (OSNSMultipleImpact)
  • primal friction contact problems (GlobalFrictionContact)
\[\begin{split}M velocity = q + reaction \\ localVelocity = H^T velocity + tildeLocalVelocity\\ reaction = H localReaction \\\end{split}\]

and \(localVelocity,localReaction\) belongs to the Coulomb friction law with unilateral contact.

With \(velocity \in R^{n}, reaction \in R^{n}, localVelocity \in R^{m}, localReaction \in R^{m}\) the unknowns, \(M \in R^{n \times n }, q \in R^{n}\). \(tildeLocalVelocity \in R^{m}\) is the modified local velocity (\(e U_{N,k}\)), \(M \in R^{n \times n }, q \in R^{n}, H \in R^{n \times m }\).

  • Generic mechanical problem (GenericMechanical)

Complete problem with bilateral equality, complementarity, impact and friction.

The Simulation process

As for Event-Driven, we introduce level index sets, with level = 0 for first order systems and level=1 for second order systems (this is related to the relative degrees but we won’t get into details about that here).

\(I_0\) is the set of all the potential UnitaryRelations (UR). For second order systems: \(I_1 = \{ ur_\alpha\in I_{0} , y^p_{\alpha} = 0 \}\). Thus, the LCP is built only for unitary relations that belongs to \(I_level\), level=0 for first order and level=1 for second order systems.

Then, the steps of a Moreau’s Time-Stepping simulation will be:

Knowing all values at the beginning of the time step \([t_i,t_{i+1}]\),

-# compute the free solutions -# for \(ur \in I_level\) formalize and solve a LCP -# update the state (according to the possibly LCP results) -# go to next time step

SP::TimeStepping s(new TimeStepping(myModel));
SP::TimeDiscretisation t(new TimeDiscretisation(timeStep,s));

s->initialize();

int N = t->getNSteps(); // Number of time steps

// --- Time loop ---
while(k < N)// for each time step ...
{
// compute xFree, or qFree,vFree
s->computeFreeStep();
// Formalize and solve a LCP
computeOneStepNSProblem("timeStepping");
// Update state, using last computed values
s->update(level); //
// transfer of state i+1 into state i and time incrementation
s->nextStep();
}

Note that all time-independent operators are computed during simulation initialisation.

Customize simulation behavior

Each time ComputeOneStepNS() function, i.e. the numerics solver, is called, it returns an int, giving some information about the convergence of the solver:

By default, when the convergence is not achieved, an exception is throwed and the process stops. Change this behavior is possible by defining a specific function of the form:

//
// your inputFile.cpp
//
void myF(int info, SP::Simulation s)
{
// do what you need ...
}

int main(int argc, char* argv[])
{
//
// ...
SP::TimeStepping your_simulation = ...
your_simulation->setCheckSolverFunction(&myF);

Then after each call to your_simulation->computeOneStepNS(…), the function myF will be called. That may be usefull to change the solver type, the tolerance or whatever is needed.