# Bio-inspired generative design of contact interfaces subjected to time-dependent loading conditions

David Hernandez-Aristizabal, Santiago Arroyave-Tobón, Kalenia Marquez-Florez, Jean-Marc Linares

► To cite this version:

David Hernandez-Aristizabal, Santiago Arroyave-Tobón, Kalenia Marquez-Florez, Jean-Marc Linares. Bio-inspired generative design of contact interfaces subjected to time-dependent loading conditions. Results in engineering, 2025, 25, pp.104340. 10.1016/j.rineng.2025.104340. hal-04955372

HAL Id: hal-04955372

https://hal.science/hal-04955372v1

Submitted on 18 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Bio-inspired generative design of contact interfaces subjected to
time-dependent loading conditions

David Hernandez-Aristizabal<sup>a,*</sup>, Santiago Arroyave-Tobón<sup>a</sup>, Kalenia Marquez-Florez<sup>a</sup>, Jean-Marc Linares<sup>a</sup>

<sup>a</sup>Aix-Marseille Univ, CNRS, ISM, Marseille, France

**Abstract**

In both biological and human-made mechanical systems, contact interfaces play a crucial role in transferring power. Although their lifespan is enhanced by uniform contact pressure distribution, designing shapes achieving this remains challenging, especially under time-dependent loads. Contrarily, synovial joints are adapted to withstand severe and time-dependent loads. This adaptation is driven by two key principles: shear stress promotes ossification, while hydrostatic pressure promotes proliferation. This study analyses the response of a generative design algorithm for contact problems based on these principles embedded into a growth rule. By varying material stiffness, initial conditions, and loading scenarios, the algorithm successfully adapts contact interfaces for both time-independent and time-dependent loads. The results demonstrate that the growth principles lead to almost uniform contact pressure under time-independent loads, up to 34.7% more uniform than reference shapes. Under time-dependent loads, it generates shapes with reduced maximum shear stress and contact pressure (up to 73%) with respect to initial shapes. Unlike other methods, such as optimisation techniques, this approach does not seek a specific contact interface where pressure is well-distributed. Instead, the development of well-distributed pressure is an intrinsic part of the growth function. Hence, the algorithm provides a family of contact interfaces exhibiting uniform contact pressure whose shape depends on the development time. Additionally, the growth rule is easier to compute than gradients for gradient-based optimisation, and, in contrast to heuristic optimisation techniques, only one direction of change is necessary. Thus, this bio-inspired approach provides a viable alternative for shape optimisation problems in contact mechanics.

**Keywords:** Bio-inspired design, Generative design, Contact interface, Bone synovial joint formation, Shape optimisation, Pressure distribution
20

# 1. Introduction

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

In biological as well as in human-made mechanical systems, contact interfaces are used to transfer power between components [1]. A crucial indicator of their performance is wear resistance. In general, lifespan performance increases as wear volume reduces or as wear distribution becomes more uniform for a given load and contact area. The latter depends on the spatial distribution of (i) material and surface properties [2], (ii) relative speed between contact surfaces, (iii) lubrication conditions, and (iv) contact pressure [3]. In bodies in translation, such as engine pistons, the relative speed is space-independent [4, 5], while, in rotating bodies, it is proportional to the radius. Yet, when there is no radius change, or when it is negligible, such as in most journal bearings [6, 7], the relative speed can be considered space-independent. In addition, contact fatigue lifespan is also related to contact pressure [8, 9], among other factors such as material strength and hardness and environmental conditions [10, 11]. Thus, wear and fatigue performance in contact interfaces are highly related to pressure uniformity, pointing out the need for designing shapes that provide uniform contact pressure for a given loading condition. However, designing such shapes is not an easy task; it is even harder when contact shapes are subjected to complex time-dependent loads, commonly the case in mechanical joints.

36

37

38

39

40

41

42

43

The first definition of a contact interface providing uniform pressure was established for bodies in contact against elastic half-spaces by Lundberg [12]. In [12], a logarithmic expression was derived for crowning cylindrical rollers to achieve uniform contact pressure—crowning involves rounding the contact profile so that it is no longer a straight line but instead follows a curve. Yet, this expression goes to infinity at roller edges. Later, Johns and Gohar [13] modified the expression to avoid this issue and introduced a method to approximate the logarithmic shape using two circles. Therefore, the modified profile is easier to manufacture. Further modifications to Lundberg's logarithmic expression can be found in [14, 15], where genetic algorithms were implemented to find

<footnote>Declarations of interest: none</footnote>

<footnote>*Corresponding author</footnote>

<footnote>E-mail: david.hernandez-aristizabal@univ-amu.fr</footnote>

<footnote>Address: 413 AVE Gaston Berger - 13625 Aix-en-Provence, France</footnote>
44

optimal parameters, and by Wu et al. [16, 17], who also modified it to account for misalignment.

45

Furthermore, crowning can also be seen as an optimisation of the initial gap between two contact surfaces [18]. Applications can be found in [19], where an artificial neural network is applied to find the optimal gap; in [20], where the gap is optimised using evolutionary optimisation; in [21], where a two-circle profile is optimised using particle swarm optimisation; and in [22], where the profile is defined by a super-ellipse. Even though crowning has found several industrial applications (e.g., roller bearings, gear teeth and railways), it is limited to time-independent loads applied in the normal direction to the contact surface.

46

47

48

49

50

51

52

Another approach to the problem of designing contact interfaces is shape optimisation. A general formulation usually reads: find $D^* \in \mathcal{D}$ such that $\Pi(\phi(D^*)) \le \Pi(\phi(D))$ for all $D \in \mathcal{D}$ such that $\phi(D) \in \mathcal{C}$; where $D^*$ is the optimal shape; $\Pi$ is a cost function; $\mathcal{D}$ is a given set of acceptable domains; and $\mathcal{C}$ is a given set of acceptable state variables (displacements, stress, energy and so on) [23]. For the case of contact pressure distribution, a reasonable cost function would be a measure of the contact pressure, e.g. its maximum value. However, the contact pressure is not differentiable with respect to changes of shape on the contact boundary [18, 24–27]. This is an important issue for gradient-based optimisation since it requires the differentiation of the cost function with respect to the design variables, in this case, the shape. In response to this, Klarbring and Haslinger [24] demonstrated that, considering an elastic body in contact against a rigid surface, the minimisation of the potential energy leads to almost uniform contact pressure. Moreover, Klarbring and Rönnqvist [28] simply ignored the non-differentiability for the optimisation of a truss structure with contact conditions. They applied a standard optimisation method (namely, the method of moving asymptotes) and found that, in some cases, a solution was found. They argued that the non-differentiability is implicitly avoided after the discretisation of the model since the non-smooth points are omitted. Another way to circumvent the non-differentiability problem was introduced by Páczelt and Szabó [29], who used a controller function to attain the optimal state of the objective function. Páczelt and Mróz [30, 31] also addressed the problem of shape optimisation by wear, modifying the objective function to minimise wear rate or friction or wear dissipation powers. We refer the reader to [18, 23] for a review of shape optimisation and to [26, 27] for a mathematical description of the problem of shape optimisation with contact. Although shape

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72
73

74

75

optimisation is more robust than crowning techniques, it is way more challenging, and, to the best of our knowledge, it is still limited to optimisation of an unknown contact surface with respect to a known one.

76

77

78

79

80

81

There are other techniques to optimise contact pressure that do not target the contact shape. Instead, they modify stiffness around the contact zone. Hilding et al. [32] presented a shape optimisation method where contact boundaries were independent of the design variables—in other words, only shape changes outside the contact zone were allowed. Further, instead of modifying shape, other alternatives include the optimisation of topology [33–40] and of material properties [41–46].

82

83

84

85

86

87

88

Nonetheless, the works mentioned above have not dealt with the case of time-dependent loads. They focus on finding an appropriate shape for a given load; yet, a change in load direction or location can rapidly lead to uneven distribution of the contact pressure, regardless of the optimisation technique or target. In practice, engineering joints such as bearings [47–50], clutches [9], and disc brakes [51], frequently experience time-varying loads and shifts in the contact region. Therefore, developing methodologies that account for time-dependent loads remains an open challenge in current research and a practical requirement in engineering applications.

89

90

91

92

93

94

95

96

97

98

100

101

The above-mentioned strategies arise from well-established mathematical concepts. Nevertheless, other possibilities are emerging from bio-inspiration to respond to the problem of shape optimisation. Inspiration from nature is particularly interesting since biological joints are subjected to time-dependent loads. A first initiation into bio-inspiration would consist of mimicking structures present in biology. This step is made after the assumption that fully developed biological structures have been functionally adapted through evolution. For example, Sysaykeo et al. [52, 53], Sanz-Idirin et al. [54] developed journal bearings inspired by the form of the mammal elbows. Subsequently, they evaluated their performance in terms of wear and load bearing. Their results showed novel shapes of contact interfaces with competitive wear performance with respect to traditional designs. A proceeding step in bio-inspiration would mimic the process of formation— in other words, it would use biological principles of morphological development. This step is based on the assumption that biological processes of formation have been finely tuned through evolution. In [55–57], a concept of structural shape optimisation based on biological growth was
studied. It mimics the mechanism of growth present in trees and aims to provide uniform stress
beneath the body surfaces. Another example of this bio-inspiration step is shown in [58, 59] where
the principles of bone adaptation were used in topology optimisation.

A promising candidate of bio-inspiration to respond to the problem of shape optimisation with
contact conditions is the process of formation of synovial joints. It is not only composed of genetic
information but also of a strategy of mechanical adaptation to the loading environment [60]. These
joints comprise bone ends, covered with articular cartilage, and a synovial capsule, mediating load
transfer and lubrication. In fact, recent works [61–63] have explored the generation of contact
interfaces based on the modelling of synovial joint formation presented in [64–68]. In [65, 67],
shapes are obtained by a process of stress adaptation, and load transfer is carried out via an inter-
mediary elastic material, simulating the synovial capsule. Applying this model, in [61], a study
of the evolution of maximum shear stress and contact pressure was presented. The results showed
a significant reduction in these values and proved the feasibility of the concept. Subsequently,
the space of parameters of the model was analysed in [62], founding optimal model parameters.
Lastly, in [63], the capsule was removed and the model was adapted to generative design under
frictionless contact. Yet, the results were limited to the generation of one of the contact surfaces
with respect to a known surface under time-independent loads. Thus, two features of the syn-
ovial joint formation were omitted: simultaneous generation of two or more contacting bodies and
adaptation to time-dependent loads. In fact, this limitation (designing only one geometry under
time-independent loads) is also present in most mathematical strategies previously mentioned [12–
17, 19–22, 24, 26–28, 33–46]. Only a few works have considered the case of optimising two parts
in contact [37, 38, 42, 44], and they are restricted to topology and material stiffness optimisation.
Moreover, to the best of our knowledge, there are no works dealing with optimisation of contact
interfaces under time-dependent loads, a common scenario in engineering.

In conclusion, the non-linearity of the contact problem results in discontinuities in gradient-
based optimisation and high computational demands in heuristic optimisation. Consequently, cur-
rent tools for designing the geometry of contact interfaces are limited to time-independent loads
and the optimisation of only one component. In contrast, synovial joints can withstand complex,
time-dependent loads, and their shapes result from a process of simultaneous formation that has
been finely tuned by evolution. Although previous works have demonstrated the suitability of applying these principles to the design of engineering contact bodies, they have not addressed multi-body formation under time-dependent loads yet.

The purpose of the present work is therefore to fill this gap. We aim to generate contact pairs adapted not only to time-independent but also to time-dependent eccentric loads. We start by studying the simultaneous generative design of two bodies in contact under a time-independent load. For this, we propose a parametric study of the generation under realistic loads—ranging from those seen in engineering to those present in synovial joints—, and we investigate the impact of parameters on pressure distribution as well as on geometric characteristics of the resulting surfaces. Then, we study the generation in cases under time-dependent eccentric loads, analysing the obtained pressure distribution along load cycles. Finally, we evaluate the performance of the generated contact interfaces in a three-dimensional framework to demonstrate the applicability of the methodology.

The rest of the paper is organised as follows. In the next section (Section 2), biological observations related to the formation of the bone-joint structure are introduced. Then, the model of bio-inspired generative design and the metrics of analysis are developed in Section 3. In Section 4, numerical examples of implementation of the bio-inspired algorithm are presented. A discussion of these example results is developed in Section 5, and the conclusions of this work are exposed in Section 6.

## 2. Biological model for generative design: bone-joint formation

In this section, we overview the formation of the bone-joint structure, emphasising the role of mechanical stimuli. We then present the state of the art related to modelling the mechanical influence on development. The insights here exposed form the basis of the generative design model discussed in the subsequent section (Section 3).

A key component of interest in the bone-joint process of formation is that mechanical stimuli have a crucial impact. They participate in the regulation of biochemical signals through mechanotransduction [60, 69–74]—a process where mechanical stimuli are translated into biochemical processes. In turn, these biochemical processes are responsible for morphological and structural
changes. Although the direct effects of mechanics on biochemical dynamics are not clear [74], it has been demonstrated that the absence of mechanical loads leads to abnormal bone and joint morphologies [69, 75]. In addition, the same process of formation is present in a large variety of bones and conserved through a broad diversity of vertebrates, suggesting that joint morphologies are inherently linked to specific movements [72]. In fact, Shwartz et al. [70] argued that mechanical regulation has probably enhanced skeletal adaptation during evolution since it allows the creation of joints and movements, and the deletion of those that become unnecessary.

At early developmental stages, the bone-joint structure consists of a cartilage template composed of chondrocytes. The cartilage template then splits and undergoes ossification, giving rise to the bone and the synovial joint [76, 77]. It has been seen that during this transition stress regulates chondrocyte activity, including proliferation, hypertrophy, and extracellular matrix production, as well as cartilage rearrangement, vascularisation, and bone circumferential width [70, 72].

The first mechanical stimuli are thought to coincide with the first muscle contractions, which appear at the same time as the primary ossification centre [75, 78]. They are necessary for the development of the whole musculoskeletal system [72] by apparently regulating cell fate [69, 70, 72]. In other words, they may regulate the differentiation programs that cells commit to (e.g. chondrogenesis), cell genetic expression (e.g. collagen II and X) and cell behaviour (e.g. proliferation, migration and hypertrophy). During this regulation, integrins and cadherins play central roles by mediating cellular responses to mechanical stimuli [79]. These molecules influence the β-catenin and the Ihh-PTHrP signalling pathways [79], which have been found to be important mediators during bone and joint formation. The β-catenin pathway is essential in cell fate maintenance [72], and the Ihh-PTHrP feedback loop has a strong influence on the either proliferative or hypertrophic activity of chondrocytes [70].

Interestingly, patterns of cyclic shear and hydrostatic pressure correlate with developmental and histological patterns [75, 78]. Carter et al. [75, 78] made evident a relation of ossification patterns and patterns of both low hydrostatic pressure and high shear stress using finite element analysis. Consequently, they suggested that low hydrostatic pressure prevents ossification while high shear stress accelerates it. They also proposed that the rate of growth due to proliferation was proportional to hydrostatic pressure, providing cartilage with a self-regulated shaping mechanism
[75, 78]. Ogura et al. [80] found that shear stress has a catabolic effect on the cartilaginous extracellular matrix and stimulates the production of bone matrix proteins. This stimulation induces hypertrophy, which subsequently produces bone elongation [76, 77] and triggers the invasion of osteoclasts, leading to cartilage ossification [76]. The change in material properties due to ossification may further concentrate stresses locally; nevertheless, the mechanical resistance of the tissue improves [81, 82], enabling it to withstand higher stresses. In contrast, hydrostatic pressure favours the maintenance of chondrocytes, enhancing the production of cartilaginous extracellular matrix proteins and facilitating the flow of nutrients [68]. Thus, ossification is influenced by an interplay between these promoting and inhibitory factors. This balance has been formalised in the concept of the osteogenic index [60], which measures the likeliness of a chondrocyte to undergo hypertrophy and lately trigger ossification.

This hypothesis on how hydrostatic pressure and shear stress influence the biochemical dynamics of the developing bone-joint structure has been tested in computational frameworks of growing domains [64–66]. In these works, the authors proposed different growth rules considering mechanical and biochemical stimuli, evaluating the ability to predict shape and ossification patterns. Heegaard et al. [64] started by introducing a growth force proportional to hydrostatic pressure and to chondrocyte density, higher near the articular surface. They used direct contact to transfer load from one developing joint to another. Their results showed that the hydrostatic-pressure-dependent term was necessary to obtain shapes congruent with biological observations. Giorgi et al. [65, 66] continued the exploration of this growth rule and studied the impact of loading conditions. While load transfer was modelled differently, using direct contact in [66] and an elastic capsule in [65], the results of both works showed shape congruency for several loading conditions, demonstrating the ability of the growth function to generate a vast diversity of shapes. Moreover, Carrera-Pinzón et al. [67] considered the case where the stress-independent growth term was neglected using an elastic capsule for load transfer. Their results were also congruent, indicating the pivotal influence of mechanical regulation. Nonetheless, Takahashi et al. [83] found from experimental observations that the effect of hydrostatic pressure on growth disappeared when hydrostatic pressure was too high. Regarding shear stress, Carter and collaborators considered that it accelerated endochondral ossification due to its correlation with ossification patterns.
Carrera-Pinzón et al. [67] also calculated the osteogenic index to predict ossification patterns. It increases with shear stress and decreases with hydrostatic pressure. Their results showed good consistency with biological joints; both ossification patterns and pathological shapes were predicted. Lastly, Sadeghian et al. [68] also included ossification by locally modifying material properties. Therefore, they also predicted the evolution in time of ossification patterns in a three-dimensional framework. In addition, they assumed a Hertzian pressure distribution between the developing joints to model load transfer.

Table 1 summarises the characteristics of these models. Notably, while most models consider a stress-independent growth term, shape congruency is primarily linked to the hydrostatic-pressure-dependent growth term. Additionally, ignoring ossification effects on material properties does not seem to affect shape congruency. Lastly, the growth rule can be applied in conjunction with various load transfer models, including direct contact, elastic capsules, or Hertzian pressure.

In conclusion, these computational models [64–68] demonstrate the accuracy of the principles of formation to reproduce bone-joint shapes. These principles read: (i) hydrostatic pressure upregulates chondrocyte proliferation, leading to the growth of the cartilage template, and (ii) shear stress upregulates chondrocyte hypertrophy, triggering local ossification, and hence inhibiting proliferation. Moreover, these models demonstrate that these principles are applicable in several modelling frameworks and loading environments. These works are therefore the foundation of our model for generative design. Particularly, we will consider an isotropic growth force proportional to hydrostatic pressure, present in all the computational models [64–68]; we will neglect the synovial capsule considering direct contact, as in [64, 66]; and we will not consider material stiffening due to ossification, as in [64–67].

### 3. Mathematical translation of bone-joint growth for generative design

In this section, we present the model of shape generation in mathematical terms. We start with a general model of growth for biological tissues with contact conditions. Later, we simplify it to minimise the number of parameters. This approach (from a general to a particular model) has two advantages: (i) it allows us to explicitly point out the assumptions made here, and (ii) it establishes a framework to envisage future works.
### 3.1. General growth model with contact interfaces

Let $\Omega_i(t) \subset \mathbb{R}^3$ be a set of $k \in \mathbb{N}$ evolving domains; where $1 \le i \le k$ and $t$ belongs to the time interval $T := [0, t_f]$, with $t_f > 0$ being a time limit. Let us divide the contour of each domain, $\Gamma_i := \partial\Omega_i$, into three separated regions: $\Gamma_i^d$, where the displacements are known, $\Gamma_i^f$, where the external loads are known, and $\Gamma_i^c$, where both displacements and loads are unknown. In biological growth, it is common to consider that inertial forces are negligible [84–86]; hence, given $\Omega_i^0 := \Omega_i(t = 0)$, the problem of growth [84–86] with frictionless contact boundary conditions [27] reads: find the displacement field $\mathbf{u}_i(\mathbf{x}(t), t) : \Omega_i(t) \times T \to \mathbb{R}^3$, such that:

$$ \nabla \cdot (\boldsymbol{\sigma}_i + \boldsymbol{\sigma}_g) = \mathbf{0} \qquad \text{in } \Omega_i(t), \qquad (1) $$

$$ \boldsymbol{\sigma}_i \cdot \mathbf{n} = \mathbf{t}_i \qquad \text{on } \Gamma_i^f(t), \qquad (2) $$

$$ \mathbf{u}_i = \mathbf{u}_i^0 \qquad \text{on } \Gamma_i^d(t), \qquad (3) $$

$$ \mathbf{u}_i \cdot \mathbf{n} - \delta \le 0 \qquad \text{on } \Gamma_i^c(t), \qquad (4a) $$

$$ (\boldsymbol{\sigma}_i \cdot \mathbf{n}) \cdot \mathbf{n} = -p_c \le 0 \qquad \text{on } \Gamma_i^c(t), \qquad (4b) $$

$$ (\mathbf{u}_i \cdot \mathbf{n} - \delta) p_c = 0 \qquad \text{on } \Gamma_i^c(t), \qquad (4c) $$

where $\boldsymbol{\sigma}$ is the stress tensor which is related to material properties; $\boldsymbol{\sigma}_g$ is the growth tensor which depends on mechanical and biochemical dynamics; $\mathbf{n}$ is the outward-pointing unit normal vector; $\mathbf{t}$ is the known external traction; $p_c$ is the magnitude of the contact pressure; $\delta$ is the initial gap between contacting surfaces (the distance separating the surfaces before the application of loads); and $\mathbf{u}^0$ are the known contour displacements. See Fig. 1 for a schematic representation of the problem.

The dynamics of biochemical species can be modelled by reaction-diffusion systems on growing domains (see for example [87]). However, since our current focus lies in examining the influence of mechanical stimuli on shape, we refrain from explicit modelling of multiple biochemical species, as it falls outside the scope of this work. Instead, we opt to consider a single biochemical species with a known distribution indicating the active growth regions.
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of the problem. The purple zones illustrate the disposition to proliferate, $a_1$, defined by $d_{lim}$ using the Sigmoid function Eq. (8).</figcaption></figure>

## 3.2. Definition of the growth function

In synovial joint formation, growth depends on biochemical and mechanical factors. The former are associated with genetic information, while the latter with loading environment. This can be represented as:

$$ \sigma_g = f_{b_1}f_m + f_{b_2}, \qquad (5) $$

where $f_{b_1}$, $f_m$, and $f_{b_2}$, represent mechanical and biochemical growth. $f_{b_1}$ represents cell disposition to proliferate due to mechanical stimuli, and $f_m$ and $f_{b_2}$ are respectively the stress-dependent and -independent growth tensors. $f_{b_1}$ and $f_{b_2}$ depend on an $n_s$-dimension vector of biochemical species, say $\mathbf{a}$, where $n_s$ is the number of such species. In computational models of formation, $\sigma_g$ is commonly considered isotropic [64–67], then:

$$ f_m = g(\sigma) \mathbf{I} \qquad (6) $$

where $\mathbf{I}$ is the identity tensor and $g$ is the function that comprises the growth rules—the adaptation to the loading environment.

As mentioned above, the adaptation process is led by two principles: (i) cell proliferation is upregulated by hydrostatic pressure, while (ii) ossification is upregulated by high shear stress. In terms of shape, these rules can be translated as: growth is promoted by hydrostatic pressure and inhibited by high shear stress. In a previous work [63], we demonstrated that, for generative
279 design without a capsule, a third rule is necessary to produce uniform contact pressure. This
280 new rule indicates that too high hydrostatic pressure also inhibits growth. This is consistent with
281 biological observations indicating that after a certain value of hydrostatic pressure, proliferation is
282 no longer promoted [83]. Let $\sigma_{hyc}$ be the hydrostatic pressure and $\sigma_v$ be the von Mises stress—we
283 can use it as a measure of shear stress.

284 For the construction of the growth rule, we need the following normalised quantities. Let $\sigma_{hyc}^*$ and $\sigma_v^*$ respectively be the normalised hydrostatic pressure and shear stress defined as:

$$ \sigma_{hyc}^* = \frac{\sigma_{hyc}}{\max(\sigma_{hyc}|_t)}, \qquad \sigma_v^* = \frac{\sigma_v}{\max(\sigma_v|_t)}. \qquad (7) $$

286 Let $S$ be the Sigmoid function, which is a smooth equivalent of the Heaviside function:

$$ S(\varphi) = (1 + \exp(\vartheta\varphi))^{-1}, \qquad (8) $$

287 where $\vartheta$ is a positive parameter that controls smoothness. Then, we can write the growth rule as:

$$ g(\boldsymbol{\sigma}) = \sigma_{hyc}^* S(\sigma_v^* - \tau_{\lim}^*) S(\sigma_{hyc}^* - \sigma_{\lim}^*), \qquad (9) $$

288 where $\tau_{\lim}^*$ controls the growth inhibition due to shear stress, $g \to 0$ as $\sigma_v^* > \tau_{\lim}^*$. Similarly, $\sigma_{\lim}^*$ controls the inhibition when hydrostatic pressure is too high, $g \to 0$ as $\sigma_{hyc}^* > \sigma_{\lim}^*$.

290 For the definition of the biochemical terms, we shall consider only one biochemical species, 291 say $a_1$, of known distribution. This is similar to what was done in [62] where it was used to 292 determine the active zones of growth—zones where growth forces are applied. We shall also 293 consider that such distribution only depends on the distance to the contact boundary. This is also a 294 common assumption in computational biology where $a_1$ represents chondrocyte density [64–68]. 295 We can hence express $a_1$ as:

$$ a_1 = S(d_{\Gamma_i^c}(\mathbf{x}(t)) - d_{\lim}), \qquad (10) $$

296 where $d_{\Gamma_i^c}(\mathbf{x}(t))$ is the minimum distance from $\mathbf{x}(t)$ to $\Gamma_i^c(t)$ and $d_{\lim}$ is a distance threshold, see 297 Fig. 1. $a_1 \to 0$ as $d_{\Gamma_i^c}(\mathbf{x}(t)) > d_{\lim}$ and $a_1 \to 1$ as $d_{\Gamma_i^c}(\mathbf{x}(t)) < d_{\lim}$. In other words, $a_1$ is 1 at $\mathbf{x}$ 298 closer to $\Gamma^c$ and 0 at $\mathbf{x}$ far from $\Gamma^c$. Therefore, $d_{\lim}$ defines the size of the active zone as seen in 299 Fig. 1.
Once we have determined $a_1$, we can choose $f_{b_1}$ to be linear, as has been done in [64–67]:

$$f_{b_1} = \alpha_g a_1, \qquad (11)$$

where $\alpha_g$ is a positive growth rate parameter. In the present work, $\alpha_g$ is constant and small enough to avoid abrupt shape changes, according to a convergence analysis realised in a previous work [63]. Yet, to further regulate the growth and ensure more controlled shape changes, $\alpha_g$ could be dependent on several state variables adjusting the growth strength. For instance, it could depend on the rate of area increase to avoid excessive or abrupt expansions. Similarly, the rate of curvature or shape change could influence $\alpha_g$. Such a formulation would also improve the computational efficiency of the algorithm by reducing the number of iterations needed while maintaining smooth evolution.

Heegaard et al. [64] numerically tested two growth conditions. In the first, growth is only dependent on biochemical factors ($f_{b_2}$ in Eq. (5)); and in the second, it depends on both biochemical and mechanical stimuli (both $f_{b_1}$, $f_m$ and $f_{b_2}$). They concluded that both of them were necessary to obtain congruent bone shapes. However, they did not consider the case where only mechanical stimuli drive growth ($f_{b_2} = 0$). This case was later tested by Carrera-Pinzón et al. [67] who showed that mechanical stimuli are sufficient to develop congruent shapes. Therefore, for this work, we neglect $f_{b_2}$ as was done in [62, 67].

### 3.3. Material mechanical behaviour

So far, we have not made any assumptions related to material behaviour. However, to solve Eqs. (1) to (4), we need to establish a specific stress-strain relation. Cartilage and bone can be modelled as materials with poro-visco-elastic properties [88, 89]. Yet, most of the computational models that have been proposed to model bone and synovial joint formation are limited to linear elasticity [64–68]. Moreover, engineering materials are likely to work in this range. Therefore, for the sake of simplicity, we can consider the linear elastic stress-strain relation. In mathematical terms:

$$\sigma_i = \frac{E_i}{1 + \nu_i} \left( \varepsilon_i + \frac{\nu_i}{1 - 2\nu_i} \text{tr}(\varepsilon_i) \right), \qquad (12)$$
where $v$ is the Poisson coefficient; $E$ is the Young's modulus; and $\varepsilon$ is the linear strain tensor. Both $v$ and $E$ can depend on the biochemical species, time and space. Such dependences would reproduce the effects of bone ossification and cartilage maturation. Nonetheless, in this work, we will assume that both of them are constant, as previously done in [61–67].

For possible extensions towards more biological behaviours, we refer the interested reader to [85, 90] for growth models applied to bone regeneration and cell migration where the material is considered poro- and visco-elastic, respectively.

### 3.4. The computational framework

For the approximation of Eqs. (1) to (4), we used finite elements with code_aster [91] and Gmsh [92] as a meshing and remeshing tool. Due to the cyclic nature of physiological loads, the model is divided into two steps, as was done in [61–63, 67, 68]. The first step represents a loaded state where both external tractions and contact loads are activated while the second one represents a released state where neither of them is present. In addition, during the first step, it is usually assumed that the growth stress is negligible with respect to external loads. A schematic representation of the algorithm is shown in Fig. 2, and the implementation is available online in: https://github.com/sarroyavet/BioDesign_ANR/tree/main/GenerativeDesign. Fig. 2 also contains the divided system of equations where $\mathbf{u}^1$ and $\mathbf{u}^2$ are the displacement fields for the loaded and the released states, respectively. As can be seen, each step requires the solution of a nonlinear (the contact step) and a linear (the growth step) problem. Nevertheless, it does not require sensitivity analyses—which is the main issue of gradient-dependent shape and topology optimisations—nor does it need multiple initial geometries—which is the case of heuristic optimisations such as those using evolutionary algorithms.

Fig. 2 also indicates the possible need for remeshing. Since we are working with an evolving domain, the mesh will progressively distort. To reduce numerical errors, we limit the maximum length of contact elements, $h_c^{\lim}$, and the minimum element quality, $e_q^{\lim}$. The latter was computed following [93]. It relates the perimeter of an element with that of a reference one with the same area. In linear elements, this quality measure has a maximum value of 1 and a minimum value of 0. For triangles and quadrilaterals, the reference elements are respectively an equilateral triangle
and a square. The element quality is hence given by:

$$e_q = \begin{cases} \left( \frac{6\sqrt{A_e}}{P_e} \right)^2, & \text{for triangles,} \\ \left( \frac{4\sqrt{A_e}}{P_e} \right)^2, & \text{for quadrilaterals,} \end{cases} \qquad (13)$$

where $A_e$ and $P_e$ are respectively the element area and perimeter.

### 3.5. Metrics for the generative design of contact interfaces

In this section, we describe the metrics needed to analyse the response of the model. We will characterise contact interfaces from both geometric and performative points of view. To do that, we will calculate the surface curvature, the conformity between two contacting surfaces and the uniformity of the contact pressure distribution.

For the curvature metric, let us consider the normalised $L_2$-norm of the mean curvature:

$$Q_{\kappa}(t) = \frac{\|\kappa(\mathbf{x}(t))\|_{L_2(\Gamma^{\mathrm{ce}}(t))}}{\sqrt{|l_c(t)|}}, \qquad (14)$$

where $\kappa(\mathbf{x}(t))$ is the mean curvature at $\mathbf{x}$; $\Gamma^{\mathrm{ce}}(t)$ is the effective contact boundary; and $l_c$ is the contact length, the size of $\Gamma^{\mathrm{ce}}$. Notice that $\Gamma^{\mathrm{ce}}$ in general differs from $\Gamma^{\mathrm{c}}$; $\Gamma^{\mathrm{ce}}$ is the part of $\Gamma^{\mathrm{c}}$ where $p_c > 0$.

For the conformity between contacting surfaces, we first need to make some assumptions about the type of contacting surfaces. We shall consider that they are two-dimensional contact interfaces that can be described as: $\Gamma_1^{\mathrm{ce}} := \{(x, y)|y = f_1(x)\}$ and $\Gamma_2^{\mathrm{ce}} := \{(x, y)|y = f_2(x)\}$ where $x \in [a, b]$ (as illustrated in Figs. 3a and 3b). The conformity between $\Gamma_1^{\mathrm{ce}}$ and $\Gamma_2^{\mathrm{ce}}$ can be measured by:

$$Q_c(t) = \frac{1}{\sqrt{a-b}} \left\| 0.5 + \frac{\kappa_1(x)\kappa_2(x)}{\kappa_1^2(x) + \kappa_2^2(x)} \right\|_{L_2([a,b])}, \qquad (15)$$

where $\kappa_1$ and $\kappa_2$ are the mean curvature values of $\Gamma_1^{\mathrm{ce}}$ and $\Gamma_2^{\mathrm{ce}}$, respectively. $Q_c$ goes to: 1 if $\kappa_2$ approaches $\kappa_1$, to 0.5 if either $\kappa_1$ or $\kappa_2$ is zero along $[a, b]$, and to 0 if $\kappa_2$ approaches $-\kappa_1$.

The uniformity of the contact pressure distribution can be calculated as:

$$Q_p(t) = 1 - \frac{\|p_{\mathrm{eq}}(x, t) - p_c(x, t)\|_{L_2(\Gamma^{\mathrm{ce}}(t))}}{\|p_c(x, t)\|_{L_2(\Gamma^{\mathrm{ce}}(t))}}, \qquad (16)$$

where $p_{\mathrm{eq}}$ is the uniform equivalent contact pressure [63]. It is given by:
$$p_{\text{eq}}(x, t) = \begin{cases} \frac{\int_{\Gamma^{\text{ce}}(t)} p_c(x, t) d\Gamma^{\text{ce}}(t)}{l_c(t)} & \text{if } p_c(x, t) > 0, \\ 0 & \text{otherwise.} \end{cases} \quad (17)$$

For the cases under cyclic eccentricity, we will need the fields of maximum, minimum, mean and amplitude of the contact pressure. For a cycle $c_j$, these fields can be defined as:

$$\begin{align} P_c^{\max}(x, c_j) &= \max_{t \in c_j} p_c(x, t), \notag \\ P_c^{\text{mean}}(x, c_j) &= \frac{1}{|c_j|} \int_{c_j} p_c(x, t) dt, \tag{18} \end{align}$$

$$\begin{align} P_c^{\min}(x, c_j) &= \min_{t \in c_j} p_c(x, t), \notag \\ P_c^{\text{ampl}}(x, c_j) &= P_c^{\max}(x, c_j) - P_c^{\min}(x, c_j). \end{align}$$

Lastly, to investigate the impact of some parameters $X_i$ on a specific metric, such as $Q_p$ or $Q_c$, we can consider a metamodel [94]. Let us consider that a metric $y$ can be modelled by two variables $X_1$ and $X_2$. We can study the linear and nonlinear impact of $X_i$ on $Y$ considering the following second-order model:

$$\hat{y}(X_1, X_2) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_{12} X_1 X_2 + \beta_{11} X_1^2 + \beta_{22} X_2^2, \quad (19)$$

where $\hat{y}$ is the approximation of $y$, and $\beta_j$ are regression parameters. For a given set of observation values $w$, Eq. (19) can be rewritten as [95]:

$$\hat{\mathbf{y}} = \mathbf{X}\beta, \quad (20)$$

and the regression parameters $\beta_j$ can be calculated using the Moore-Penrose pseudo-inverse $\mathbf{X}^\dagger$ [96] as:

$$\beta = \mathbf{X}^{\dagger} \mathbf{w} = (\mathbf{X}^{\mathrm{T}} \mathbf{X})^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{w}. \quad (21)$$

The covariance matrix of $\beta$ is given by [95]:

$$\operatorname{cov}(\beta) = (\mathbf{X}^{\mathrm{T}} \mathbf{X})^{-1} \operatorname{var}(\mathbf{r}_{\mathrm{e}}), \quad (22)$$
where var($\mathbf{r}_e$) is the variance of the estimated residuals, $\mathbf{r}_e = \mathbf{w} - \hat{\mathbf{y}}$. Hence, the variation of $\beta$ can be estimated through the Student t-statistic with $n - p$ degrees of freedom as:

$$ \Delta\beta_j = t_{n-q}\sqrt{\text{cov}(\beta_j, \beta_j)}, \qquad (23) $$

where $t_{n-q}$ is the t-value—here, it is considered with a confidence level of 0.95—; and $n$ and $p$ are the number of samples and model parameters, respectively. The model uncertainty can be measured by the root-mean-square error as:

$$ \text{RMSE}(\mathbf{r}_e) = \sqrt{\frac{1}{n-p} \sum_{i} r_{e_i}^2}. \qquad (24) $$

Finally, the accuracy of the metamodel can be measured by the coefficient of determination $r^2$, given by [95]:

$$ r^2 = \frac{\sum_i r_{e_i}^2}{\sum_i (\bar{w} - w_i)^2}, \qquad (25) $$

where $w_i$ are the real observed values used to compute $\beta$ in Eq. (21), $\bar{w}$ the average value, and $y_i$ the approximated values using the metamodel. $r^2 \in [0, 1]$ indicates how well the metamodel explains the variation of the variable of analysis, $y$. For instance, if $r^2 = 0.9$, the metamodel describes 90% of the variation of $y$.

## 4. Results

To explore the response of the proposed generative design algorithm, let us consider the study case illustrated in Fig. 3a. As can be seen, a two-dimensional framework is considered where the initial geometry is a semi-circle in contact with a flat block. The present case is modelled under plane strain assumptions. Since the width, $w$, is unitary in plane mechanics, the initial geometry is defined by the radius, $R$, of the upper body. This configuration has a level of conformity $Q_c = 0.5$ (see Eq. (15)). Fig. 3a also illustrates the loading conditions of the study case: an external force, $F$, with an eccentricity, $e$. Notice that, in Fig. 3a, static equilibrium is only ensured when $e = 0$. Therefore, for $e \neq 0$, the shape must already have been modified, see Figs. 3b and 4. We shall also evaluate the contact pressure distribution of a circular and a flat shape as control groups. Notice that the circular shape can only be used when $e = 0$; otherwise, there is no static equilibrium.
We shall use the parameter values for the mechanical contribution of growth defined in [63]:

$\sigma_{\text{lim}}^* = 0.5$, $\tau_{\text{lim}}^* = 0.5$, $\vartheta = 7.1$ and $\alpha_g/E = 0.00625$. With such parameters in Eq. (9), the study case is characterised by: $|\mathbf{F}|, e, R, E, \nu, d_{\text{lim}}$ and $t_f$. We can reduce it by considering the following dimensionless variables:

$$S_f = \frac{1000|\mathbf{F}|}{Ewl_{cf}}, \quad B_e = \frac{l_{c0}}{l_{cf}}, \quad M_f = \frac{2e}{l_{cf}}, \quad D_g = \frac{d_{\text{lim}}}{l_{cf}}, \quad (26)$$

where $l_{c0}$ and $l_{cf}$ are respectively the initial and the final contact length. $S_f$ is the stiffness factor which relates the pressure at the contact interface with the material rigidity; $B_e$ is the normalised initial size which indicates the amount of growth time required—this time augments as the initial contact length gets smaller with respect to the final one—; $M_f$ is the eccentricity factor relating the eccentricity with the maximum allowed eccentricity—since the contact boundary is responsible for the equilibrium state, we cannot consider an eccentricity higher than $l_{cf}/2$, see Fig. 4—; and $D_g$ is the normalised size of the growth region relating the size of the growth region with the size of the final contact interface. For the definition of the radius of the initial geometry, we can use the theory of Hertz for a cylinder in contact with a half-space [13]. Using the defined dimensionless numbers:

$$\frac{R}{w} = \frac{125\pi l_{cf} B_e^2}{4(1 - \nu)^2 w S_f} \quad (27)$$

Lastly, let us take Poisson's coefficient of steel since it is largely used in engineering contexts; hence, $\nu = 0.3$ hereafter. Therefore, the study case is now determined by $S_f$, $B_e$, $M_f$ and $D_g$.

### 4.1. Response of the model: geometric characteristics and contact pressure distribution

Let us begin by analysing the response of the generative design model when the force is centred and there is no growth distance threshold. In other words, we start by the study of the influence of the stiffness factor ($S_f$) and the normalised initial size ($B_e$) when there is no eccentricity ($M_f = 0$) and the growth region covers the whole domain ($D_g \to \infty$). For that, we can use a metamodel to explicitly describe $Q_p$, $Q_c$, $Q_{\kappa_1}$ and $Q_{\kappa_2}$ as functions of $S_f$ and $B_e$.

To capture linear and nonlinear effects, let us consider a second-order metamodel Eq. (19), with $y = Q_{p,c,\kappa_1,\kappa_2}$, $X_1$ related to $B_e$ and $X_2$ related to $S_f$. Then, $\beta_1, \beta_2, \beta_{12}, \beta_{11}$ and $\beta_{22}$ are related
to $B_e$, $S_f$, $B_e S_f$, $B_e^2$ and $S_f^2$, respectively. To define $\beta = [\beta_0, \beta_1, \beta_2, \beta_{12}, \beta_{11}, \beta_{22}]^T$, we need to fit the model using a set of sampling. In other words, we need to execute the generative design algorithm for a set of $S_f$ and $B_e$ combinations and evaluate the response of the system. To establish the sampling points, we can use a full factorial design for the parameter space considering the minimum, maximum and centred values of each variable [97]. This type of design ensures that $X_1$ and $X_2$ vary independently in the numerical experiments [97]. We first need to define the domains of $S_f$ and $B_e$.

The domain for $S_f$ can be defined from the values observed both in engineering and biology, see Table 2. In engineering cases [47–50, 98], this value varies between 0.319 and 8.183. Yet, it can be as low as 0.001 [7]. For biological cases, the projected pressure in elbows is around 2.4MPa for a wide variety of species [99]. In addition, the Young modulus at a single epiphysis varies in space between 10 and 3000MPa [81], suggesting a stiffness factor varying between 0.8 and 240. However, given the spatial variability of the Young modulus, we can limit the order of $S_f$ to 10. Therefore, let us consider as domain of interest $S_f \in [0.1, 10.0]$ which considers most engineering cases seen in literature [47–50] as well as biological loading conditions [81, 99]. This domain reflects loading conditions where the contact pressure is 10000 to 100 smaller than the material rigidity.

For the normalised initial size, we need to consider the impact on the computational time; the smaller $B_e$ the longer the time required to achieve the desired size. Based on this constrain, let us consider [0.2, 0.8] as the acceptable domain for $B_e$. $B_e > 0.8$ represents cases where the development time is too short; consequently, the final shape is not completely adapted to the loading conditions. In contrast, $B_e < 0.2$ represents cases where the development time is too long which would make the simulations expensive in terms of the computational time.

Both domains can be studied as logarithmic ranges. Therefore, the samples are combinations of $S_f \in \{0.1, 1, 10\}$ and $B_e \in \{0.2, 0.4, 0.8\}$, and are normalised into the reference axis as $X_1, X_2 \in \{-1, 0, 1\}$.

For the domain discretisation, Table 3 shows a convergence analysis of $h_c^{\lim}/w$ with respect to the contact pressure at $t = 0$. The first row of the table (error 1) corresponds to the present case. From this result, we can take $h_c^{\lim}/w = 0.0025$ for the rest of the simulations. For this size, the error
is less than 2.5%. Additionally, for this mesh size, the algorithm took around 38.39 seconds by iteration, with an average number of nodes of 31876; hence, 0.0012 seconds by iteration by node. Moreover, around 76.8% of the iteration time corresponded to the solution of the contact problem, 21.3% to the computation of the growth displacement and 1.9% to other internal computations. The simulations were run in an 8-core 11th Gen Intel ® Core™ i7-1185G7 at 3.00GHz with 32 GB RAM with Ubuntu 20.04.6 LTS.

Fig. 5 shows the results of the metamodels describing the relation among the stiffness factor ($X_1$ related to $S_f$), the normalised initial size ($X_2$ related $B_e$) (both in the scale [-1, 1]), the contact pressure quality ($Q_p$), the conformity between the contacting surfaces ($Q_c$), and the normalised curvature norm ($Q_{\kappa1}$)—the behaviour of $Q_{\kappa2}$ was similar to $Q_{\kappa1}$. Although we can compute $Q_p$, $Q_c$ and $Q_{\kappa1}$ at each time step, we use only the final values to construct the metamodels. In other words, we consider that the final state is the result of the process of formation. The absolute value of the components of $\beta$ indicates the impact of the variable on the response. As can be seen in Fig. 5a, the metamodels were refined removing less significant variables. The uncertainty was computed using Eq. (24) in an iterative fashion after removing the less significant variables. The refinement was stopped once the uncertainty started to augment. Hence, as we can see in Fig. 5a, $Q_p$ depends mainly on the square of $B_e$ and the product of $B_e$ and $S_f$; $Q_c$ is not linearly affected by $S_f$; and $Q_{\kappa1}$ is not affected by the square of $S_f$. The impact of each variable can also be seen in Fig. 5b which is based on the values of $\beta$ in Table 4, the up and down arrows indicate whether the relation is direct or inverse, respectively. As can be seen, $Q_p$ and $Q_c$ are mainly dominated by nonlinear effects while, $Q_{\kappa1}$ is mainly dominated by linear interactions. Further, Table 4 also indicates the coefficient of determination of each metamodel, $r^2$, see Eq. (25): for $Q_p$, $Q_c$ and $Q_{\kappa1}$ $r^2 = 0.721, 0.877$ and $0.977$, respectively. From these values, we can say that the second-order metamodel is appropriate to predict the curvature and the conformity metrics. The combined impact of linear and quadratic effects in the metamodels can more clearly be seen in Fig. 5c, where the reddish surfaces indicate the average behaviour and the gray ones the model variation. Therefore, from Fig. 5c, higher curvatures are obtained with lower $B_e$—which means longer development times—and higher $S_f$—which means higher contact pressure or lower material stiffness—, and from Fig. 5c, the conformity between contacting surfaces increases when the development time is
long and the stiffness factor is small. Lastly, although 27.9% of the behaviour of $Q_p$ is not capture
by the metamodel, the graph indicates that the quality of the contact pressure is mainly affected by
$B_e$ and is maximum around $B_e = 0.5$. This can be explained by the fact that as $B_e$ increases, the
time of development is shorter leading to a lack of development, a sort of immature state. On the
other hand, as $B_e$ decreases the development time augments and the contact pressure quality enters
an oscillatory state. This is shown in Fig. 7c and explained below. Nevertheless, the results of the
metamodel support the choice of the growth parameters as $Q_p$ is only slightly affected by $B_e—Q_p$
varied between 0.93 and 0.97 indicating that the generated shapes have almost uniform distribu-
tions of pressure despite the variations of $S_f$ and $B_e$. In contrast, the Hertzian contact pressure
distribution using a circular shape provides gives $Q_p = 0.72$, and for flat profiles, $Q_p = 0.78$.
Hence, the algorithm generates contact interfaces with a pressure distribution between 19.2% and
34.7% better than the reference shapes.

Fig. 6 shows the generated shapes. They can be classified into three categories: (i) convex $f_1$ with concave $f_2$ (for example, $S_f = 10$ and $B_e = 0.2$), (ii) concave $f_1$ with convex $f_2$ (for example, $S_f = 0.1$ and $B_e = 0.4$), and (iii) both $f_1$ and $f_2$ have concave and convex parts (for example, $S_f = 1$ and $B_e = 0.2$). In addition, since both generated surfaces are curved, they may also withstand horizontal forces.

Fig. 7 depicts the evolution of the generative process. It corresponds to the case with $S_f = 0.1$ and $B_e = 0.4$. In Fig. 7a, we can see the evolution of the contact pressure; it begins with the Hertzian shape, and as the algorithm goes along, it flattens and the contact length augments. This augmentation is seen in Fig. 7b with respect to time. Interestingly, it shows an initial fast increase of the contact length followed by a slower continuation. This behaviour is also present in the other cases where $B_e < 0.8$ and $S_f \le 1$, and the point of change is also that where the contact length doubles the initial length. This change in the rate of length increase is associated with the category of contact pair. In fact, the initial geometry is of category I, and as long as the contact pair remains in this category, the rate of length increase is high. If the contact pair transitions from category I to II or III, the rate of length increase is reduced, as seen in Fig. 7b. Lastly, Fig. 7c shows the variation of the pressure quality. It increases rapidly from $Q_p = 0.72$ to 0.90, and then, it continues to increase slowly with some oscillations until 0.95. In addition, we can also remark that
when the contact length doubles the initial length, the oscillation of the pressure quality increases.
Nevertheless, $Q_p$ remains high, indicating an almost uniform contact pressure.

## 4.2. Response of the system under time-dependent loads

In mechanical engineering, loads are commonly misaligned with respect to the joint centre. For example, in Figs. 3a and 3b, the eccentricity, $e$, represents a misalignment from the centred position. The lifetime and the performance of mechanical elements can be compromised if the contact interface is not adapted to it. For instance, even though roller and journal bearings are usually crowned to ensure an almost uniform contact pressure distribution under centred load, due to load misalignment the resulting distribution is highly concentrated on small regions [16]. In turn, these regions are overexposed leading to premature failure. Moreover, such regions are difficult to predict since they depend on the misalignment magnitude which may vary with time.

As synovial joints are subjected to time-dependent loads, it is reasonable to evaluate the ability of the proposed generative design method to provide shapes adapted to misalignment conditions. In other words, we shall test the response of the system when the eccentricity factor is not zero, $M_f \neq 0$. Nevertheless, we cannot use the same initial geometry as the previous case because it cannot support eccentric loads. Instead, we can utilise a geometry generated when $M_f = 0$, as those shown in Fig. 6. Since in the previous case we identified three shape categories, we can use a generated geometry from each category. For instance, let us consider the shapes resulting from $S_f = 10$ and $B_e = 0.2$, $S_f = 0.1$ and $B_e = 0.4$, and $S_f = 1$ and $B_e = 0.2$ from categories I, II and III, respectively, see Fig. 6. Additionally, to extend the analysis to time-dependent states, we can also consider a varying $M_f$ as: $M_f = M_f(\sin(\omega t))$ with $\omega$ being the oscillation frequency. In the following, the maximum value of $M_f$ is 0.5 which corresponds to an eccentricity of half the distance from the centre to the edge of the contact zone. Although higher eccentricity values can be used, they would lead to a strong decrease in the contact area as this is necessary for static equilibrium, see Fig. 4.

In this case, the convergence analysis for the discretisation is limited to the size of the mesh used in the previous section. The second row of Table 3 (error 2) shows the error of the contact pressure as $h_c^{\text{lim}}/w$ approaches 0.0025. As can be seen, the error is less than 2.5% for a size of
0.02, an order of magnitude larger than in the previous case. In fact, for the same size, 0.0025, the error is lower than 1%. The reduction of the magnitude of the error is associated with the size of the contact length. In fact, since the geometry is larger, there are more elements in contact. Nevertheless, even though we could use a larger size, as indicated by the analysis, we shall keep $h_c^{\lim}/w = 0.0025$. As we will see, the eccentricity leads to a reduction in the contact length; hence, the number of elements in contact is also reduced.

Fig. 8 shows the evolution of the pressure distribution when the generative design is under constant eccentric loads corresponding to $M_f = 0.5$. Regarding the value of $S_f$, it is the same as the one used in the previous generation. Notice that the generated profiles can be scaled by modifying $f_1$ and $f_2$ if $S_f$ is modified [13, 62]. As can be seen in Fig. 8, due to the eccentric load, the pressure distribution at $t_0$ is no longer uniform. For the category-I shape (developed with $S_f = 10$ and $B_e = 0.2$), it is concentrated at the interior of the contact boundary while, for category-II ($S_f = 0.1$ and $B_e = 0.4$) and -III ($S_f = 1$ and $B_e = 0.2$) shapes, it is concentrated at the edges. In addition, such a concentration is more important in shapes of categories II and III. This difference can also be seen in the distribution of the shear stress as shown in the top row of Fig. 9. It illustrates the distribution of the shear stress normalised with respect to the maximum initial value. As can be seen, the reddish region is larger for the category-I shape than for the category-II one—this indicates that the former shape has a better distribution of the shear stress at $t = 0$. Moreover, Figs. 8 and 9 also show a significant reduction of the surface that is effectively in contact, the right-hand-side part. This is in fact an expected behaviour due to the load eccentricity and the tendency to uniform contact pressure. Since the static equilibrium must be ensured, the centroid of the contact pressure must be aligned with the eccentricity. Therefore, the effective length of contact is reduced to enable the development of a uniform contact pressure whose centroid is aligned with the point of load application. This reduction is schematically illustrated in Fig. 4.

Fig. 8 also presents the impact of the normalised size of the growth region ($D_g$) when a fixed eccentric load is applied. The initial pressure concentration tends to disappear when $D_g = 0.025$, while it remains present—yet, at a lower magnitude—when $D_g \ge 0.25$. We attribute this difference to the activation of growth at zones far from the contact region when $D_g \ge 0.25$. This activation did not have a negative impact in the previous section because the stress distribution was symmetric,
and hence, growth was symmetrically balanced. In contrast, in the present case, such a symmetry is lost due to the eccentricity of the load, as can be seen in the bottom row of Fig. 9. It presents the distribution of the growth stress at times $t = 0$ and $t = 200$ for $D_g = 0.025$ and $D_g \to \infty$. As can be seen, although the growth stress becomes symmetric around the contact zone, it remains asymmetric far from it when $D_g \to \infty$. In consequence, it generates an unbalanced growth. On the other hand, with $D_g = 0.025$, this issue is diminished, if not removed. Lastly, Fig. 9 also shows that the stress is better distributed as the contact pressure concentrations disappear.

Fig. 10 presents the response of the generative algorithm when the load eccentricity oscillates between $M_f = -0.5$ and $M_f = 0.5$. It is presented in terms of the evolution of the fields of maximum and mean contact pressure, see Eq. (18). As shown in Fig. 10a, the bio-inspired algorithm makes the maximum pressure quality increase; yet, it makes the mean pressure quality decrease.

The quality of the maximum pressure goes from 0.64 to 0.87, while the quality of the mean pressure goes from 0.84 to 0.64. In contrast, for the reference flat shape, the quality of the mean and maximum pressures are respectively 0.77 and 0.58. Consequently, the algorithm provides a better contact pressure distribution at every developmental state. Regarding the mean pressure, it is only better for the initial cycles of development, before it goes below 0.77. Yet, as can be seen in Fig. 10d both the maximum and mean pressure are highly concentrated at the edges for the flat shape. As can be seen, even though $Q_p$ of the mean pressure is higher for the developed shape at the last cycle, it is more concentrated in the flat shape.

The evolution of the maximum and mean pressure can also be seen in Figs. 10b and 10c. The highly concentrated pressure at the edges of the initial cycle disappears in both fields. In consequence, the contact pressure at the interior of the contact zone increases. The increase in the mean pressure at the centre seen in Fig. 10b translates into an augmentation of the rate of wear at this location; yet, it is diminished at the edges. Moreover, the decrease in the maximum contact pressure seen in Fig. 10c translates into a drop of the maximum shear stress by cycle. In turn, the fatigue damage slows down.

This can be more clearly seen in Fig. 11 which shows the evolution of the contact pressure distribution by cycle during the process of generation. Fig. 11a shows the maximum mean pressure $\max P_c^{\text{mean}}$ and the maximum pressure amplitude $\max P_c^{\text{ampl}}$ by cycle, where the initial cycle is set
as a reference, the star point. As the algorithm goes along, both maximum values decrease until the maximum mean pressure reaches its minimum (triangle symbol). Then, the maximum mean pressure starts to increase while the amplitude continues to decrease. At the circle point, the maximum mean pressure is equivalent to that of the initial cycle, but the maximum amplitude is lower. Later, the maximum amplitude reaches its minimum limiting value (square symbol). Thereafter, the maximum mean pressure keeps growing and the maximum amplitude remains constant. Figs. 11b and 11c show the pressure variation along the contact zone for the circle and square cycles compared to the initial cycle. As can be seen in Figs. 11a and 11c, the pressure variation at the square cycle has been reduced; the maximum pressure amplitude was reduced around 73% while the maximum mean pressure only augmented around 0.08%. This translates into a less demanding fatigue-loading environment.

### 4.3. From 2D contact interfaces to 3D cylindrical joints

In real engineering problems, the assumption of two-dimensional behaviour is not always suitable. In particular, cylindrical joints—whose surfaces are conformal—require the analysis in three dimensions. In this section, we test the performance of the shapes generated in the previous sections in a three-dimensional framework. We study the pressure distribution considering centred and eccentric radial loads. In addition, we also test the performance under axial forces.

For the construction of the cylindrical joint, let us start by considering a two-dimensional contact interface with contours $\Gamma_1^{c_e}$ (described by $f_1$) and $\Gamma_2^{c_e}$ (described by $f_2$) and of contact length $l_{cf}$, as those generated in the previous sections. The cylindrical joint is defined by the surfaces of revolution generated with these contours. The revolution is done with a mean diameter of $w$. A schematic representation of this three-dimensional framework is shown in Fig. 12, which also illustrates the loading conditions where the load $F_r$ is given by:

$$F_r = \frac{S_f E w l_{cf}}{2000}, \qquad (28)$$

provided that $w = l_{cf}$.

In terms of the size of the contact elements, the third row of Table 3 (error 3) shows the convergence with respect to the contact pressure at the reference contact line in the plane where
loads are applied (orange line in Fig. 12). As can be seen, with $h_c^{\lim}/w = 0.02$, the numerical error
is less than 4%. Since the domain does not change in this case, we shall use this mesh size.

As an example, Fig. 13a shows the distribution of the contact pressure using the category-II
geometry generated with a normalised initial size $B_e = 0.4$ and a stiffness factor $S_f = 0.1$ for
the application of radial and axial forces. As can be seen, when $F_a = 0$, the pressure is well
distributed in the contact region. This can also be seen in Fig. 13b which displays the pressure
on the reference contact line where loads are applied (orange line in Fig. 12) as $F_r$ is applied.
In contrast, as $F_a$ grows, the pressure accumulates at one side of the contact zone and increases
significantly, see Fig. 13c. Although the contact pressure is not uniform after the application of
the axial load, this figure demonstrates the potential of the generated geometries to support both
radial and axial forces.

Fig. 14a displays the contact pressure using the category-II geometry generated with $B_e = 0.4$ and $S_f = 0.1$ and a cyclic eccentric load with $M_f$ between -0.5 and 0.5. As can be seen, the initial contact zone (when $M_f = 0$) is smaller than that of the previous case. This is due to the additional adaptation to load eccentricity. The figure also shows that the contact zone becomes asymmetric when the eccentricity is applied, $M_f = 0.5$. Further, the red zone becomes larger indicating an augmentation of the contact pressure. Nevertheless, as can be seen in Fig. 14b, such an augmentation is less than 30% when the eccentricity is maximum. This figure also shows that in the three-dimensional framework the behaviour or the contact interface is similar to what was observed in the two-dimensional case: a field of mean pressure concentrated at the centre (compared with Fig. 10b) and a field of maximum pressure without sharp edges (compared with Fig. 10c).

## 5. Discussion

The results presented in this work reveal the robustness of the bio-inspired generative algorithm
for designing contact pairs with almost uniform contact pressure under time-independent loads.
The same set of growth parameters was able to provide such a contact pressure with different sizes,
stress states, and loading conditions. The results demonstrate that the obtained contact pressure
quality is independent of the stress states. Moreover, compared to a circular and a flat shape, the
contact pressure distribution is up to 34.7% more uniform. In addition, the results also show the versatility of the algorithm to deal with different loading conditions. Therefore, the bio-inspired growth rules here presented represent an implicit direction of growth, of geometric change, that searches a uniform distribution of the contact pressure and the augmentation of the contact region. This direction of growth can also be seen as an easier-to-compute alternative to the gradient in gradient-based optimisation. In addition, it is also different from heuristic optimisation methods in the sense that such methods explore several directions of change and the presented one has a unique direction.

The results also demonstrate the ability of the algorithm to deal with time-dependent eccentric loads. Since the time-dependent load is eccentric, the centroid of the resultant contact pressure shifts accordingly. Consequently, the pressure distribution changes at each instant. We identified two possible criteria for evaluating the pressure distribution over a cycle of load: the mean pressure and the maximum pressure at each point throughout the cycle. We found that the model privileges a uniform maximum over a uniform mean contact pressure, as seen in Fig. 10. The central zone of the contact interface is more likely to support load consistently throughout the cycle, whereas the edges are less active. For instance, when the eccentricity is on the right, the left edge supports less load, while the middle zone continues to support it. A uniform mean pressure would then suggest that pressure concentrates at the edges during moments of higher eccentricity, compensating for instances of lower pressure. In contrast, a uniform maximum pressure implies a more consistent pressure distribution at each instant, although a distribution shifts over time. In general, the change of the mean pressure leads to a reduction of the wear damage at the edges, but an acceleration of it at the middle. Regarding the change of the maximum pressure, it leads to a reduction on the fatigue damage. Thus, it is possible to say that the bio-inspired growth privileges the resistance to fatigue over the resistance to wear.

Additionally, the results also expose the possibility to utilise the generated profiles in a three-dimensional framework. Figs. 13 and 14 show that the behaviour of the contact interfaces seen in the two-dimensional framework remains. Therefore, the proposed bio-inspired algorithm can be used as an *in silico* tool to generate geometries adapted to their loading environment. A logical next step would be to explore suitable manufacturing techniques, such as additive manufacturing.
The results also show that the generated geometries may support both radial and axial loads. Yet, the distribution of the contact pressure was severely damaged when axial loads were considered. Thus, future work should study the adaptability to more complex loading conditions and the generation directly in a three-dimensional framework.

The presented algorithm also differs from other approaches, such as crowning or optimisation strategies, in that the contact interface is not the unique solution to a specific problem. Instead, the algorithm provides a family of contact interfaces whose contact pressure is almost uniform. In fact, the problem to be solved, Eqs. (1) to (4), is a rate of growth, but the domain change occurs in such a way that the contact pressure becomes and remains almost uniform. Moreover, the types of shapes seen in Fig. 6 also illustrate the variability of the response. This shape diversity is supplied by the fact that both contact shapes are subjected to growth and that the growth function directs the development towards a contact interface with almost uniform contact pressure. In [63], only one of the shapes was allowed to grow; hence, the generated profile was constrained by the fixed one. This is also the case of crowning and optimisation methods where one of the geometries is always considered fixed [13–15, 18, 24, 26, 32, 40]. In [62], both shapes were allowed to grow, but the problem was to obtain one optimal shape.

Although it is arguable to consider that the method should provide a unique shape, the possibility of obtaining a diversity of shapes allows the exploration of novel interfaces. Moreover, such diversity is largely seen in biological joints. For instance, interphalangeal joints are composed of a convex and a concave surface (shapes of category I and II) while elbow joints are composed of surfaces that have both convex and concave regions (shapes of category III). It may then be interesting to explore the advantages of each type of contact interface in future work.

## 6. Conclusion

In this work, we have considered a bio-inspired generative algorithm for the design of contact interfaces in engineering applications. It is inspired by the formation of endochondral bones and synovial joints. The algorithm is based on the modelling of evolving domains in contact with each other. Its principal characteristic is that the growth term is derived from the current biological theory of shape formation. The results presented here demonstrate that the growth function directs
the changes towards contact interfaces with almost uniform pressure in time-independent cases.
701
Therefore, it is a possible alternative to the gradient direction used in gradient-based optimisation
702
and to heuristic directions of change.
703

Moreover, the generative algorithm presented differentiates from other methods, such as crown-
ing and optimisation, in that the problem to be solved does not search for a specific contact in-
terface where pressure is well distributed. Instead, the development of well-distributed pressure
704
705
706
is an intrinsic part of the growth function. Therefore, the algorithm generates a family of shapes
707
adapted to the loading conditions.
708

In time-dependent cases, the presented generative design method leads to contact interfaces
709
with reduced variability of the contact pressure. In consequence, the generated shapes are adapted
710
to the loading environment as occurs in biological joints. This adaptation leads to a better distri-
bution of wear and fatigue damage.
711
712

Nonetheless, there are still limitations to be overcome and questions to be answered. Firstly,
713
the generative algorithm was developed in a two-dimensional framework where geometries could
714
grow freely. Although an equivalent geometry was proposed for three-dimensional cylindrical
715
joints, it is limited to cases where the diameter is equal to the length. Hence, an extension of it
716
into three-dimensional frameworks or enclosed spaces should be dealt with in future works. In
717
addition, this would allow us to consider higher eccentricity values since the static equilibrium
718
would be given by the outer joint. Secondly, since the model can provide different types of contact
719
interfaces, it is still an open question what advantages each of them has. Especially, given the
720
diversity of synovial joint shapes seen in nature. It is still not clear what impact could have the
721
inclusion of biochemical factors and their dynamics in the generated shapes as well as in the
722
process of generation itself. Therefore, an exploration of the role of biochemical factors could
723
be envisaged. For instance, the mechanical model can be coupled to a reaction-diffusion system
724
describing proliferative chondrocytes whose concentration would vary with the domain size as
725
well as with the loading environment. This could provide a mechanism of growth regulation.
726
Additionally, the current model considers that the mechanical stimuli have an immediate impact on
727
the growth force; yet, it is more plausible that growth is a delayed response to mechanical stimuli
728
from previous time instants. By incorporating a reaction-diffusion system, this delay between
729
730 stimulus and effect could be accurately modelled. Further, the intrinsic change of mechanical
731 properties seen in formation during ossification has not been taken into account either. Finally,
732 the results presented here correspond to numerical simulations, and future work can also focus on
733 manufacturability.

734 Declaration of competing interest

735 The authors declare that they have no known competing financial interests or personal rela-
736 tionships that could have appeared to influence the work reported in this paper.

737 Acknowledgements

738 This research was supported by the French Research National Agency (ANR) Grant No. ANR-
739 20-CE10-0008.

740 References

[1] J. A. Tichy, D. M. Meyer, Review of solid mechanics in tribology, International Journal of Solids and Structures 37 (2000) 391–400. doi:10.1016/s0020-7683(99)00101-8.

[2] R. Hjelm, L. Everlid, E. Bergseth, F. Reinle, B. Brodmann, M. Tu, L. Bard, J. Wahlström, A multi-perspective method for gear efficiency and contact analysis, Results in Engineering 20 (2023) 101582. doi:10.1016/j.rineng.2023.101582.

[3] P. Põdra, S. Andersson, Simulating sliding wear with finite element method, Tribology International 32 (1999) 71–81. doi:10.1016/s0301-679x(99)00012-2.

[4] J. Zhang, F. Lyu, B. Xu, W. Huang, W. Wu, Z. Guo, H. Xu, X. Huang, Simulation and experimental investigation on low wear rate surface contour of piston/cylinder pair in an axial piston pump, Tribology International 162 (2021) 107127. doi:10.1016/j.triboint.2021.107127.

[5] D. Teng, J. Wang, C. Li, X. Sa, Investigation of friction and wear behavior of cast aluminum alloy piston skirt with graphite coating using a designed piston skirt test apparatus, Materials 15 (2022) 4010. doi:10.3390/ma15114010.

[6] P. Bergmann, F. Grün, F. Summer, I. Gódor, Evaluation of wear phenomena of journal bearings by close to component testing and application of a numerical wear assessment, Lubricants 6 (2018) 65. doi:10.3390/lubricants6030065.

[7] A. A. Schmidt, T. Schmidt, O. Grabherr, D. Bartel, Transient wear simulation based on three-dimensional finite element analysis for a dry running tilted shaft-bushing bearing, Wear 408–409 (2018) 171–179. doi:10.1016/j.wear.2018.05.008.

[8] K. Nakazawa, N. Maruyama, T. Hanawa, Effect of contact pressure on fretting fatigue of austenitic stainless steel, Tribology International 36 (2003) 79–85. doi:10.1016/s0301-679x(02)00135-4.

[9] K. A. Dutt, S. J. Joshi, D. B. Shah, S. B. Soni, D. Prajapati, Experimental investigation of fatigue for one way power transmission clutches using accelerated life tests, Results in Engineering 22 (2024) 102136. doi:10.1016/j.rineng.2024.102136.

[10] R. G. Budynas, J. K. Nisbett, Fatigue Failure Resulting from Variable Loading, ninth edition ed., Mc Graw Hill, 2011, pp. 265–356.

[11] R. L. Norton, Fatigue failure theory, fourth edition ed., Pearson, 2011, pp. 303–418.
768

769

770

771

772

773

774

775

776

777

778

779

780

781

782

783

784

785

786

787

788

789

790

791

792

793

794

795

796

797

798

799

800

801

802

803

804

805

806

807

808

809

810

811

812

813

814

815

816

817

818

[12] G. Lundberg, Elastische berührung zweier halbräume, Forschung auf dem Gebiete des Ingenieurwesens 10 (1939) 201–211. doi:10.1007/bf02584950.

[13] P. Johns, R. Gohar, Roller bearings under radial and eccentric loads, Tribology International 14 (1981) 131–136. doi:10.1016/0301-679x(81)90058-x.

[14] K. S. Kumar, R. Tiwari, P. V. V. N. Prasad, An optimum design of crowned cylindrical roller bearings using genetic algorithms, Journal of Mechanical Design 131 (2009). doi:10.1115/1.3116344.

[15] L. Cui, Y. He, A new logarithmic profile model and optimization design of cylindrical roller bearing, Industrial Lubrication and Tribology 67 (2015) 498–508. doi:10.1108/ilt-01-2015-0006.

[16] Z. Wu, Y. Xu, S. Deng, K. Liu, Study on logarithmic crowning of cylindrical roller profile considering angular misalignment, Journal of Mechanical Science and Technology 34 (2020) 2111–2120. doi:10.1007/s12206-020-0130-7.

[17] Z. Wu, Y. Xu, K. Liu, Study on logarithmic crowning of tapered roller profile considering angular misalignment, Journal of Tribology 142 (2020). doi:10.1115/1.4047106.

[18] D. Hilding, A. Klarbring, J. Petersson, Optimization of structures in unilateral contact, Applied Mechanics Reviews 52 (1999) 139–160. doi:10.1115/1.3098931.

[19] J. Cho, S. Shin, W. Yoo, Crown shape optimization for enhancing tire wear performance by ann, Computers & Structures 83 (2005) 920–933. doi:10.1016/j.compstruc.2004.11.011.

[20] W. Li, Q. Li, G. P. Steven, Y. M. Xie, An evolutionary shape optimization procedure for contact problems in mechanical designs, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 217 (2003) 435–446. doi:10.1243/095440603321509711.

[21] M. Najjari, R. Guilbault, Formula derived from particle swarm optimization (pso) for optimum design of cylindrical roller profile under ehl regime, Mechanism and Machine Theory 90 (2015) 162–174. doi:10.1016/j.mechmachtheory.2015.03.016.

[22] M. Duchemin, C. Tugui, V. Collee, Optimization of contact profiles using super-ellipse, SAE International Journal of Materials and Manufacturing 10 (2017) 234–244. doi:10.4271/2017-01-1349.

[23] B. D. Upadhyay, S. S. Sonigra, S. D. Daxini, Numerical analysis perspective in structural shape optimization: A review post 2000, Advances in Engineering Software 155 (2021) 102992. doi:10.1016/j.advengsoft.2021.102992.

[24] A. Klarbring, J. Haslinger, On almost constant contact stress distributions by shape optimization, Structural Optimization 5 (1993) 213–216. doi:10.1007/bf01743581.

[25] E. A. Fancello, J. Haslinger, R. A. Feijóo, Numerical comparison between two cost functions in contact shape optimization, Structural Optimization 9 (1995) 57–68. doi:10.1007/bf01742646.

[26] E. Bretin, J. Chapelat, P.-Y. Outtier, Y. Renard, Shape optimization of a linearly elastic rolling structure under unilateral contact using nitsche's method and cut finite elements, Computational Mechanics 70 (2022) 205–224. doi:10.1007/s00466-022-02164-z.

[27] E. Bretin, J. Chapelat, C. Douanla-Lontsi, T. Homolle, Y. Renard, Shape sensitivity analysis of an elastic contact problem: Convergence of the nitsche based finite element approximation, Nonlinear Analysis: Real World Applications 72 (2023) 103836. doi:10.1016/j.nonrwa.2023.103836.

[28] A. Klarbring, M. Rönnqvist, Nested approach to structural optimization in nonsmooth mechanics, Structural Optimization 10 (1995) 79–86. doi:10.1007/bf01743534.

[29] I. Páczelt, T. Szabó, Optimal shape design for contact problems, Structural Optimization 7 (1994) 66–75. doi:10.1007/bf01742507.

[30] I. Páczelt, Z. Mróz, On optimal contact shapes generated by wear, International Journal for Numerical Methods in Engineering 63 (2005) 1250–1287. doi:10.1002/nme.1321.

[31] I. Páczelt, Z. Mróz, Optimal shapes of contact interfaces due to sliding wear in the steady relative motion, International Journal of Solids and Structures 44 (2007) 895–925. doi:10.1016/j.ijsolstr.2006.05.027.

[32] D. Hilding, B. Torstenfelt, A. Klarbring, A computational methodology for shape optimization of structures in frictionless contact, Computer Methods in Applied Mechanics and Engineering 190 (2001) 4043–4060. doi:10.1016/s0045-7825(00)00310-8.

[33] Y. Wei, Y. Qin, R. Balendra, Q. Jiang, Fe analysis of a novel roller form: a deep end-cavity roller for roller-type bearings, Journal of Materials Processing Technology 145 (2004) 233–241. doi:10.1016/s0924-0136(03)
[34] A. Myśliński, Level set method for optimization of contact problems, Engineering Analysis with Boundary Elements 32 (2008) 986-994. doi:10.1016/j.enganabound.2007.12.008.

[35] A. Myśliński, Piecewise constant level set method for topology optimization of unilateral contact problems, Advances in Engineering Software 80 (2015) 25-32. doi:10.1016/j.advengsoft.2014.09.020.

[36] A. Myśliński, M. Wróblewski, Structural optimization of contact problems using cahn-hilliard model, Computers & Structures 180 (2017) 52-59. doi:10.1016/j.compstruc.2016.03.013.

[37] C. Niu, W. Zhang, T. Gao, Topology optimization of continuum structures for the uniformity of contact pressures, Structural and Multidisciplinary Optimization 60 (2019) 185-210. doi:10.1007/s00158-019-02208-8.

[38] C. Niu, W. Zhang, T. Gao, Topology optimization of elastic contact problems with friction using efficient adjoint sensitivity analysis with load increment reduction, Computers & Structures 238 (2020) 106296. doi:10.1016/j.compstruc.2020.106296.

[39] F. Fernandez, M. A. Puso, J. Solberg, D. A. Tortorelli, Topology optimization of multiple deformable bodies in contact with large deformations, Computer Methods in Applied Mechanics and Engineering 371 (2020) 113288. doi:10.1016/j.cma.2020.113288.

[40] H. Kristiansen, K. Poulios, N. Aage, Topology optimization for compliance and contact pressure distribution in structural problems with friction, Computer Methods in Applied Mechanics and Engineering 364 (2020) 112915. doi:10.1016/j.cma.2020.112915.

[41] Q. Lin, N. Yang, J. Hong, L. Liu, Y. Zhang, An optimization design of contact interface material stiffness for improving the uniformity of the contact pressure, Journal of Tribology 142 (2020). doi:10.1115/1.4046172.

[42] Q. Lin, Y. Zhou, N. Yang, J. Hong, C. Wang, Design of component structure in assemblies for simultaneously regulating contact pressure distribution and natural frequencies, European Journal of Mechanics - A/Solids 94 (2022) 104557. doi:10.1016/j.euromechsol.2022.104557.

[43] Y. Zhou, Q. Lin, J. Hong, N. Yang, Bidirectional evolutionary optimization design of material stiffness for the uniformity of the contact stress, European Journal of Mechanics - A/Solids 89 (2021) 104288. doi:10.1016/j.euromechsol.2021.104288.

[44] Y. Zhou, Q. Lin, X. Yang, J. Hong, N. Zhang, F. Zhao, Material stiffness optimization for contact stress distribution in frictional elastic contact problems with multiple load cases, International Journal of Mechanics and Materials in Design 17 (2021) 503-519. doi:10.1007/s10999-021-09544-y.

[45] Y. Zhou, Q. Lin, W. Li, T. Yue, J. Yan, J. Hong, Optimization of elastoplastic behavior of contact interface for improved contact stress distribution, Mechanics of Advanced Materials and Structures 30 (2022) 803-813. doi:10.1080/15376494.2021.2024926.

[46] Y. Zhou, Q. Lin, C. Wang, J. Guo, J. Yan, J. Hong, Material stiffness optimization for homogenizing contact stress distribution based on particle swarm optimization using elite opposition-based learning mutation, Mechanics of Advanced Materials and Structures (2023) 1-13. doi:10.1080/15376494.2023.2283809.

[47] H. Strand, Journal bearing housing design—a statistical study with fem, Tribology International 40 (2007) 665-671. doi:10.1016/j.triboint.2005.11.018.

[48] B. K. Prasad, Sliding wear response of a bronze bushing: Influence of applied load and test environment, Journal of Materials Engineering and Performance 21 (2012) 2155-2164. doi:10.1007/s11665-012-0139-x.

[49] H. Jia, J. Li, J. Wang, K. Xiao, Y. Han, A lubricated wear model for determining wear surface geometry on journal-bearing surfaces, Surface Topography: Metrology and Properties 8 (2020) 035002. doi:10.1088/2051-672x/ab9fd2.

[50] W. Zhou, Y. Zhao, H. Yuan, J. Liu, X. Wang, Research on the contact pressure calculation method for the misaligned elastomeric journal bearing, Journal of Marine Science and Engineering 10 (2022) 141. doi:10.3390/jmse10020141.

[51] J. Wahlström, A pin-on-disc tribometer study of friction at low contact pressures and sliding speeds for a disc brake material combination, Results in Engineering 4 (2019) 100051. doi:10.1016/j.rineng.2019.100051.

[52] D. Sysaykeo, E. Mermoz, T. Thouveny, Clearance and design optimization of bio-inspired bearings under off-center load, CIRP Annals 69 (2020) 121-124. doi:10.1016/j.cirp.2020.03.006.

[53] D. Sysaykeo, J.-M. Linares, E. Mermoz, Wear behavior of a bio-inspired bearing for off-center loads, Journal of Bionic Engineering 17 (2020) 1251-1262. doi:10.1007/s42235-020-0107-3.
[54] A. Sanz-Idirin, S. Arroyave-Tobon, J.-M. Linares, P. J. Arrazola, Load bearing performance of mechanical joints inspired by elbow of quadrupedal mammals, Bioinspiration & Biomimetics 16 (2021) 046025. doi:10.1088/1748-3190/abeb57.

[55] C. Mattheck, S. Burkhardt, A new method of structural shape optimization based on biological growth, International Journal of Fatigue 12 (1990) 185-190. doi:10.1016/0142-1123(90)90094-u.

[56] C. Mattheck, D. Erb, K. Bethge, U. Begemann, Three-dimensional shape optimization of a bar with a rectangular hole, Fatigue & Fracture of Engineering Materials & Structures 15 (1992) 347-351. doi:10.1111/j.1460-2695.1992.tb01276.x.

[57] J. Chen, W. Tsai, Shape optimization by using simulated biological growth approaches, in: 4th Symposium on Multidisciplinary Analysis and Optimization, American Institute of Aeronautics and Astronautics, 1992, pp. 895-901. doi:10.2514/6.1992-4802.

[58] M. Bagge, A model of bone adaptation as an optimization process, Journal of Biomechanics 33 (2000) 1349-1357. doi:10.1016/s0021-9290(00)00124-x.

[59] A. Andrade-Campos, A. Ramos, J. A. Simões, A model of bone adaptation as a topology optimization process with contact, Journal of Biomedical Science and Engineering 05 (2012) 229-244. doi:10.4236/jbise.2012.55030.

[60] D. R. Carter, G. S. Beaupré, M. Wong, R. L. Smith, T. P. Andriacchi, D. J. Schurman, The mechanobiology of articular cartilage development and degeneration, Clinical Orthopaedics & Related Research 427 (2004) S69-S77. doi:10.1097/01.blo.0000144970.05107.7e.

[61] S. Arroyave-Tobon, K. Marquez-Florez, P. Heymann, J.-M. Linares, Generative design of joint contact surfaces inspired by biological morphogenesis, CIRP Annals 71 (2022) 125-128. doi:10.1016/j.cirp.2022.04.006.

[62] K. Marquez-Florez, S. Arroyave-Tobon, J.-M. Linares, From biological morphogenesis to engineering joint design: A bio-inspired algorithm, Materials & Design 225 (2023) 111466. doi:10.1016/j.matdes.2022.111466.

[63] D. Hernandez-Aristizabal, S. Arroyave-Tobon, D. A. Garzon-Alvarado, J.-M. Linares, Stress-adaptive design of 2d contact interfaces with uniform pressure: A bio-inspired approach, International Journal of Solids and Structures 270 (2023) 112238. doi:10.1016/j.ijsolstr.2023.112238.

[64] J. H. Heegaard, G. S. Beaupré, D. R. Carter, Mechanically modulated cartilage growth may regulate joint surface morphogenesis, Journal of Orthopaedic Research 17 (1999) 509-517. doi:10.1002/jor.1100170408.

[65] M. Giorgi, A. Carriero, S. J. Shefelbine, N. C. Nowlan, Mechanobiological simulations of prenatal joint morphogenesis, Journal of Biomechanics 47 (2014) 989-995. doi:10.1016/j.jbiomech.2014.01.002.

[66] M. Giorgi, A. Carriero, S. J. Shefelbine, N. C. Nowlan, Effects of normal and abnormal loading conditions on morphogenesis of the prenatal hip joint: application to hip dysplasia, Journal of Biomechanics 48 (2015) 3390-3397. doi:10.1016/j.jbiomech.2015.06.002.

[67] A. F. Carrera-Pinzón, K. Márquez-Flórez, R. H. Kraft, S. Ramtani, D. A. Garzón-Alvarado, Computational model of a synovial joint morphogenesis, Biomechanics and Modeling in Mechanobiology 19 (2019) 1389-1402. doi:10.1007/s10237-019-01277-4.

[68] S. M. Sadeghian, F. D. Shapiro, S. J. Shefelbine, Computational model of endochondral ossification: Simulating growth of a long bone, Bone 153 (2021) 116132. doi:10.1016/j.bone.2021.116132.

[69] J. Kahn, Y. Shwartz, E. Blitz, S. Krief, A. Sharir, D. A. Breitel, R. Rattenbach, F. Relaix, P. Maire, R. B. Rountree, D. M. Kingsley, E. Zelzer, Muscle contraction is necessary to maintain joint progenitor cell fate, Developmental Cell 16 (2009) 734-743. doi:10.1016/j.devcel.2009.04.013.

[70] Y. Shwartz, E. Blitz, E. Zelzer, One load to rule them all: Mechanical control of the musculoskeletal system in development and aging, Differentiation 86 (2013) 104-111. doi:10.1016/j.diff.2013.07.003.

[71] A. D. Berendsen, B. R. Olsen, Bone development, Bone 80 (2015) 14-18. doi:10.1016/j.bone.2015.04.035.

[72] N. Felsenthal, E. Zelzer, Mechanical regulation of musculoskeletal system development, Development 144 (2017) 4271-4283. doi:10.1242/dev.151266.

[73] E. Kozhemyakina, A. B. Lassar, E. Zelzer, A pathway to bone: signaling molecules and transcription factors involved in chondrocyte development and maturation, Development 142 (2015) 817-831. doi:10.1242/dev.105536.
[74] R. S. Decker, Articular cartilage and joint development from embryogenesis to adulthood, Seminars in Cell & Developmental Biology 62 (2017) 50–56. doi:10.1016/j.semcdb.2016.10.005.

[75] D. R. Carter, M. Wong, The role of mechanical loading histories in the development of diarthrodial joints, Journal of Orthopaedic Research 6 (1988) 804–816. doi:10.1002/jor.1100060604.

[76] V. S. Salazar, L. W. Gamer, V. Rosen, Bmp signalling in skeletal development, disease and repair, Nature Reviews Endocrinology 12 (2016) 203–221. doi:10.1038/nrendo.2016.12.

[77] G. L. Galea, M. R. Zein, S. Allen, P. Francis-west, Making and shaping endochondral and intramembranous bones, Developmental Dynamics 250 (2020) 414–449. doi:10.1002/dvdy.278.

[78] D. R. Carter, T. E. Orr, D. P. Fyhrie, D. J. Schurman, Influences of mechanical stress on prenatal and post-natal skeletal development, Clinical Orthopaedics and Related Research 219 (1987) 237–250. doi:10.1097/00003086-198706000-00034.

[79] M. P. Dieterle, A. Husari, B. Rolauffs, T. Steinberg, P. Tomakidi, Integrins, cadherins and channels in cartilage mechanotransduction: perspectives for future regeneration strategies, Expert Reviews in Molecular Medicine 23 (2021). doi:10.1017/erm.2021.16.

[80] T. Ogura, T. Minas, A. Tsuchiya, S. Mizuno, Effects of hydrostatic pressure and deviatoric stress on human articular chondrocytes for designing neo-cartilage construct, Journal of Tissue Engineering and Regenerative Medicine (2019). doi:10.1002/term.2863.

[81] E. F. Morgan, G. U. Unnikrisnan, A. I. Hussein, Bone mechanical properties in healthy and diseased states, Annual Review of Biomedical Engineering 20 (2018) 119–143. doi:10.1146/annurev-bioeng-062117-121139.

[82] A. K. Williamson, A. C. Chen, K. Masuda, E. J. M. A. Thonar, R. L. Sah, Tensile mechanical properties of bovine articular cartilage: Variations with growth and relationships to collagen network components, Journal of Orthopaedic Research 21 (2003) 872–880. doi:10.1016/s0736-0266(03)00030-5.

[83] K. Takahashi, T. Kubo, K. Kobayashi, J. Imanishi, M. Tkigawa, Y. Arai, Y. Hirasawa, Hydrostatic pressure influences mrna expression of transforming growth factor-β1 and heat shock protein 70 in chondrocyte-like cell line, Journal of Orthopaedic Research 15 (1997) 150–158. doi:10.1002/jor.1100150122.

[84] M. A. Lewis, J. D. Murray, Analysis of stable two-dimensional patterns in contractile cytosol, Journal of Nonlinear Science 1 (1991) 289–311. doi:10.1007/bf01238816.

[85] L. Murphy, A. Madzvamuse, A moving grid finite element method applied to a mechanobiochemical model for 3d cell migration, Applied Numerical Mathematics 158 (2020) 336–359. doi:10.1016/j.apnum.2020.08.004.

[86] A. Stephanou, A mathematical model for the dynamics of large membrane deformations of isolated fibroblasts, Bulletin of Mathematical Biology 66 (2004) 1119–1154. doi:10.1016/j.bulm.2003.11.004.

[87] A. Madzvamuse, A. J. Wathen, P. K. Maini, A moving grid finite element method applied to a model biological pattern generator, Journal of Computational Physics 190 (2003) 478–500. doi:10.1016/s0021-9991(03)00294-8.

[88] S. Chiravarambath, N. K. Simha, R. Namani, J. L. Lewis, Poroviscoelastic cartilage properties in the mouse from indentation, Journal of Biomechanical Engineering 131 (2008). doi:10.1115/1.3005199.

[89] C. Sandino, D. D. Mcerlain, J. Schipilow, S. K. Boyd, The poro-viscoelastic properties of trabecular bone: a micro computed tomography-based finite element study, Journal of the Mechanical Behavior of Biomedical Materials 44 (2015) 1–9. doi:10.1016/j.jmbbm.2014.12.018.

[90] K. Márquez-Flórez, D. Garzón-Alvarado, C. Carda, M. Sancho-Tello, Computational model of articular cartilage regeneration induced by scaffold implantation in vivo, Journal of Theoretical Biology 561 (2023) 111393. doi:10.1016/j.jtbi.2022.111393.

[91] Électricité de France. Finite element code_aster, Analysis of Structures and Thermomechanics for Studies and Research, 1989–2024. URL: www.code-aster.org, published: open source.

[92] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering 79 (2009) 1309–1331. doi:10.1002/nme.2579.

[93] R. Durand, B. Pantoja-Rosero, V. Oliveira, A general mesh smoothing method for finite elements, Finite Elements in Analysis and Design 158 (2019) 17–30. doi:10.1016/j.finel.2019.01.010.
972

973

[94] R. Myers, D. Montgomery, C. Anderson-Cook, The Analysis of Second-Order Response Surfaces, Wiley, 2016, pp. 273–324.

974

975

[95] J. P. Kleijnen, Design and Analysis of Simulation Experiments, Springer International Publishing, 2015. doi:10.1007/978-3-319-18087-8.

976

977

978

[96] N. Castro-González, F. M. Dopico, J. M. Molera, Multiplicative perturbation theory of the moore–penrose inverse and the least squares problem, Linear Algebra and its Applications 503 (2016) 1–25. doi:10.1016/j.laa.2016.03.027.

979

980

[97] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, S. Tarantola, Global Sensitivity Analysis. The Primer, Wiley, 2007. doi:10.1002/9780470725184.

981

982

983

[98] F. König, A. Ouald Chaib, G. Jacobs, C. Sous, A multiscale-approach for wear prediction in journal bearing systems – from wearing-in towards steady-state wear, Wear 426–427 (2019) 1203–1211. doi:10.1016/j.wear.2019.01.036.

984

985

[99] K. Marquez-Florez, S. Arroyave-Tobon, L. Tadrist, J.-M. Linares, Elbow dimensions in quadrupedal mammals driven by lubrication regime, Scientific Reports 14 (2024). doi:10.1038/s41598-023-50619-x.
Nomenclature

<table><tr><td>α<sub>g</sub></td><td>Growth rate parame- ter.</td><td>D<sub>g</sub></td><td>Normalised size of the growth region.</td><td>P<sub>c</sub><sup>ampl</sup></td><td>Field of the amplitude of the contact pressure during a cycle.</td></tr><tr><td>β</td><td>Regression parame- ters.</td><td>d<sub>Γ<sub>i</sub><sup>c</sup>(t)</sub></td><td>Minimum distance from a point to Γ<sub>i</sub><sup>c</sup>(t).</td><td>P<sub>c</sub><sup>max</sup></td><td>Field of maximum contact pressure dur- ing a cycle.</td></tr><tr><td>Γ<sup>c</sup></td><td>Contact boundary.</td><td>d<sub>lim</sub></td><td>Reference size of the growth zone.</td><td>P<sub>c</sub><sup>mean</sup></td><td>Field of mean contact pressure during a cy- cle.</td></tr><tr><td>δ</td><td>Initial gap between contacting surfaces.</td><td>E</td><td>Young modulus.</td><td>P<sub>e</sub></td><td>Element perimeter.</td></tr><tr><td>ε</td><td>Linear strain tensor.</td><td>e</td><td>External force eccen- tricity.</td><td>p<sub>eq</sub></td><td>Uniform equivalent contact pressure.</td></tr><tr><td>ϑ</td><td>Sigmoid function smoothness.</td><td>e<sup>lim</sup><sub>q</sub></td><td>Minimum element quality.</td><td>Q<sub>κ</sub></td><td>Normalised curvature norm.</td></tr><tr><td>κ</td><td>Mean curvature value.</td><td>F</td><td>External force.</td><td>Q<sub>c</sub></td><td>Conformity between two contacting sur- faces.</td></tr><tr><td>ν</td><td>Poisson's coefficient.</td><td>f<sub>b1</sub></td><td>Cell disposition to grow due to mechani- cal stimuli.</td><td>Q<sub>p</sub></td><td>Contact pressure qual- ity.</td></tr><tr><td>σ</td><td>Stress tensor.</td><td>f<sub>b2</sub></td><td>Basal growth.</td><td>r<sub>e</sub></td><td>Metamodel residual vector.</td></tr><tr><td>σ<sub>g</sub></td><td>Growth tensor</td><td>f<sub>m</sub></td><td>Mechanical contribu- tion to growth.</td><td>RMSE</td><td>Root mean square er- ror.</td></tr><tr><td>σ<sup>*</sup><sub>hyc</sub></td><td>Normalised hydro- static stress.</td><td>g</td><td>Growth rule.</td><td>S</td><td>Sigmoid function.</td></tr><tr><td>σ<sup>*</sup><sub>lim</sub></td><td>Hydrostatic stress threshold.</td><td>h<sup>lim</sup><sub>c</sub></td><td>Maximum contact ele- ment size.</td><td>S<sub>f</sub></td><td>Stiffness factor.</td></tr><tr><td>σ<sup>*</sup><sub>v</sub></td><td>Normalised shear stress.</td><td>l<sub>c</sub></td><td>Contact length.</td><td>t</td><td>Time.</td></tr><tr><td>τ<sup>*</sup><sub>lim</sub></td><td>Shear stress threshold.</td><td>M<sub>f</sub></td><td>Eccentricity factor.</td><td>t</td><td>External traction.</td></tr><tr><td>a</td><td>Vector of biochemical species.</td><td>n</td><td>Outward-pointing unit normal vector.</td><td>u</td><td>Displacement field.</td></tr><tr><td>A<sub>e</sub></td><td>Element area.</td><td>p<sub>c</sub></td><td>Magnitude of the con- tact pressure.</td><td></td><td></td></tr></table>
<table><caption>Table 1: Computational models of bone-joint formation. SIG: stress-independent growth. HDG: hydrostatic-pressure dependent growth. LT: type of load transfer between growth joints. OI: osteogenic index. Oss: ossification. SC: shape congruence.</caption><thead><tr><th>Ref.</th><th>SIG</th><th>HDG</th><th>LT</th><th>OI</th><th>Oss</th><th>SC</th></tr></thead><tbody><tr><td>[64]</td><td>✓</td><td>✗</td><td>Direct contact</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>[64]</td><td>✓</td><td>✓</td><td>Direct contact</td><td>✗</td><td>✗</td><td>✓</td></tr><tr><td>[65]</td><td>✓</td><td>✓</td><td>Elastic capsule</td><td>✗</td><td>✗</td><td>✓</td></tr><tr><td>[66]</td><td>✓</td><td>✓</td><td>Direct contact</td><td>✗</td><td>✗</td><td>✓</td></tr><tr><td>[67]</td><td>✗</td><td>✓</td><td>Elastic capsule</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>[68]</td><td>✓</td><td>✓</td><td>Hertzian pressure</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

<table><caption>Table 2: Stiffness factors in engineering and biology. For engineering cases, the pressure at the contact interface can be calculated dividing the applied load by the projected area of contact. In other words, if a load *F* is applied on a bushing of diameter *D* and length *L*, the projected pressure is given by *F*/*(DL)*, and *S*<sub>*f*</sub> = 1000*F*/*(DLE)*.</caption><thead><tr><th>Ref.</th><th>S<sub>f</sub><sup>min</sup></th><th>S<sub>f</sub><sup>max</sup></th></tr></thead><tbody><tr><td>[48]</td><td>0.319</td><td>0.435</td></tr><tr><td>[98]</td><td>0.526</td><td>0.526</td></tr><tr><td>[49]</td><td>0.499</td><td>0.499</td></tr><tr><td>[7]</td><td>0.001</td><td>0.001</td></tr><tr><td>[47]</td><td>0.476</td><td>5.238</td></tr><tr><td>[50]</td><td>8.183</td><td>8.183</td></tr><tr><td>[81, 99]</td><td>0.8</td><td>240</td></tr><tr><td>This work</td><td>0.1</td><td>10</td></tr></tbody></table>

<table><caption>Table 3: Convergence of the mesh size of the two-dimensional framework. Error 1 is the error for the case of centred load, error 2 for eccentric load and error 3 for the three-dimensional framework. The error is computed as $||p_c^{\text{lim}/2} - p_c^{\text{lim}}||_{L_2(\Gamma_{\text{ce}})} / ||p_c^{\text{lim}/2}||_{L_2(\Gamma_{\text{ce}})}$. For the two-dimensional framework, the size used was 0.0025, with an error less than 2.5%, and for the three-dimensional framework, it was 0.02, with an error less than 4%. Notice that for error 2, a shape obtained in the case of centred load was used; hence 0.0025 is the minimum size allowed.</caption><thead><tr><th>h<sub>c</sub><sup>lim</sup>/w</th><th>Error 1</th><th>Error 2</th><th>Error 3</th></tr></thead><tbody><tr><td>0.08</td><td>0.284</td><td>0.057</td><td>0.107</td></tr><tr><td>0.04</td><td>0.157</td><td>0.050</td><td>0.075</td></tr><tr><td>0.02</td><td>0.083</td><td>0.015</td><td>0.040</td></tr><tr><td>0.01</td><td>0.065</td><td>0.012</td><td></td></tr><tr><td>0.005</td><td>0.035</td><td>0.009</td><td></td></tr><tr><td>0.0025</td><td>0.025</td><td></td><td></td></tr><tr><td>0.00125</td><td>0.010</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Figure 2: Algorithm to approximate the solution of Eqs. (1) to (4), using finite elements with code_aster [91] (blue thick contours) and the meshing tool Gmsh [92] (black thick contours).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Schematic representations of the geometry and boundary conditions: (a) represents the initial geometry (a semi-circle of radius *R* in contact with a flat block) and (b) the final geometry. *F* and *e* are respectively the applied force and the eccentricity. The upper right graph of each figure illustrates the distribution of the contact pressure. The lower right graph shows the contact profiles as *y* = *f*₁(*x*) and *y* = *f*₂(*x*), where *a* and *b* are the limits of the contact boundary in *x*. *l*<sub>c0</sub> and *l*<sub>cf</sub> indicate the initial and the final length of the contact zone, respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Free body diagram of a schematic representation of the upper geometry. The eccentricity or misalignment *e* must be lower than a threshold to ensure the existence of an equilibrium state.</figcaption></figure>

<table><caption>Table 4: Parameters of the metamodels for *Q*<sub>*p*</sub>, *Q*<sub>*c*</sub> and *Q*<sub>*κ*₁</sub>, see Eq. (19).</caption><thead><tr><th></th><th>Q<sub>p</sub></th><th>Q<sub>c</sub></th><th>Q<sub>κ₁</sub></th></tr></thead><tbody><tr><td>β<sub>0</sub></td><td>0.96 ± 0.01</td><td>0.99 ± 0.12</td><td>0.56 ± 0.11</td></tr><tr><td>β<sub>1</sub> → B<sub>e</sub></td><td>0</td><td>-0.13 ± 0.07</td><td>-0.30 ± 0.08</td></tr><tr><td>β<sub>2</sub> → S<sub>f</sub></td><td>0</td><td>0</td><td>0.35 ± 0.08</td></tr><tr><td>β<sub>12</sub> → B<sub>e</sub>S<sub>f</sub></td><td>0.01 ± 0.01</td><td>0.09 ± 0.08</td><td>-0.17 ± 0.10</td></tr><tr><td>β<sub>11</sub> → B<sub>e</sub><sup>2</sup></td><td>-0.02 ± 0.01</td><td>-0.11 ± 0.11</td><td>-0.16 ± 0.14</td></tr><tr><td>β<sub>22</sub> → S<sub>f</sub><sup>2</sup></td><td>0</td><td>-0.08 ± 0.11</td><td>0</td></tr><tr><td>r<sup>2</sup></td><td>0.72</td><td>0.88</td><td>0.98</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 5: Results of the metamodels. (a) indicates the refinement of the metamodels removing less significant variables, (b) illustrates the impact of each variable on the response, the up and down arrows indicate whether the relation is direct or inverse, respectively, and (c) shows $Q_p$, $Q_c$ and $Q_{k1}$ as functions of the stiffness factor ($S_f$) and of the normalised initial size ($B_e$), for the model with no eccentricity ($M_f = 0$) and a growth region covering the whole domain ($D_g \to \infty$). The ranges of $S_f$ and $B_e$ are logarithmic and normalised to the standard reference axis [-1, 1]. The sampling points are defined with $S_f \in \{0.1, 1, 10\}$ and $B_e \in \{0.2, 0.4, 0.8\}$. The scatter points represent the real values of the variables and the surfaces represent the fitted model (Eq. (19)). The reddish surfaces indicate the average behaviour and the gray ones the variation.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Generated shapes varying $B_e$ and $S_f$ and considering that $M_f = 0$ and $D_g \to \infty$. Each contact pair is classified into one of the following categories: (I) convex $f_1$ and concave $f_2$, (II) concave $f_1$ and convex $f_2$, and (III) $f_1$ and $f_2$ have convex and concave regions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Characteristic evolution of (a) the contact pressure distribution, (b) the contact length and (c) the contact pressure quality. This case corresponds to $B_e = 0.4$, $S_f = 0.1$ and $M_f = 0$ and $D_g \to \infty$. In (a) Each line represents the contact pressure at a given time: the yellow and lightest line is the initial pressure distribution while the red and darkest line is the final one. The abscissa contains the arc length of the contact boundary. The ordinate is the normalised contact pressure using as reference pressure the maximum of the initial contact pressure. In (b) and (c) the contact length is presented normalised with respect to the initial contact length.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Impact of the normalised size of the growth region ($D_g$) on the evolution of the contact pressure distribution under the application of a fixed eccentric load corresponding to an eccentricity factor $M_f = 0.5$. The pressure is normalised using its maximum value at $t = 0$. Categories I, II and III refer to the shape of the initial geometry taking as examples the shapes presented in Fig. 6 for $S_f = 10$ and $B_e = 0.2$, $S_f = 0.1$ and $B_e = 0.4$, and $S_f = 1$ and $B_e = 0.2$, respectively. For all the presented cases, the shown geometries are taken from the last iteration, at $t_f = 200$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Distribution of the shear stress (top row) and of the growth stress (bottom row) at $t = 0$ and $t = 200$ with an eccentricity factor $M_f = 0.5$. The initial geometry ($t = 0$) corresponds to the result of the generation without eccentricity: in the left, the case with $S_f = 10$ and $B_e = 0.2$ was used (category I), and in the right, the case with $S_f = 0.1$ and $B_e = 0.4$ was used (category II), see Fig. 6. Additionally, the results on the left correspond to the generation under a normalised size of the growth region $D_g = 0.025$ and the ones in the right to $D_g \to \infty$. In the bottom row, $g$ refers to the growth term Eq. (9).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) $P_c^{\text{mean}}$ (—) and $P_c^{\text{max}}$ (---)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) $P_c^{\text{mean}}$ (—) and $P_c^{\text{max}}$ (---)</figcaption></figure>

Figure 10: Evolution of the contact pressure under a time-dependent eccentric load, with $M_f$ between -0.5 and 0.5. Fig. (a) shows the quality of the fields of (—) mean, $P_c^{\text{mean}}$, and (---) maximum, $P_c^{\text{max}}$, pressure by cycle. Fig. (b) depicts the field of mean pressure by cycle normalised by its maximum value during the first cycle, $p_{\text{ref}}^{\text{max}} = \max(P_c^{\text{mean}}(c_0))$; $c_0$ and $c_f$ in the colour bar refer respectively to the initial and the final cycle. Fig. (c) depicts the field of maximum pressure by cycle, $P_c^{\text{max}}$, normalised by its maximum value during the first cycle, $p_{\text{ref}}^{\text{max}} = \max(P_c^{\text{max}}(c_0))$. In this case, the category-III shape was used as the initial geometry and each eccentricity cycle was divided into 32 discrete steps. Lastly, Fig. (d) shows the fields of (—) mean, $P_c^{\text{mean}}$, and (---) maximum, $P_c^{\text{max}}$, pressure by cycle for the flat reference shape. Both fields are normalised by $p_{\text{ref}}^{\text{max}} = \max(P_c^{\text{max}})$, the maximum contact pressure during the cycle.
<figure><img src="image_15.png" /></figure>

Figure 11: Fields of contact pressure distribution by cycle using as initial geometry the one generated without eccentricity with $S_f = 1$ and $B_e = 0.2$ (category III) and with $D_g \to \infty$. Fig. (a) shows the normalised value of the maximum mean ($\max P_c^{\text{mean}}$) and maximum amplitude ($\max P_c^{\text{ampl}}$) of each cycle. Both fields are normalised such that the maximum mean and amplitude of the pressure is 1 at the first cycle. In the figure, the star indicates the initial cycle; the triangle indicates the cycle where $\max P_c^{\text{mean}}$ is minimum; the circle refers to the cycle where $\max P_c^{\text{mean}}$ is equivalent to the initial cycle but $\max P_c^{\text{ampl}}$ is lower; and the square refers to the cycle where $\max P_c^{\text{ampl}}$ reaches its minimum. Figs. (b) and (c) compare the contact pressure distribution of the initial cycle (star) with the cycles marked by the circle and the square, respectively.

<figure><img src="image_16.png" /></figure>

Figure 12: Three-dimensional scheme. Both geometries are obtained from a revolution of a two-dimensional contact interface. The revolution is made such that the mean diameter is $w$. $F_r$ and $F_a$ are respectively radial and axial forces applied in the xy-plane. A reference contact line which coincides with the xy-plane is shown in orange.
<figure><img src="image_17.png" /><figcaption>Figure 13: Pressure distribution under radial and axial loads in the three-dimensional framework. The geometry was obtained from a revolution of the contact pair generated with $B_e = 0.4$, $S_f = 0.1$ and $M_f = 0$ in the two-dimensional framework. Fig. (a) is a contour plot of the contact pressure distribution when $F_a = 0$ and $F_r$ is computed using Eq. (28). The red ◦ symbols represent the position of $F_r$. Figs. (b) and (c) depict the pressure distribution along the reference contact line where loads are applied, orange line in Fig. 12. In (b), $F_a = 0$ and in (c), $F_a$ grows until reaching the value of $F_r$. In both cases, the pressure is normalised by the maximum pressure when $F_a = 0$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 14: Pressure distribution under eccentric loads in the three-dimensional framework. The geometry was obtained from a revolution of the contact pair generated with $B_e = 0.4$, $S_f = 0.1$ and cyclic eccentricity, $M_f$ between -0.5 and 0.5, in the two-dimensional framework. Fig. (a) is a contour plot of the contact pressure showing the distribution when $M_f = 0$ and $M_f = 0.5$. $M_f$ determines the location of $F_r$ which is computed using Eq. (28) and is represented by the red ◦ symbols. The pressure is normalised by the maximum pressure when $M_f = 0.5$. Fig. (b) shows the fields of mean ($P_c^{\text{mean}}$) and maximum ($P_c^{\text{max}}$) pressure along the reference contact line where loads are applied, orange line in Fig. 12. The fields are normalised by the maximum contact pressure when $F_r$ is maximum and $M_f = 0$, $P_{\text{ref}}^{\text{max}} = \max(p_c|_{F_r^{\text{max}}, M_f=0})$.</figcaption></figure>