Including Surface Roughness in Contact via a Modified Gap Function in GetFEM

[Djoumessi]

Hello everyone,

I am currently exploring how to model surface roughness in contact mechanics using GetFEM, and I would like to know if anyone has already attempted something similar.

My idea is to incorporate the roughness directly into the definition of the normal gap. Instead of considering the obstacle surface as perfectly smooth, I would compute an “effective gap” of the form:

gn​(x)=(u_2​−u_1^*​)⋅n−h(x)

where h(x) represents the roughness profile mapped onto the contact surface.
The function h(x) could come from:

• a point cloud obtained from surface measurements (e.g., profilometry), or
• a synthetic/fractal roughness function generated numerically.

Once the effective gap is defined, the idea is to apply any of the standard contact formulations (penalty, augmented Lagrangian, Nitsche, etc.) without modifying the core solver—only the gap computation changes.

Before going further, I would like to know:

• Has anyone already implemented contact with rough surfaces in GetFEM?
• Is there an existing example or recommended workflow for injecting a custom obstacle function h(x) into the contact gap computation?
• Any advice on potential pitfalls or good practices when working with roughness fields on the contact boundary?

Any insight or references would be greatly appreciated.
Thank you!

[Konstantinos.Poulios]

Augmented Lagrangian and other exact contact methods are not relevant for modeling of rough contact. Rough contact modelling is basically equivalent to penalty contact or barrier function contact with some tailored law for the pressure vs gap function.

Rough contact is therefore easier to implement than non-rough contact. You do not need to redefine the gap function g, you just need to add a contact potential W(g) to your model, with units of energy per area. All the rough contact modelling goes into defining the W(g) function. One characteristic of this function is that it should go to infinity for g going towards 0. So it really needs to be a barrier function in order to represent the real system.

If you look at the old papers by Archard you can probably find good inspiration for defining this function.

[Djoumessi]

Hello,

Thank you very much for your insightful explanation and for pointing me toward the contact potential W(g) approach.
This clarifies a lot and gives me a clear direction to model rough contact more effectively without redefining the gap function itself.

I will definitely look into Archard’s papers for inspiration on suitable barrier functions and start implementing a penalty/barrier function based contact energy in GetFEM.

Thanks again for sharing your expertise—it’s really helpful!

Best regards,

[Djoumessi]

Hi Konstantinos,

I am coming back to you with a few additional questions in order to better understand how to integrate a rough surface into a contact model.

From our previous discussion, I understood that the simplest approach is to implement a statistical or fractal-type roughness function, such as a Greenwood–Williamson formulation, directly within the model. While searching for existing implementations in GetFEM, I came across your article published in Wear in 2013:
Implementation and applications of a finite-element model for the contact between rough surfaces - ScienceDirect

I had a few questions regarding its implementation. Based on our earlier exchange, I initially thought that a penalty-type approach would be the easiest to implement. However, in your article, you adopt a Lagrange multiplier formulation. I would like to know whether this approach is actually simpler to implement within GetFEM, or whether it mainly depends on the type of contact problem considered.

In addition, I did not fully understand how the surface topography used for the roughness was incorporated into the finite element model. Once the rough surface was generated, did you introduce it into the model as an external field, or in some other way? A clarification on this aspect would be extremely helpful.

Finally, if the implementation were to be done in a more complex setting for example, contact between two cylinders would the same approach remain valid?

Any additional guidance you could provide would be greatly appreciated.
Thank you in advance for your help.

Best regards,

[Konstantinos.Poulios]

well, in the case of the paper it was possible to model the roughness directly because we were only interested in a small representative sample, something like 0.5mm x 0.5mm. If you have some larger surfaces in contact, you cannot afford discretizing the roughness topography with your finite element mesh.

For larger surfaces in contact you should instead model the stiffness of the roughness as part of the contact formulation as I described before.

What is the actual application you are dealing with?

[Djoumessi]

In my case, I am working on a larger-scale problem involving two cylinders in contact along their lateral surfaces, as well as situations of self-contact. Given the size of the surfaces, it is not feasible to discretize the roughness topography directly with the finite element mesh.

For this reason, I am focusing on a Greenwood–Williamson–type approach, based on a statistical representation of surface roughness. My main difficulty is understanding how to incorporate this roughness law into GetFEM, within the contact formulation.

I also plan to evaluate this approach in a large deformation setting later on, but for now I would like to clarify how to integrate the roughness model.

[Konstantinos.Poulios]

you need to know some special syntax in order to define the gap function. I will provide here directly the version for large deformations:

md.add_macro("W(s)","......s") # some function of s that is derived from a rough contact model
md.add_raytracing_transformation("raytrace", release_dist) # release_dist is how far from the surface to search for contact
md.define_variable_group("u", "u1", "u2") # if you have different displ. fields for each body
md.add_master_contact_boundary_to_raytracing_transformation\
  ("raytrace", m1, "u", MASTER_RG)
md.add_slave_contact_boundary_to_raytracing_transformation\
  ("raytrace", m2, "u", SLAVE_RG)
md.add_macro("normal", "Normalized(Inv(F(u2)').Normal)") # normal on the slave surface
md.add_macro("gap", "(Interpolate(X,raytrace)+Interpolate(u,raytrace)-X-u2).normal")
md.add_nonlinear_term(mim_slave, "Interpolate_filter(raytrace, W(gap), 1)", SLAVE_RG)

this is not tested code, just so that you get the idea conceptually

[Djoumessi]

Hi Konstantinos,

Thank you very much for this code snippet. I’ll test it and get back to you.
This is exactly the line that was missing in my approach:

md.add_nonlinear_term(mim_slave, "Interpolate_filter(raytrace, W(gap), 1)", SLAVE_RG);

GetFEM is really powerful.
I’ll come back to you very soon.

Best regards,
René

[Konstantinos.Poulios]

@Djoumessi

the line that was missing in my approach:

this is the line which is not tested, it will probably give you an error, but we can fix it by deriving the corresponding 1st order expressions from this 0-order expression by hand. It is not very difficult.

[Konstantinos.Poulios]

for example, this should theoretically work

md.add_macro("W(gap)", "1/2*c*sqr(neg_part(gap))")
md.add_nonlinear_term(mim_slave, "Interpolate_filter(raytrace, W(gap), 1)",  SLAVE_RG) 

but it does not because neg_part is not twice differentiable, but the W function is actually twice differentiable.

For this reason we need to provide the 1st order expression instead:

md.add_nonlinear_term(mim_slave,
  "Interpolate_filter(raytrace, -c*neg_part(gap)*(Interpolate(Test_u,raytrace)-Test_u2).normal2, 1)",
  SLAVE_RG) 

you can see how it works in practice in this example

#!/usr/bin/env python
# -*- coding: UTF8 -*-
#
############################################################################
import sys
import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

π = np.pi

# Input data
R = 50.     # ring radius
r = 5.      # cross section radius

NC = 26     # Number of elements along the ring circumference
Nr1 = 2     # Number of elements in the cross section core
Nr2 = 2     # Number of elements in the cross section peel

E = 1e5      # Young's modulus
nu = 0.45    # Poisson's ratio

mult_fem_order = 1  # contact multipliers finite element order
integration_degree = 5         # 27 gauss points per hexa
integration_contact_degree = 7 # 16 gauss points per quad

c_penal = 1e-2*E # penalty parameter
release_dist = r

# auxiliary constants
EXTERIOR_RG = 1
HANDLE1_RG = 2
HANDLE2_RG = 3

# Mesh
m_ = gf.Mesh("import", "structured_ball",
              "GT='GT_QK(2,2)';ORG=[0,0];SIZES=[%g];NSUBDIV=[%i,%i]"
               % (r, Nr1, Nr2))
m0 = gf.Mesh("prismatic", m_, NC, 2) # second order

pts = m0.pts()
for i in range(pts.shape[1]):
  x = pts[0,i]
  z = pts[2,i]
  pts[0,i] = (R+x)*np.cos(2*π*z)
  pts[2,i] = (R+x)*np.sin(2*π*z)
m0.set_pts(pts)

m1 = gf.Mesh("empty", 3) 
m1.merge(m0, tol=1e-3) # this trick is used just for gluing the open ends of the ring in m0

m1.set_region(EXTERIOR_RG, m1.outer_faces())
m1.set_region(HANDLE1_RG, m1.outer_faces_in_box([-R-1.1*r,-1.1*r,-0.3*R],
                                                [-R+1.1*r, 1.1*r, 0.3*R]))
m1.set_region(HANDLE2_RG, m1.outer_faces_in_box([R-1.1*r,-1.1*r,-0.3*R],
                                                [R+1.1*r, 1.1*r, 0.3*R]))
m1.optimize_structure()
m1.export_to_vtu("mesh1.vtu")

m2 = gf.Mesh("clone", m1)
m2.translate([-R,2.5*r,0])
m2.export_to_vtu("mesh2.vtu")

# FEM
mf1 = gf.MeshFem(m1, 3)
mf1.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))

mf2 = gf.MeshFem(m2, 3)
mf2.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))

mfmult1 = gf.MeshFem(m1, 3)
mfmult1.set_classical_fem(mult_fem_order)

mfmult2 = gf.MeshFem(m2, 3)
mfmult2.set_classical_fem(mult_fem_order)

mim1 = gf.MeshIm(m1, integration_degree)
mim2 = gf.MeshIm(m2, integration_degree)
mim_contact1 = gf.MeshIm(m1, integration_contact_degree)
mim_contact2 = gf.MeshIm(m2, integration_contact_degree)

# Model
md = gf.Model("real")

md.add_fem_variable("u1", mf1)
md.add_fem_variable("u2", mf2)

# Elasticity
md.add_initialized_data("K", E/(3*(1-2*nu)))
md.add_initialized_data("G", E/(2*(1+nu)))
md.add_macro("F(u)", "Id(3)+Grad_u")
md.add_macro("tau(u)", "K*log(Det(F(u)))*Id(3)+G*pow(Det(F(u)),-2/3)*Deviator(F(u)*F(u)')")
md.add_nonlinear_term(mim1, "(tau(u1)*Inv(F(u1)')):Grad_Test_u1")
md.add_nonlinear_term(mim2, "(tau(u2)*Inv(F(u2)')):Grad_Test_u2")

# Lagrange multipliers for handles
md.add_variable("q11",3) # ring 1 handle 1
md.add_variable("q12",3) # ring 1 handle 2
md.add_variable("q21",3) # ring 2 handle 1
md.add_variable("q22",3) # ring 2 handle 2
md.add_initialized_data("p11", [0,0,0]) # we will keep this fixed
md.add_initialized_data("p12", [0,0,0]) # we will keep this fixed only at the beginning
md.add_initialized_data("p21", [0,0,0]) # we will drive this only at the beginning
md.add_initialized_data("p22", [0,0,0]) # we will drive this till the end

md.add_linear_term(mim1, "(u1-p11).q11", HANDLE1_RG)
md.add_linear_term(mim1, "(u1-p12).q12", HANDLE2_RG)
md.add_linear_term(mim2, "(u2-p21).q21", HANDLE1_RG)
md.add_linear_term(mim2, "(u2-p22).q22", HANDLE2_RG)

md.add_initialized_data("penal", 0.)
md.add_linear_term(mim1, "1e4*(u1-p11).Test_u1", HANDLE1_RG) # add some elastic support to handle 1 of ring 1 to avoid twisting
md.add_linear_term(mim1, "penal*(u1-p12).Test_u1", HANDLE2_RG) # add some elastic support to handle 2 of ring 1 to avoid twisting
md.add_linear_term(mim2, "penal*(u2-p21).Test_u2", HANDLE1_RG) # add some elastic support to handle 1 of ring 2 to avoid twisting
md.add_linear_term(mim2, "1e4*(u2-p22).Test_u2", HANDLE2_RG) # add some elastic support to handle 2 of ring 2 to avoid twisting

# Contact condition
md.add_initialized_data("c", c_penal)
md.add_raytracing_transformation("raytrace", release_dist)
md.define_variable_group("u", "u1", "u2")
md.add_master_contact_boundary_to_raytracing_transformation("raytrace", m1, "u", EXTERIOR_RG)
md.add_slave_contact_boundary_to_raytracing_transformation("raytrace", m2, "u", EXTERIOR_RG)

md.add_macro("normal2", "Normalized(Inv(F(u2)').Normal)")
md.add_macro("gap2", "(Interpolate(X,raytrace)+Interpolate(u,raytrace)-X-u2).normal2")
#md.add_macro("W(gap)", "1/2*c*sqr(neg_part(gap))") # dW/dgap = -c*neg_part(gap)
md.add_nonlinear_term(mim_contact2,
                      "Interpolate_filter(raytrace, -c*neg_part(gap2)*(Interpolate(Test_u,raytrace)-Test_u2).normal2, 1)",
                      EXTERIOR_RG) 

print("Model dofs: %i" % md.nbdof())
print("Displacement fem dofs: %i" % (mf1.nbdof()+mf2.nbdof()))

for step in range(90):
  if step <= 30: # initial phase
    t = step/30
    md.set_variable("p21", [1.8*R*(max(t,min(4*t,0.2))),
                            -5*r*min(t/0.4,1), #-5*r*(np.cos(π*t)-1)/2,
                            0])
    md.set_variable("p22", [1.8*R*max(0,(t-0.4)/0.6)**2,
                            -5*r*min(t/0.4,1),
                            0])
  elif step <= 38: # second phase, releasing handles 21 and 22
    if step == 31:
      md.disable_variable("q12")
      md.disable_variable("q21")
      md.set_variable("q12", [0,0,0])
      md.set_variable("q21", [0,0,0])
    md.set_variable("penal", (1e2,3e1,1e1,3e0,1e0,1e-1,1e-3,0.)[step-31])
  else: # final phase
    t = 1+(step-38)/100
    if t > 1.25: # slow down
      t = 1.25 + 0.2*(t-1.25)
    md.set_variable("p22", [1.8*R*max(0,(t-0.4)/0.6)**2,
                            -5*r,
                            0])
  print("Step", step, "t =", t)
  md.solve("noisy", "max_iter", 20, "max_res", 1e-8,
           "lsearch", "simplest", "alpha max ratio", 3., "alpha min", 0.02, "alpha mult", 0.6)
  mf1.export_to_vtu(f"ring1_{step}.vtu", mf1, md.variable("u1"), "Displacements")
  mf2.export_to_vtu(f"ring2_{step}.vtu", mf2, md.variable("u2"), "Displacements")

This example implements the standard penalty approach which is the blue curve in the graph below. Your task is to replace the blue curve with one that looks like the red one:

[Djoumessi]

Hi Konstantinos,

Thank you very much for your explanations. They allowed me to finalize my first example, which for now only implements a standard penalty formulation. I am currently working on introducing surface roughness using the Greenwood–Williamson model of the form

F_{n}(h) \equiv \frac{1}{\sqrt{2\pi}} \int_{h}^{\infty} (s - h)^{n} \, e^{- \tfrac{1}{2} s^{2}} \, ds, \qquad n = 0,\ \tfrac{1}{2},\ 1,\ \tfrac{3}{2}

I will come back to you if I have any additional questions.
Thank you in advance.

Best regards,
Thierry

[Djoumessi]

Hello Konstantinos,

I have made some progress with this model.

Instead of using the classical Greenwood–Williamson (GW) model directly, I implemented the discrete approach proposed by Ciavarella in the paper:

A “re-vitalized” Greenwood and Williamson model of elastic contact between fractal surfaces.

As a first step, I generated a rough surface composed of asperities using a Mandelbrot–Weierstrass (MW) fractal model. I relied on an existing code snippet for this part, which I can share if needed.

Then, each asperity is treated individually, and the contact force is computed using the classical Hertzian contact law. The total contact response is obtained by discrete summation, following Ciavarella’s approach, and is finally injected into a GetFEM macro to be used within the contact formulation.

Below, I include the main code excerpts related to:

• the generation of the fractal roughness,

#----- using Fractal model and discret integral idee from re-vitalized GW -----#

    # Grille
    L = 1e-3    # size (m)
    N = 256     # resolution
    x = np.linspace(0, L, N)
    y = np.linspace(0, L, N)
    X, Y = np.meshgrid(x, y)

    # Parameters fractals (Weierstrass-like)
    A0 = 100e-9     # amplitude  (m)
    b = 2.0         # frequency factor (>1)
    H = 0.8         # dimension of Hausdorff approx.
    M = 8           

    # surface
    heights = np.zeros((N, N))
    for n in range(M):
        kx = (2*np.pi/L) * (b**n)
        ky = (2*np.pi/L) * (b**n)
        An = A0 * (b**(-n*H))
        phase_x = np.random.uniform(0, 2*np.pi)
        phase_y = np.random.uniform(0, 2*np.pi)
        heights += An * (np.cos(kx*X + phase_x) + np.cos(ky*Y + phase_y))

    # Normalization 
    heights -= heights.mean()
    sigma_target = heights.std()  

    # Visualisation 3D 
    fig = plt.figure(figsize=(6,5))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(X*1e3, Y*1e3, heights*1e9, cmap='viridis', linewidth=0, antialiased=False)
    ax.set_xlabel("x (mm)")
    ax.set_ylabel("y (mm)")
    ax.set_zlabel("aperities heigh (nm)")
    #plt.title("Surface rugueuse fractale (Mandelbrot–Weierstrass)")
    plt.tight_layout()
    plt.show()
    # Histogramm to obtain z_i, w_i
    Nbins = 64  # number of bins 
    hist, edges = np.histogram(heights.flatten(), bins=Nbins, density=True)
    z_i = 0.5*(edges[1:] + edges[:-1])        # bins center
    dz = np.diff(edges)                        # 
    w_i = hist * dz
    # Mechanical parameter
    Eeff = c_penal #E/(2*(1+nu))         
    Rasp = 10e-6         # radius of asperities should be modify (not constant)
    eta  = (Nbins / L**2)  # asparaties density
    #  global Constant GW (Hertz)
    Cgw = (4.0/3.0) * Eeff * np.sqrt(Rasp) * eta

• the computation of the contact forces,
• and their integration into the GetFEM macro.

terms = []
    for zi_val, wi_val in zip(z_i, w_i):
        # Each term: w_i * pow(max(neg_part(g) - z_i,0), 1.5)
        #terms.append(f"{wi_val:.8e}*pow(max(neg_part(g)-{zi_val:.8e},0),1.5)") # max is not continuous
        terms.append(f"{wi_val:.8e}*pow(pos_part(neg_part(g)-{zi_val:.8e}),1.5)")

    sum_str = " + ".join(terms) if terms else "0"
    pGW_str = f"{Cgw:.8e}*({sum_str})"



    md.add_macro("pGW(g)", pGW_str)
    md.add_nonlinear_term(mim_contact2,
                          "Interpolate_filter(raytrace, -pGW(gap2)*(Interpolate(Test_u,raytrace)-Test_u2).normal2, 1)",
                           3)

I would appreciate any feedback, in particular regarding:

• the consistency of this discrete approach within the GetFEM framework,
• and possible numerical or structural improvements to the implementation.

Thank you in advance for your comments.

Best regards,
René

[Konstantinos.Poulios]

Well, you do not need to couple the rough contact analysis model and the finite element model.

The rough contact analysis is done offline, it is always the same as long as you have the same roughness on the entire surface. The outcome for rough contact analysis, whether it is GW or Ciavarella or whaterver is a simple function in the form stress(gap), fitted to the results of your offline rough contact analysis.
Then you add

md.add_macro("dPsi_dg(gap)", "....some_short_and_simple_function_of_gap....")
md.add_nonlinear_term(mim_contact2,"...dPsi_dg(gap)*(Interpolate(Test_u,raytrace)-Test_u2)..."

I do not think it is a good idea to put an entire rough contact model inside the dPsi_dg macro. Such an online computation of rough contact is not necessary and GWFL is not meant for so complicated expressions.

[Djoumessi]

Hi Konstantinos,

Thank you for your analysis. Indeed, I have tried the fitting approach you mention, but I was concerned that the fitted model might not capture all the details and could introduce additional errors.

Here is the macro I implemented for this approach, which seems to work well so far:

md.add_initialized_data("Kgw", float(K));
md.add_initialized_data("b1gw", float(b1));
md.add_initialized_data("b2gw", float(b2));
md.add_macro("pGW(g)", "Kgw*pow(neg_part(g),1.5)*(1 + b1gw*neg_part(g) + b2gw*sqr(neg_part(g)))");
md.add_nonlinear_term(mim_contact2,
                      "Interpolate_filter(raytrace, -pGW(gap2)*(Interpolate(Test_u,raytrace)-Test_u2).normal2, 1)",
                      3);

[Konstantinos.Poulios]

@Djoumessi

Thank you for your analysis. Indeed, I have tried the fitting approach you mention, but I was concerned that the fitted model might not capture all the details and could introduce additional errors.

don’t worry, the modeling error in all these models compared to reality is much larger than the fitting error if you do a reasonable fitting.

[Konstantinos.Poulios]

@Djoumessi

md.add_macro("pGW(g)", "Kgw*pow(neg_part(g),1.5)*(1 + b1gw*neg_part(g) + b2gw*sqr(neg_part(g)))");

this function does not seem to agree with the slope of the red curve here

probably you do not use a consistent definition of the gap when doing the GW analysis. The gap should always be positive, and the contact stress should go to infinite when gap approaches 0.

[Konstantinos.Poulios]

unfortunately I cannot point you to a paper describing this in detail, but the only consistent way to define a gap function is as a difference in far field displacements between a rough and a smooth model. See if you can decrypt the following illustration.

[Djoumessi]

Thank you very much for the illustration. I will try another approach following this guideline and will get back to you with the results I obtain.

Best regards,

[Djoumessi]

Hi Konstantinos,

I tried the fitting-based approach.

I generated a fractal rough surface, then performed an offline fitting to obtain an analytical expression for the contact law. The fitted parameters are then injected into the GetFEM model. I also took inspiration from papers by Peter Wriggers.

Here is the relevant part of the implementation:

K_fit, lam_fit = offline()
print("K_fit =", K_fit, "lam_fit =", lam_fit)

md.add_initialized_data("Kgw", float(K_fit))
md.add_initialized_data("lamgw", float(lam_fit))

# g = local gap (in practice neg_part(gap2)), p(g) = Wg(g)
md.add_macro(
    "pGW(g)",
    "Kgw*pow(neg_part(g),1.5)*(1 - exp(-neg_part(g)/lamgw))"
)

md.add_nonlinear_term(
    mim_contact2,
    "Interpolate_filter(raytrace, "
    "-pGW(gap2)*(Interpolate(Test_u,raytrace)-Test_u2).normal2, 1)",
    3
)

This part seems to work correctly, although I still need to further calibrate the parameters.

However, I am now trying to evaluate and export the contact pressure on surface 3 in order to visualize it. I attempted to do:

mfmult2 = gf.MeshFem(m2, 1)
#mfmult2.set_classical_fem(0)
mfmult2.set_fem(gf.Fem("FEM_QK(3,1)"))
contact_pressure = md.interpolation("pGW(gap2)", mfmult2, 3)

but I obtain the following error:

RuntimeError: (Getfem::InterfaceError)
Error in getfem_generic_assembly_compile_and_exec.cc, line 1294
virtual int getfem::ga_instruction_copy_interpolated_small_vect::exec():
Invalid vector size: 3 != 0

I am not sure why this happens.
Do you have any idea what might be causing this issue or how I should correctly interpolate/export this quantity?

Thank you in advance for your help.

Best regards,
René

[Konstantinos.Poulios]

@Djoumessi

Invalid vector size: 3 != 0

It seems like a small bug, probably the normal vector (Normal) used in the definition of the gap function is not precomputed on nodes, only on Gauss points. The thing is, the normal vector is not well defined on nodes, depending on which of the neighbor elements is used to compute the normal, a different vector will be obtained. In that case, some adhoc rule like computing the average would be necessary.

It is in general better to use a finite element space projection than a pointwise interpolation

mfout_ = gf.MeshFem(m2, 1)
mfout_.set_classical_fem(1)
mfout = gf.MeshFem("partial", mfout_, RG_CONTACT)
PRESS = gf.asm_generic(mim, 1, "pGW(gap2)*Test_t", R_CONTACT, md, 
                       "t", True, mfout, np.zeros(mfout.nbdof()), "select_output", "t")
mass_mat = gf.asm_mass_matrix(mim, mfout, mfout, RG_CONTACT)
PRESS = gf.linsolve_mumps(mass_mat, PRESS).T