Friction is a complex physical phenomenon that involves the characteristics of the surface such as surface roughness, temperature, normal stress, and relative velocity. An example of this complexity is that quite often in contact problems neutral lines develop. This means that along a contact surface, the material flows in one direction in part of the surface and in the opposite direction in another part of the surface. Such neutral lines are, in general, not known a priori.
The actual physics of friction and its numerical representation continue to be topics of research. Currently, in CivilFEM the modeling of friction has basically been simplified to Coulomb friction model.
Coulomb Friction
The Coulomb model can be characterized by:
and
where
σt is the tangential (friction) stress
σn is the normal stress
μ is the friction coefficient
t is the tangential vector in the direction of the relative velocity:
Similarly, the Coulomb model can also be written in terms of nodal forces instead of stresses:
and
where
is the tangential (friction) force
is the normal force
In CivilFEM the friction models is based on forces.
With the Coulomb friction model the integration point stresses are first extrapolated to the nodal points and then transformed, so that a direct component is normal to the contacted surface. Given this normal stress and the relative sliding velocity, the tangential stress is then evaluated and a consistent nodal force is calculated.
For a given normal stress or normal force, the friction stress or force has a step function behavior based upon the value of the relative sliding velocity 
or the tangential relative incremental displacement 
as outlined in figure below for a 2-D case, where the relative velocity and incremental displacement are scalar values.
Since this discontinuity in the friction value may easily cause numerical difficulties, different approximations of the step function have been implemented. They are graphically represented in figure below and bilinear model will be discussed as is the one available in CivilFEM.
The the bilinear model is based on relative tangential displacements. Instead of defining special constraints to enforce sticking conditions, the bilinear model assumes that the stick and slip conditions correspond to reversible (elastic) and permanent (plastic) relative displacements, respectively. The clear resemblance with the theory of elasto-plasticity will be used to derive the governing equations.
First, Coulomb’s law for friction is expressed by a slip surface φ:
The stick or elastic domain is given by φ < 0 , while φ > 0 is physically impossible.
Next, the rate of the relative tangential displacement vector is split into an elastic (stick) and a plastic (slip) contribution according to:
and the rate of change of friction force vector is related to the elastic tangential displacement by:
in which matrix D is given by:
With 
the slip threshold or relative sliding displacement below which sticking is simulated.
The value of is by default determined by the program as 0.0025 times the average edge length of the finite elements defining the deformable contact bodies.
Now attention is paid to the case that, given a tangential displacement vector, the evolution of would result in aphysically impossible situation, so φ > 0. This implies that the plastic or slip contribution must be determined.
It is assumed that the direction of the slip displacement rate is given by the normal to the slip flow potential , given by:
So that, by indicating the slip displacement rate magnitude as 
:
with the slip surface, φ , different from the slip flow potential, 
an analogy to nonassociative plasticity can be observed.
Since a ‘force point’ can never be outside the slip surface, it is required that:
In this way, the magnitude of the slip rate can be determined. To this end, the equations above can be combined to:
or:
Utilizing this result, the final set of rate equations reads:
Similar to non-associative plasticity, matrix 
will generally be non-symmetric. In CivilFEM, a special procedure has been used which results in a symmetric matrix, while maintaining sufficient numerical stability and rate of convergence.
The bilinear model also uses an additional check on the convergence of the friction forces, which has been achieved if the following equation is fulfilled:
Where 
is the current total friction force vector (the collection of all nodal contributions), 
is the total friction force vector of the previous iteration and is the friction force tolerance, which has a default value of 0.05.
If a node comes into contact and the structure is still stress-free, then the friction stiffness matrix according to the derivation above will still be zero. This could result in an ill-conditioned system during the next solution of the global set of equations. To avoid this problem, the initial friction stiffness will be based on the average contact body stiffness (following from the trace of the matrix defining the material behavior), which is determined during increment 0 of the analysis.
