David Hernandez-Aristizabal
email: david.hernandez-aristizabal@univ-amu.fr

Santiago Arroyave-Tobón

Diego A Garzon-Alvarado

Jean-Marc Linares

Stress-adaptive design of 2D contact interfaces with uniform pressure: a bio-inspired approach

Keywords: bio-inspired design, generative design, joint morphogenesis, contact interface, stress adaptation

In this article, a bio-mimetic approach to the generative design of contact interfaces with uniform pressure distribution is explored. During the morphogenesis process, biological joints grow depending on the states of stress generating a shape adapted to the mechanical loads. This adaptation is driven by two main rules: shear stress inhibits growth, whereas cyclic hydrostatic compressive stress promotes it. In this work, we demonstrate that the stress-dependent growth rules of synovial joint morphogenesis can be applied to generate contact interfaces with uniform pressure in engineering applications. For that, we present a mathematical model that comprises a contact formulation and a bio-inspired growth function; the model is solved numerically using finite element methods. We analyse the impact of the growth rules of synovial joint development on the contact pressure distribution of two-dimensional contact interfaces. We study the parameter space of the bio-inspired growth function to fine-tune the model parameters. The model is tested in several cases with different boundary conditions and material properties. The results show that the proposed generative design process leads to contact interfaces that provide uniform contact pressure.

A quantitative metric of the uniformity of the contact pressure is also defined. This metric indicates that the bio-inspired process generates geometries as good as those obtained with other methods present in the literature without the need of an initial shape close to the final one. Thus, this work demonstrates that synovial joint morphogenesis can be adapted for generative design in engineering.

Introduction

Synovial joints are characterised by remarkable wear resistance. This performance is attributed to a complex structure that comprises tribological, material, and geometric properties [START_REF] Neville | Synovial joint lubrication -does nature teach more effective engineering lubrication strategies?[END_REF]. During morphogenesis, such a structure is formed by a growth process which is influenced not only by chemical but also by mechanical cues [START_REF] Galea | Making and shaping endochondral and intramembranous bones[END_REF][START_REF] Shea | The importance of foetal movement for co-ordinated cartilage and bone development in utero[END_REF][START_REF] Felsenthal | Mechanical regulation of musculoskeletal system development[END_REF]. For example, in chondrocytes and other cells, mechanical loads can activate integrins-a class of transmembrane proteins that attach to the extracellular matrix-which later activate intracellular signalling pathways that modify the cell activity [START_REF] Docheva | Integrin signaling in skeletal development and function[END_REF]. In fact, it has been shown that the absence of foetal movement leads to bone and cartilage abnormalities [START_REF] Shea | The importance of foetal movement for co-ordinated cartilage and bone development in utero[END_REF][START_REF] Felsenthal | Mechanical regulation of musculoskeletal system development[END_REF]. Hence, synovial joints adapt to mechanical loads as they grow.

This load adaptation is driven by two main rules. Cycles of compressive stress promote chondrocyte proliferation and cartilage repair and maintenance [START_REF] Sah | Biosynthetic response of cartilage explants to dynamic compression[END_REF][START_REF] Palmoski | Effects of static and cyclic compressive loading on articular cartilage plugs in vitro[END_REF]. In contrast, shear stress promotes chondrocyte ossification and cartilage degradation [START_REF] Beaupré | Mechanobiology in the development, maintenance, and degeneration of articular cartilage[END_REF][START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF]. This means that cyclic compressive stress promotes growth, while shear stress inhibits it and promotes material hardening. The accuracy of these rules has been proven computationally in several works [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF][START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF][START_REF] Sadeghian | Computational model of endochondral ossification: Simulating growth of a long bone[END_REF] where the authors reproduce synovial joint formation.

In mechanical engineering, joints are also common mechanisms for load transfer. They usually undergo dynamic contact loads which expose them to wear. In consequence, material is constantly lost on the surface and the lifetime of these mechanisms is shortened [START_REF] Tichy | Review of solid mechanics in tribology[END_REF]. Further, if the wear rate is not uniformly distributed, certain regions end up being overexposed to material lost and the failure occurs prematurely [START_REF] Kato | Classification of wear mechanisms/models[END_REF]. Hence, it is convenient to ensure a uniform wear rate.

Wear rate depends on the contact pressure, the material properties and the sliding speed of the bodies in contact [START_REF] Tichy | Review of solid mechanics in tribology[END_REF][START_REF] Popov | Wear[END_REF]. At imposed speeds, contact pressure uniformity is therefore a key factor. Nonetheless, the contact analysis is highly non-linear and presents unexpected phenomena.

For example, the contact between finite cylinders leads to infinite contact pressure at the edges-a phenomenon called edge effect.

In theory, this phenomenon can be reduced by crowning the cylinder along its axial direction.

Lundberg [START_REF] Lundberg | Elastische berührung zweier halbraüme[END_REF] investigated this phenomenon and found an analytical expression that leads to uniform contact pressure. Although it is inaccurate close to the edges [START_REF] Johnson | Contact mechanics[END_REF], it has been widely used in engineering and is an important reference profile. For example, it has been used in the shape design of bearing rollers [START_REF] Johns | Roller bearings under radial and eccentric loads[END_REF][START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF] and gear teeth [START_REF] Bergseth | Logarithmical Crowning for Spur Gears[END_REF]. Nonetheless, since Lundberg's profile is not adequate for many applications-for example when the contact interface does not comprise any flat surface-other strategies have been developed to address this problem, as described below.

With the aim of improving contact pressure distribution several approaches have been investigated. Among them, we can find developments in material properties [START_REF] Liu | Influence of material parameters on the contact pressure characteristics of a multi-disc clutch[END_REF][START_REF] Lin | An optimization design of contact interface material stiffness for improving the uniformity of the contact pressure[END_REF][START_REF] Zhou | Optimization of elastoplastic behavior of contact interface for improved contact stress distribution[END_REF], on topology [START_REF] Wei | FE analysis of a novel roller form: a deep end-cavity roller for roller-type bearings[END_REF][START_REF] Kristiansen | Topology optimization for compliance and contact pressure distribution in structural problems with friction[END_REF][START_REF] Lin | Design of component structure in assemblies for simultaneously regulating contact pressure distribution and natural frequencies[END_REF], on surface texture [START_REF] Tala-Ighil | Effects of surface texture on journal-bearing characteristics under steady-state operating conditions[END_REF][START_REF] Xu | Research on surface topography wear of textured work roll in cold rolling[END_REF][START_REF] Zhao | Effect of micro-textures on lubrication characteristics of spur gears under 3d line-contact EHL model[END_REF] and on shape [START_REF] Hilding | A computational methodology for shape optimization of structures in frictionless contact[END_REF][START_REF] Páczelt | Optimal shapes of contact interfaces due to sliding wear in the steady relative motion[END_REF][START_REF] Ou | A direct shape optimization approach for contact problems with boundary stress concentration[END_REF][START_REF] Najjari | Formula derived from particle swarm optimization (PSO) for optimum design of cylindrical roller profile under EHL regime[END_REF][START_REF] Duchemin | Optimization of contact profiles using super-ellipse[END_REF][START_REF] Sysaykeo | Clearance and design optimization of bio-inspired bearings under off-center load[END_REF]. These works are based on optimisation algorithms following either gradient-based or heuristic strategies. The former requires rigorous mathematical formulations and the computation of derivatives [START_REF] Upadhyay | Numerical analysis perspective in structural shape optimization: A review post 2000[END_REF]. Although they converge fast, the computation of the gradient may constrain the problem. For example, in [START_REF] Hilding | A computational methodology for shape optimization of structures in frictionless contact[END_REF] the shape of the contour in contact must remain constant, and the optimisation is achieved by modifying the shape of the contours that are not in contact. In addition, gradient-based algorithms are particularly sensitive to the initial guess since they can end up trapped in a local minimum [START_REF] Upadhyay | Numerical analysis perspective in structural shape optimization: A review post 2000[END_REF]. Concerning heuristic approaches, they can provide almost global solutions, but they usually require a more extensive evaluation of the design space [START_REF] Upadhyay | Numerical analysis perspective in structural shape optimization: A review post 2000[END_REF]. In addition, for the optimisation of the contact pressure distribution, it is sometimes necessary to impose the general shape of the profile, as done in [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF][START_REF] Duchemin | Optimization of contact profiles using super-ellipse[END_REF]. Further, when the design variable is the contact profile, the modifications are restricted to local variations [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF][START_REF] Ou | A direct shape optimization approach for contact problems with boundary stress concentration[END_REF][START_REF] Najjari | Formula derived from particle swarm optimization (PSO) for optimum design of cylindrical roller profile under EHL regime[END_REF][START_REF] Duchemin | Optimization of contact profiles using super-ellipse[END_REF][START_REF] Sysaykeo | Clearance and design optimization of bio-inspired bearings under off-center load[END_REF]. In consequence, the designer needs to have a certain knowledge of the form of the target shape. Another approach was proposed in [START_REF] Arroyave-Tobon | Generative design of joint contact surfaces inspired by biological morphogenesis[END_REF][START_REF] Marquez-Florez | From biological morphogenesis to engineering joint design: A bio-inspired algorithm[END_REF] where it was shown the feasibility of applying cartilage growth rules in engineering design. Nonetheless, they did not evaluate the uniformity of the contact pressure distribution and they focused on a particular case.

In this work, we demonstrate that the stress-dependent growth rules of synovial joint morphogenesis can be applied to generate contact interfaces with uniform pressure in engineering applications. We present a mathematical model that comprises a contact formulation and a bio-inspired growth function. We define a quantitative measure for pressure uniformity to evaluate the obtained contact interfaces. We study the impact of the growth rules of synovial joint development on the contact pressure distribution of two-dimensional contact interfaces. Then, we adapt the growth function for the generative design of contact interfaces and analyse the parameter space to fine-tune the model parameters. The model is tested in several cases with different boundary conditions and material properties.

The rest of the document is organised as follows. The bio-inspired growth model is presented in Section 2 and its numerical approximation in Section 3. In Section 4, we establish an analysis of the parameters to study the behaviour of the growth function and in Section 5 we set up examples to demonstrate the versatility of the methodology. Finally, in Section 6, we examine and discuss the results and in Section 7, we provide the conclusions of the study.

Bio-inspired stress-dependent growth function

Force balance equilibrium for bodies in contact

Synovial joints are formed in a growth process highly dependent on chemical and mechanical stimuli [START_REF] Galea | Making and shaping endochondral and intramembranous bones[END_REF]. As a result, they are adapted to both internal and external stimuli. To use this process for industrial applications, we first need to explore the growth forces that generate this adaptation.

In particular, the growth forces that are dependent on mechanical cues.

We shall begin with the properties of the desired design. Let us consider a set of k ∈ N continuous domains:

Ω i ∈ R 3 (where 1 ≤ i ≤ k) with closed boundaries: Γ i = ∂Ω i ∈ R 3 .
Let us also consider that these domains might be: in contact with each other, under external forces and constrained at given points. Thus, each

Γ i = Γ c i ∪ Γ f i ∪ Γ d i where Γ c i , Γ f i and Γ d i are
respectively: the contact boundaries, the boundaries where external forces are applied and the constrained boundaries. The force balance equation is then given by [START_REF] Wriggers | Contact kinematics[END_REF]:

∇ • σ i (x) = 0 ∀x ∈ Ω i , (1)
σ i (x) • n i (x) = f i (x) ∀x ∈ Γ f i , (2)
σ i (x) • n i (x) = p i (x) ∀x ∈ Γ c i , (3)
u i (x) = u 0 i (x) ∀x ∈ Γ d i , (4)
for 1 < i < k. σ i are stress tensors; n i are outward unit vectors; f i and p i are traction stress vectors related to external forces and contact pressures, respectively; and u 0 i are fields of fixed displacements.

Our objective is to define Ω i such that each p i has a uniform distribution-in other words, ∂p i /∂t i = 0 along Γ c i , where t i are tangent unit vectors to the surfaces Γ i : t i • n i = 0. To solve this problem, we will use the theory of synovial joint morphogenesis to obtain Ω i by a process of load adaptation.

Before continuing, let us establish some useful notation. The hydrostatic stress is defined as:

σ hyd = tr(σ)/3.
For the shear stress let us use the von Mises stress σ v = √ 3J 2 -where J 2 is the second invariant of the deviatoric stress tensor. It will also be useful to write the compressive stress as: σ hyc which is 0 when σ hyd > 0 and -σ hyd when σ hyd < 0.

The growth function

As we will use the theory of synovial joint morphogenesis to define Ω i , let us see the biological growth model. For cellular processes of growth, we can assume that the inertial forces are negligible with respect to viscous and elastic forces [START_REF] Lewis | Analysis of stable two-dimensional patterns in contractile cytogel[END_REF]. Then, the equilibrium of forces can be written as [START_REF] Lewis | Analysis of stable two-dimensional patterns in contractile cytogel[END_REF][START_REF] Stéphanou | A mathematical model for the dynamics of large membrane deformations of isolated fibroblasts[END_REF][START_REF] Murphy | A moving grid finite element method applied to a mechanobiochemical model for 3d cell migration[END_REF]:

∇ • σ i (x) + σ g (x) = 0 ∀x ∈ Ω i , (5)
where σ g is a stress tensor of growth forces. This tensor is then in charge of the shape adaptation to external loads. In synovial joint formation, this tensor is likely to promote oriented growth as chondrocytes and collagen fibres have different properties and orientations depending on the proximity to the contact region [START_REF] Ge | Osteoarthritis and therapy[END_REF]. Near the surface, they are oriented tangentially to the surface. Farther from the surface, they are oriented perpendicularly to the surface [START_REF] Ge | Osteoarthritis and therapy[END_REF]. In fact, it has been suggested that such preferential orientation might be related to the direction of principal strains [START_REF] Rakhsha | Simulation of surface strain in tibiofemoral cartilage during walking for the prediction of collagen fibre orientation[END_REF]. Thus, growth orientation not only influences the shape but clearly also the mechanical properties of the cartilage, particularly, its anisotropic nature. Nonetheless, in biology, some authors have considered a hydrostatic nature for the growth tensor in computational models and have obtained shapes similar to those seen in nature [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF]. Therefore, for this work, we will also assume hydrostatic growth stress for simplicity. We can then let σ g = σ g I where I is the identity matrix.

Research in joint formation indicates that, apart from biochemical cues, growth is controlled by compressive and shear stresses following the next two laws [START_REF] Sah | Biosynthetic response of cartilage explants to dynamic compression[END_REF][START_REF] Palmoski | Effects of static and cyclic compressive loading on articular cartilage plugs in vitro[END_REF][START_REF] Beaupré | Mechanobiology in the development, maintenance, and degeneration of articular cartilage[END_REF][START_REF] Carter | The role of mechanical loading histories in the development of diarthrodial joints[END_REF]: (i) proliferation takes places when chondrocytes are subjected to cyclic compressive stress and (ii) ossification is triggered by shear stress. If we neglect the difference in material properties before and after ossification, we can translate proliferation into growth and ossification into growth inhibition. In short, we can say that cyclic compressive stress prompts growth, while shear stress inhibits it.

To obtain the desired growth function, let us define different combinations of these two basic rules:

• Scenario I: expansion proportional to compressive stress. In this situation, we can study the effect of the first biological law in the absence of the second.

• Scenario II: contraction proportional to compressive stress. In this situation, we can see the behaviour of the growth process in opposition to the first biological law. This will shed light on how the contact interfaces driven by compressive stress stimuli can be.

• Scenario III: expansion proportional to compressive stress and inhibited by high shear stress. In this situation, we can test the evolution of the growth process following the two biological laws.

• Scenario IV: expansion proportional to compressive stress and inhibited by high shear or high compressive stresses. In this situation, we propose to enhance the inhibition of the growth. Inhibition will take place not only in the regions with high shear stress but also in regions with high compressive stress. As the distributions of shear and compressive stresses are not the same, this scenario will help us understand which regions should be discouraged from growing.

Since in biology, inhibition is usually imposed by threshold parameters, the following equation describes the proposed scenarios:

σ g (x) = α g σ hyc (x)H(τ lim -σ v (x))H(σ lim -σ hyc (x)), (6)
where α g is a parameter referring to the strength of the growth force, τ lim and σ lim are respectively the shear and the compressive stress thresholds, and H(ϕ) is the Heaviside function defined as:

H(ϕ) =              1 for ϕ ≥ 0, 0 for ϕ < 0. (7)
Although this type of equation has been successfully used in biology [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF], in [START_REF] Arroyave-Tobon | Generative design of joint contact surfaces inspired by biological morphogenesis[END_REF][START_REF] Marquez-Florez | From biological morphogenesis to engineering joint design: A bio-inspired algorithm[END_REF] it was indicated that a discontinuous σ g as that in Eq. ( 6) leads to discontinuous growth. Further, this discontinuity also affects the contact pressure distribution. To make Eq. ( 6) continuous, instead of the Heaviside function, let us use the Sigmoid function which is defined as:

S (ϕ) = 1 1 + exp (-υϕ) , ( 8
)
where υ is a positive parameter. yields:

σ g (x) = α g σ hyc (x)S (τ lim -σ v (x))S (σ lim -σ hyc (x)). ( 9
)
This function is now continuous; however, the stress magnitude is highly dependent on the magnitude of the external loads. As we want to define a system able to modify the contact shape independent of the stress magnitude, it is necessary to define the growth function in terms of normalised hydrostatic and shear stresses. We propose a normalisation based on the maximum stress values: let σ * hyc (x) and σ * v (x) be normalised compressive and shear stresses:

σ * hyc (x) = σ hyc (x) max(σ hyc ) , σ * v (x) = σ v (x) max(σ v ) . As σ * hyc (x), σ * v (x) ∈ [0, 1],
S(0.5 -ϕ, v) v → 0 v = 1 v = 10 v = 100 v → ∞
ϕ ref = λ ϕ ref = 0 ϕ ref = 0.5 ϕ ref = 1.0 ϕ ref = γ (b) Sigmoid function varying ϕ ref .
Figure 1: Sigmoid function response. As υ approaches zero, S becomes the constant function 0.5. In contrast, as υ approaches infinity, S becomes the Heaviside function. Further, for a given υ, there exists a minimum value γ such that for all

ϕ ref ≥ γ, S (ϕ ref -ϕ) ≈ 1 for all ϕ ∈ [0, 1].
Similarly, there also exists a maximum value λ such that for all

ϕ ref ≤ λ, S (ϕ ref -ϕ) ≈ 0 for all ϕ ∈ [0, 1].
equation can then be rewritten as:

σ g (x) = α g σ * hyc (x)S τ * lim -σ * v (x) S σ * lim -σ * hyc (x) β g , (10)
where τ * lim and σ * lim are now dimensionless parameters. β g has units of pressure to make Eq. ( 10) compatible with Eq. ( 5), for this work β g = 1GPa. :::: The :::::::: strength ::: of ::: σ g :: is ::::: now ::::: given ::: by :::: the :::::::: product ::

of ::: β g :::: and :::: α g . :::::: Thus, ::: we :::: can :::: let ::: α g ::: be : a ::::::::::::::: dimensionless :::::::::: parameter :::: that ::::::::: provides :::: the :::::::::: magnitude ::::: and ::

β g :::: the :::::::::: parameter ::::: that :::::::: contains :::: the :::::: units, ::: in :::: this ::::: case :::::::::::: β g = 1GPa. : The parameter space for each scenario is shown in Table 1.

Table 1: Parameter space for each scenario of study. In scenario I, we consider uninhibited expansion proportional to hydrostatic compressive stress; in scenario II, uninhibited contraction proportional to hydrostatic compressive stress; in scenario III, expansion proportional to hydrostatic compressive stress and inhibited by high shear stress; and in scenario IV, expansion proportional to hydrostatic compressive stress and inhibited by both high shear and high hydrostatic compressive stresses.

Scenario

α g τ * lim σ * lim I (0, +∞) [γ, +∞) [γ, +∞) II (-∞, 0) [γ, +∞) [γ, +∞) III (0, +∞) (0, γ) [γ, +∞) IV (0, +∞) (0, γ) (0, γ)

The material constitutive law

As the purpose of this work is to study adaptation in an engineering context, it is reasonable to consider a linear elastic model. Therefore, the stress-strain relation is given by:

σ(x) = E 1 + ν ε(x) + ν 1 -2ν trε(x) , (11)
where E and ν are Young's modulus and Poisson's ratio, respectively, and

ε(x) = 0.5 ∇u(x) + (∇u(x))
is the linear strain tensor.

It would be possible to use the properties of biological materials. Although this would be closer to synovial joint development, the adaptation might be biased by stress fields that do not represent the engineering context. In addition, even in biology, linear elastic models have been used to reproduce synovial joint and bone formation [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF][START_REF] Giorgi | Mechanobiological simulations of prenatal joint morphogenesis[END_REF][START_REF] Sadeghian | Computational model of endochondral ossification: Simulating growth of a long bone[END_REF].

A bio-inspired shape-design methodology

We propose an iterative design process that allows us to replicate cyclic compressive stresses.

Each cycle is divided into a stress computation step and a growth step. For the stress computation, we need to solve Eqs. ( 1) to (4)and, for the growth step, Eqs. ( 4) and ( 5) together with the following boundary condition:

u i (x) • n i (x) = 0, ∀x ∈ Γ f i . (12)
In other words, during the growth step, the displacement along Γ f i is purely tangential. Then, we redefine Ω i by applying the computed displacements.

To objectively compare the response of the algorithm using different values of the parameters in Table 1, let us consider a time-like variable, tm . As it can be seen in Eq. ( 10), σ g is proportional to α g . This means that the velocity of growth is also proportional to α g . Therefore, we consider:

tm = m|α g |, (13)
where m ∈ N is the number of a given iteration, as a time-like measure. For example, if we have two cases, one with α g = 1 and the other with α g = 10, one iteration of the latter is equivalent to ten iterations of the former.

Algorithm 1 Bio-inspired shape-design algorithm. Input: Geometry: Ω i (0) for 1 ≤ i ≤ k. Input: Boundary conditions: f i and u 0 i . Input: Material properties: E i and ν i . Input: Growth function parameters: α g , τ * lim and σ * lim . Input: Total simulation time: tf . Output: Geometry adapted to the loading conditions:

Ω i (t f ). 1: m ← 1 2: tm ← m|α g | 3: while tm < tf do 4:
Compute the stress state that stimulates growth solving Eqs. ( 1) to (4).

5:

Compute the growth stress using Eq. ( 10).

6:

Find the displacement field u(x( tm )) that satisfies Eqs. ( 4), ( 5) and (12).

7:

x( tm+1 ) ← x( tm ) + u(x( tm ))

8:

m ← m + 1 9:
tm ← m|α g | 10: end while Additionally, Eq. ( 10) also depends on the parameters of inhibition, υ, τ * lim and σ * lim . Therefore, if the inhibition region is larger with one set of parameters than with another, the geometry will change more slowly with the former set than with the latter. In order to compare the model response with different parameters, we will make the L 2 -norm of the growth function equal to |α g | at each iteration. In other words, we will use σ α instead of σ g defined as:

σ α = |α g | ||σ g || L 2 (∪ i Ω i ) σ g (14)
where

|| • || L 2 (∪ i Ω i ) refers to the L 2 -norm in ∪ i Ω i .
The use of Eqs. ( 13) and ( 14) allows us to objectively compare the response of the algorithm using different parameters.

Algorithm 1 presents the algorithm to mimic the biological growth process. As can be seen, computing the displacement field (steps 4-6) is quite simple compered to gradient-based algorithmsthere are no derivatives to calculate nor surface-dependent variables. The implementation was done using Code Aster [START_REF] De | Finite element code aster, analysis of structures and thermomechanics for studies and research[END_REF] and GMSH [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and postprocessing facilities[END_REF]. The former was used to solve the equations by finite element methods and the latter was used to control the evolution of the mesh. The computational implementation is available in: https://github.com/sarroyavet/BioDesign_ANR/ tree/main/StressAdaptiveDesignOf2D

Impact of the parameters on the contact interface

In this section, we will set up a reference case of study to analyse the impacts of the model parameters on the uniformity of the contact pressure. In addition, we will define metrics to measure the performance of the algorithm and the generated contact interfaces. As we will explore different growth scenarios and parameters, we will consider a two-dimensional framework and frictionless contact. In the context of contact analysis, the theory of semi-infinite elastic bodies is widely used [START_REF] Johnson | Contact mechanics[END_REF]. In consequence, we will consider the contact between a half-space and an infinite cylinder that satisfies the state of plane-strain. Further, the consideration of frictional effects significantly augments the computational time. For example, in [START_REF] Mattei | Frictionless vs. Frictional Contact in Numerical Wear Predictions of Conformal and Non-conformal Sliding Couplings[END_REF], the authors showed that wear simulations considering friction were up to three times longer than their frictionless counterparts.

Even though we are not modelling wear, the process is similar: a calculation of the contact pressure followed by a geometry update. Thus, using a two-dimensional frictionless contact framework will significantly reduce computational time.

We will also limit the growth to only one domain. Thus, the degrees of freedom of the design problem will be reduced to one domain, facilitating the analysis of the adaptation process. This assumption is also reasonable in the engineering context as it would be appropriate to define an easy-to-manufacture geometry and only generate an adapted and complex shape. For example, in the case of roller bearings, it is common to crown the roller profile while letting flat that of the race. Nonetheless, it would be interesting to explore the advantages of allowing both domains grow in future works.

To sum up, in this section, the model and the algorithm will be simplified by the following assumptions: (i) two-dimensional framework where plane-strain conditions are applied, (ii) frictionless contact and (iii) growth applied to only one domain.

Description of the case of study

Let us consider the contact between two elastic semi-infinite bodies: a cylinder (Ω 1 ) and a half-space (Ω 2 ). Let us also apply a load on the cylinder in the normal direction to the surface of the half-space. Since this case satisfies the plane-strain conditions, it can be reduced to a twodimensional framework, a schematic view is depicted in Fig. 2a. To analyse the adaptation of one

Ω 1 (0) Ω 1 (0) Ω 2 F p (a) Initial geometry. F Ω 1 (t f ) Ω 2 p Ω 1 (t f ) (b) Final geometry.
Figure 2: Two-dimensional setting of the case of study. Here p refers to the contact pressure along the contact interface and F to the applied force.

Table 2: Geometric parameters of the initial geometries. R, l i are illustrated in Fig. 3 for both types of geometry.

Geometry In other words, let σ g (x) = 0 for all x ∈ Ω 2 .

R[mm] l 1 [mm] l 2 [mm] l 3 [mm] l 4 [mm]
The goal of the growth function is to yield a shape to obtain a constant contact pressure. Thus, after the design process, we expect to have a geometry adapted to the load that generates a uniform contact pressure as depicted in Fig. 2b.

Geometry, material properties and loading conditions

First, let Young's modulus and Poisson's ratio be 210GPa and 0.3, respectively, and let F = 10kN. Regarding the initial cross section of Ω 1 , we can define two types depending on the distribution of the initial contact pressure: (i) a higher pressure in the middle of the contact area or (ii) a higher contact pressure at the edges of the contact area. We can use a circle for the former and a rectangle with rounded edges for the latter. Fig. 3 shows a schematic representation of the geometry and Table 2 establishes the values of the dimensions. The load, the material properties and the initial size of the geometries were adapted from the case of study presented in [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF].

Performance evaluation

To evaluate the performance of the procedure, let us define a measure to qualify the contact pressure distribution. In the best situation, the contact pressure is uniformly distributed and, therefore, all the contact area is properly exploited. Thus, let us define a measure to compare the pressure distribution of a given contact area with an equivalent uniform pressure. Such a measure should be able to tell how close a given pressure distribution is to a uniform distribution. For each tm , let a c ( tm ) be the contact area and let p e (a c ( tm )) be an equivalent constant pressure given by:

p e (a c ( tm )) = F a c ( tm ) . (15)
Now, let us define the pressure distribution quality as:

Q p ( tm ) = 1 - ||p e (a c ( tm )) -p( tm )|| L 2 (Γ c 1 (t m )) ||p e (a c ( tm ))|| L 2 (Γ c 1 (t m )) + ||p( tm )|| L 2 (Γ c 1 (t m )) , (16)
where

|| • || L 2 (Γ c 1 (t m )) refers to the L 2 -norm in Γ c 1 ( tm )
. By means of the triangle inequality [START_REF] Axler | Hilbert spaces[END_REF], we know that 0 ≤ Q p ≤ 1. Additionally, Q p → 1 as p → p e and Q p → 0 as p concentrates at a specific point. We then look for high values of Q p . Table 3 indicates the value of Q p based on the reconstruction of results presented in other works.

Apart from this evaluation at each step, we can also compute the average of Q p as well as the [18] Fig. 4 0.93 Optimisation and approximation of a logarithmic profile. [START_REF] Duchemin | Optimization of contact profiles using super-ellipse[END_REF] Fig. 11 0.95 Optimisation of a super-elliptic profile. [START_REF] Ou | A direct shape optimization approach for contact problems with boundary stress concentration[END_REF] Fig. 15 0.95 Optimisation of the contact profile. The present work 0.96 Bio-inspired generative design.

[19] Fig. 7 0.99 Optimisation of a logarithmic profile.

Other design variables.

[30] Fig. 6 0.86 Optimisation of the non-in-contact contour.

[24] Fig. 9 0.94 Optimisation of a predefined topology. [START_REF] Zhou | Optimization of elastoplastic behavior of contact interface for improved contact stress distribution[END_REF] Fig. 5a 0.94 Optimisation of the material properties. [START_REF] Lin | Design of component structure in assemblies for simultaneously regulating contact pressure distribution and natural frequencies[END_REF] Fig. 6b 0.99 Topology optimisation.

average of the maximum contact pressure along the simulation as:

Qp = m=m f m=1 Q p ( tm ) m f , pmax = m=m f m=1 p max ( tm ) m f , with t f = |α g |m f ,
where m f refers to the last iteration and t f to the total time. Thus, we can evaluate the consistency of the design process in a general manner.

Exploration of the design process

As we aim to understand how the growth rules of joint formation can be applied to our case of study, we propose the following pathway. First, we will see the impact of α g in terms of convergence in both scenarios I and II. From this study, we will define an appropriate value of α g .

Then, we will explore the response of the design process in scenarios III and IV and qualify the contact pressure distribution varying τ * lim , σ * lim and υ. Following this exploration, we will select an appropriate combination of parameters and examine the evolution of the design process utilising a rectangle as the initial geometry.

Examples of application

In this section, we will present other examples to evaluate the applicability of the bio-inspired growth model after the definition of the parameters. We will consider plane-stress and axisymmet-

F Ω 1 (0) Ω 2 R l n Figure 4
: Initial geometry of the axisymmetric case of study. The coordinate system l and n will be used later to define the surface of Ω 1 relative to the surface of Ω 2 .

ric conditions and we will change the material properties, the magnitude of the load and the fixed geometry (Ω 2 ). Further, we will also consider an additional rule to satisfy an additional design requirement.

Plane-stress and softer material

For this test, let us consider the contact between two thin plates. We shall consider the analogous plane-stress case to the plane-strain case presented in Section 4.1 (see Fig. 3). Nonetheless, let us consider a softer material (for both Ω 1 and Ω 2 ) with E = 2GPa and ν = 0.4-these elastic properties are close to those of some engineering polymers-, R = 10mm and F = 100N.

Axial symmetry and geometries with distinct material properties

In this case, let us consider the contact between a soft material for Ω 1 (E = 2GPa and ν = 0.4) and a rigid material for Ω 2 (E = 250GPa and ν = 0.2) and the initial geometry presented in Fig. 4 with R = 14mm and F = 3000N. This case represents the loading conditions and material properties of hip-joint prostheses under static analyses [START_REF] Ruggiero | Implementation of a finite element deformation model within an elastohydrodynamic lubrication numerical solver for a ball in socket tribopair[END_REF][START_REF] Girard | Femoral head diameter considerations for primary total hip arthroplasty[END_REF].

Additional design requirement: maximum von Mises stress

So far, the design process makes the domain grow indefinitely. This means that we need to include an additional requirement to constrain the growth. It could be a geometric restriction (such as a maximum contact area) or a structural constraint (such as a maximum shear stress).

We can combine the growth function in Eq. ( 10)-whose main objective is to provide a uniform contact pressure-with a restriction function. Let us then rewrite Eq. ( 5) as:

∇ • σ(x) + η g σ g (x) + η r σ r (x) = 0, (17)
where η g and η r are the strengths of the contact pressure improvement and of the restriction, respectively. We can use Q p to define η g as:

η g = 1 -Q p . (18)
For the definition of σ r , let us impose a maximum von Mises stress. First, we can define a von Mises stress quality as:

Q v = 1 - max(σ v (x)) -τ ref max(σ v (x)) + τ ref , (19)
where τ ref is the target von Mises stress of the restriction. As it can be seen,

Q v → 1 as max(σ v ) → τ ref , Q v < 1 when max(σ v ) > τ ref and Q v > 1 when max(σ v ) < τ ref .
We can also consider σ r = I (hydrostatic) and

η r = α r (1 -Q v ), (20)
where α r is a positive parameter. As we can see, σ r is space-independent and produces global expansion or contraction. For this case, we will consider the plane-strain case presented in Section 4.

Results and discussion

Convergence

To ensure numerical precision, the contact algorithm was tested with different contact element lengths. For that, the contact pressure was calculated using the circular geometry described above.

This case can be analytically solved using Hertz theory [START_REF] Johnson | Contact mechanics[END_REF] in cases where the cylinder diameter is significantly larger than the contact area width. In consequence, we applied a load of 100N so that, according to Hertz theory, the maximum contact pressure is 903.4MPa, and the contact area width is 0.141 [mm]. Thus R a c and the conditions for Hertz theory are satisfied.

) = Γ c f (x) 2 dΓ c / Γ c dΓ c and in (b) e Γ c ×[0,t f ] ( f ) = t f 0 Γ c f (x, t) 2 dΓ c dt/ t f 0 Γ c dΓ c
dt are the L 2 -norms of a given function f . In (a) f = (p Hertzp FEM )/p FEM max ) which is the normalised error between the contact pressure calculated by Hertz theory and by the finite element method. In (b) f = (p α/2p α )/p α/2 max (t) which is the normalised contact pressure error using as reference a simulation using α/2. Fig. 5a shows the L 2 -norm of the normalised difference between the results using Hertz theory and finite elements varying the element length. The normalisation was made with respect to the maximum contact pressure calculated by Hertz theory. It can be seen that as the mesh is refined, the numerical approximation approaches the analytical solution. In addition, for this case, an element length of 0.01mm is sufficient to make the error norm lower than 0.05. Thus, this value was selected for the other simulations.

As indicated in Eq. ( 13), α g works as a time step variable. Therefore, the convergence of the algorithm with α g was also tested. For that, t f was set to 500-which is large enough to see the behaviour of the design process and to test the convergence. The L 2 -norm of the normalised contact pressure between simulations was computed. The normalisation was made by means of the maximum contact pressure at each step. Fig. 5b shows the relative difference among several simulations varying α g from 0.3125 to 10 for each scenario. We can see that for scenarios II, III and IV, as α g decreases, such a difference also decreases. Thus, for these scenarios the algorithm is also convergent in terms of α g -we will take a closer look for scenario I in the following section.

We selected α g = 1.25 to study the response of the design process in scenarios III and IV as, for such a value, the error norm is lower than 0.05. For scenarios I and II, we used the simulations with α g = 0.3125, the lowest value used in the convergence analysis.

The response of the design process in scenarios I and II

Fig. 6a shows the evolution of the contact pressure distribution (from blue to red lines) for scenario I. For this case, it was necessary to add a condition to the mesh such that the number of nodes in contact was always above a given threshold-we used a minimum of 50 nodes. This additional condition was necessary since the contact area was reduced during the simulation. We can see that growing (from blue to red distributions) proportionally to the compressive stress makes a pressure distribution highly concentrated on the initial point of maximum contact pressure. This is a reasonable result given the fact that the maximum compressive stress is located on the surfacemore precisely, at the point of maximum contact pressure. Thus, this region grows faster than its neighbours enlarging the curvature. In addition, this curvature increase reduces the contact area and, in consequence, the maximum contact pressure rises. We can say that in the scenario I, the shape tends to have an infinite contact pressure and a zero contact area. This singularity explains the behaviour of the convergence seen in Fig. 5b.

On the other hand, Fig. 6b shows that contraction (from blue to red distributions) in proportion to compressive stress makes the pressure concentrate at the edges. As the pressure changes from being concentrated in the middle to the edges, the profile flattens and the maximum pressure lowers. However, once the edge effect is formed, the maximum pressure rises . :::

and : Additionally, it can also be seen ::: the ::::::::: pressure :::::::::::: distribution :::::: tends :: to ::::::::: become :::::::: singular ::: at :::: the :::::: edge. ::::: Yet, ::: in :::: the :::: used ::::::::::: timescale, :::: this ::::::::: tendency ::::: does :::: not :::::: affect ::: the ::::::::::::: convergence :::: with :::::::: respect :: to :::: α g , :: as :::: can ::: be ::::: seen ::: in ::::::: Fig. 5b. : :::::: Lastly, :::: we :::: can :::: also :::: see that this contraction also enlarges the contact area.

Response of the design process in scenarios III and IV

We now know that while growing proportionally to compressive stress decreases the contact area and concentrates the pressure at the middle, contraction proportional to compressive stress increases the contact area and concentrates the pressure at the edges. The former generates geometries that tend to have an infinite curvature at the contact zone, while the latter makes geometries which tend to be flat (zero curvature). Since both scenarios are not desirable, we continued the exploration of the growth rules. constant function, Eq. ( 10) reduces to σ g (x) = cσ * hyc (x), where c is a constant, and no inhibition is present. In fact, we can see in Fig. 7b that the region of higher Qp is also the region of higher contact pressure. This is exactly what happened in the uninhibited growth scenario.

We can also see that, in scenario III, we can not obtain an appropriate result. In both graphs of Fig. 7, the red regions specify the desired results-high pressure quality and low contact pressure.

As can be seen, these regions exclude each other-the red region in Fig. 7a intersects with the blue region in Fig. 7b and the other way around-; this justifies the need to modify the growth function.

To address this situation, an additional inhibition was added. The results in scenario I indicate that unrestricted growth proportional to the compressive stress yields extremely high contact pressure. Since the maximum shear stress is not located in the same region as the maximum compressive stress, the region of maximum compressive stress is not always necessarily inhibited.

Therefore, the growth force keeps enlarging the curvature around the point of maximum contact pressure. Thus, in scenario IV, the growth force was inhibited by high compressive stress. As a consequence, two new parameters appear σ * lim and another υ. Yet, we let σ * lim = τ * lim and let the same υ for both Sigmoid functions of Eq. [START_REF] Carrera-Pinzón | Computational model of a synovial joint morphogenesis[END_REF]. Fig. 8 shows the response of the design process with the new growth function modifying υ and the stress limits. In Fig. 8a, we see that Qp is significantly greater (up to 0.975) than in the previous scenario. In addition, the zone of higher pressure quality intersects the zone of lower maximum contact pressure-for σ * lim = τ * lim = 0.5 and υ > 5. Hence, we see that the added inhibition leads to the desired contact pressure distribution.

Fig. 8 also indicates the best cases by (they correspond to υ = 6.4, 7.1, 7.8, 8.5). In addition, in Fig. 8a, there is a second local maximum value around υ = 3.6 indicated by . In order to better understand this situation, Fig. 9 shows the evolution of the contact pressure distribution with τ * lim = σ * lim = 0.5 and υ = 3.6 and υ = 7.1. In both cases, as the iterations go on, the contact area augments, the maximum contact pressure decreases and the contact pressure distribution flattens.

In addition, the pressure quality is high-above 0.96 for υ = 3.6 and above 0.97 for υ = 7.1.

Nonetheless, we can notice an important difference. For υ = 3.6, the maximum pressure remains in the centre, while for υ = 7.1, it goes to the edges. Thus, the combined inhibition of shear and compressive stresses balances the expansion and contraction present in scenarios I and II. With υ = 3.6, the effects in scenario I slightly dominate, while with υ = 7.1 the effects in scenario II are somewhat more present. This type of behaviour can also be seen in logarithmic profiles. For example, Cui and He [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF] presented a modified Lundberg profile with similar contact pressure distributions. Depending on the combination of parameters, the pressure tends to be concentrated either in the middle or at the edges.

Geometric properties of the contact interface

In terms of geometry, we can see both the profile and the curvature of the last shape obtained with υ = 7.1 in Fig. 10. Fig. 10a shows that the algorithm generates a profile whose slope increases slowly from the centre to a region close to the end of the contact area. Then, it rises rapidly. This can be more easily seen in Fig. 10b where the curvature is plotted. There is a peak of maximum curvature around the end of the contact area-the edge of the contact interface-that suggests the presence of a sharp edge. This can be related to the fact that the junction between the adapted profile and a residual of the initial profile is at a c . This figure also shows the profiles and curvatures of logarithmic profile (constructed based on Lundberg's equations [18-20, 51, 52]) and of a circular profile. Again, we can see similarities between our profile and logarithmic profiles. As we can see, the Lundberg profile curvature is also low at the middle and with a sharp peak at the edges. The contact interface using υ = 3.6 also presents this affinity.

Influence of the initial geometry

To further investigate the response of the design process, the initial geometry was changed.

More precisely, the initial contact pressure distribution. Fig. 11 shows the evolution of the contact geometry, in reality, the influence of the growth force is the same. For υ = 3.6, the point of initial maximum pressure defines the region where the pressure will be concentrated. Thus, a protuberance around this region appears. In contrast, since for υ = 7.1 the edge effect slightly dominates, the domain grows as in Fig. 9b. Therefore, for applications where a uniform contact pressure is desirable, it is better to set the initial geometry such that the initial contact pressure distribution has its maximum centred.

Change of material properties and deformation model

To study the applicability of the growth function, we defined other cases of study modifying the material properties, the shape of Ω 2 , the magnitude of the load, and the two-dimensional model.

We tested such modifications using the growth function with σ * lim = τ * lim = 0.5 and υ = 7.1.

For the plane-stress case, in Fig. 12, we can see the evolution of the contact pressure with σ * lim = τ * lim = 0.5 and υ = 7.1. As we can see, the contact area is enlarged while the contact pressure becomes uniform, reaching a value of pressure distribution quality of Q p = 0.98. Similarly to the previous case (Fig. 9b), the maximum pressure is located at the edge. With respect to the axisymmetric case, in Fig. 13, we can see the results of the growth process with σ * lim = τ * lim = 0.5 and υ = 7.1. Fig. 13a shows the shape of Ω 1 in the coordinate system defined by l and n drawn in Fig. 4. Interestingly, the profile is also similar to Lundberg's profile as can be seen in Fig. 10a. In addition, as seen in Fig. 13b, the contact area was also enlarged, and the contact pressure distribution was improved. In addition, we obtained a high value of pressure distribution quality at the last iteration, Q p = 0.97.

Although the number of examples presented is not enough to definitely conclude about the response of the model, they show a tendency to generate contact interfaces with uniform pressure.

This is an important insight since we consider that the efficacy of the bio-inspired generative design cannot be sensitive to the parameters of the growth function. In other words, we consider that the response of the design process seen in Fig. 8 is representative for situations different to the one described in Section 4.

In terms of material properties, the model can be extended to improve the growth phenomenon or to use more complex materials. For the former, it is possible to add a viscous stress tensor similar to the one present in [START_REF] Murphy | A moving grid finite element method applied to a mechanobiochemical model for 3d cell migration[END_REF] to model cell migration. Such a tensor would directly introduce time dependence. Nonetheless, it could represent an unnecessary cost for a methodology to generate shapes adapted to their loading environments. Further, additional parameters-shear and bulk viscosities-would be introduced. Regarding more complex materials, it would be interesting to study the process response adding plasticity. In addition, in line with bio-inspiration from synovial joints and considering the increasing development of cartilage-inspired materials with poroelastic properties [START_REF] Raske | Finite element investigations of the fluidsolid behaviour in a bio-inspired poroelastic bearing[END_REF], the use of poroelastic constitutive laws may be of interest.

Response of the design process after adding another design requirement

In the previous results, since there were no restrictions for the design process, only the evolution of the contact pressure distribution was evaluated regardless of the profile length. In contrast, in this section a restriction to the maximum shear stress is also present; thus, the contact length is fixed. :::::: length ::: of ::: the :::::::: contact ::::::: profile ::::::: reaches :: a :::::: stable ::::::: value. the restriction and makes the contact pressure uniform. This result demonstrates that the process of synovial joint morphogenesis can be adapted to engineering applications Nonetheless, the shape of the contact pressure is slightly different (see Figs. 9b and14d).

Q p Q v (c) (d)
The edge effect is rounded by applying the restriction term. This means that σ r also affects the distribution of the contact pressure. In addition, σ g also modifies the contact area. Thus, it would be interesting to separate the role of each term. In other words, it would be convenient to define σ g and σ r such that the former would only affect the contact pressure distribution while the latter would only modify the contact area.

Comparison with other works

The results presented show that the proposed model generates high quality contact interfaces (in terms of contact pressure distribution). We obtained values of Q p between 0.96 and 0.98. As can be seen in Table 3, only the works presented in [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF][START_REF] Lin | Design of component structure in assemblies for simultaneously regulating contact pressure distribution and natural frequencies[END_REF] are superior in terms of Q p (see Table 3). Nonetheless, in [START_REF] Cui | A new logarithmic profile model and optimization design of cylindrical roller bearing[END_REF], the methodology is defined only for contact against flat geometries.

The work presented in [START_REF] Lin | Design of component structure in assemblies for simultaneously regulating contact pressure distribution and natural frequencies[END_REF] can be applied to curved geometries. However, its design variable is not the shape of the contact interface but the material distribution; hence, the potential engineering applications are not the same than those of the present work.

Conclusion

In this work, we studied the impact of applying the rules of synovial joint morphogenesis on the generative design of contact interfaces. For that, we implemented them in an algorithm that mimics synovial joint growth under engineering conditions in an iterative manner. The output of the generative design process is a contact interface adapted to the loading conditions.

The obtained contact interfaces provided a uniform pressure distribution. Such interfaces were achieved following this stress-adaptive rule: growth is proportional to compressive stress, but it is inhibited at regions with high shear stress or high compressive stress. The results showed that the rule works well under different two-dimensional situations (plane-strain, plane-stress and axisymmetric conditions), with different material properties and under different loading conditions.

Further, the results also showed that the adaptation process can be applied against curved surfaces.

Furthermore, it is possible to control the size of the contact interface by an additional rule that promotes a global expansion (contraction) when the maximum shear stress is higher (lower) than a given requirement.

The proposed model differs from those present in the literature since it does not require the initial geometry to be close to the final solution, which is an important characteristic in terms of generative design.

Nonetheless, we acknowledge that our work has some limitations. The present work focuses on the contact interface and does not deal with the shape of the free boundary. In addition, as, so far, the current algorithm has only been implemented for frictionless and linear elastic conditions, we cannot conclude about the direct applicability of the model outside these conditions.

Future work is required to continue developing the proposed generative design process. It is necessary to determine appropriate criteria to define the shape of the free boundaries. It would be interesting to study the behaviour of the bio-inspired growth process in applications with dynamic loads or with friction and lubrication. Additionally, more biological processes could be mimicked.

For example, it might be interesting to explore the advantages of mimicking the hardening effect during the ossification process or the definition of an anisotropic growth tensor-for instance, one taking into account the direction of the principal strains.

Finally, the results presented here demonstrate that synovial joint morphogenesis can be adapted to generative design in engineering.

Fig. 1

 1 shows the behaviour of the S (ϕ ref -ϕ) varying υ and an offset along the ϕ-axis, ϕ ref , for ϕ ∈ [0, 1]. As it can be seen in Fig. 1a, as υ approaches zero, S becomes the constant function 0.5. In contrast, as υ approaches infinity, S becomes the Heaviside function. Further, Fig. 1b shows that: for a given υ, there exists a minimum value γ such that for all ϕ ref ≥ γ, S (ϕ ref -ϕ) ≈ 1 for all ϕ ∈ [0, 1]. Similarly, there also exists a maximum value λ such that for all ϕ ref ≤ λ, S (ϕ ref -ϕ) ≈ 0 for all ϕ ∈ [0, 1]. Hence, changing the H(•) by S (•) in Eq. (6)

  the adaptation becomes independent of the stress magnitude-in other words, it only depends on the stress distribution. Using the normalised stresses, the growth force

  (a) Sigmoid function varying υ.

  external stimuli untouched, let now the growth be only applied to the cylinder.

Ω 1

 1

Figure 3 :

 3 Figure 3: Geometry of the case of study.

  Hertz max ) (a) Convergence in h.

  Convegence in α g with t f = 500.

Figure 5 :

 5 Figure 5: Convergence of the model in mesh size and equivalent time step. In (a) e Γ c ( f ) = Γ c f (x) 2 dΓ c / Γ c dΓ c and in (b) e Γ c ×[0,t f ] ( f ) =

Figure 6 :

 6 Figure 6: Evolution (from blue to red lines) of the contact pressure along the contact line in scenarios I (a) and II (b) with |α g | = 0.3125. l refers to arc length and the red shapes refer to the geometry of last iteration.

Figure 7 :

 7 Figure 7: Response of the design process in scenario III with t f = 2000 in terms of the average of the pressure quality (a) and the average of the maximum contact pressure (b) along the simulation.

Fig. 7 Figure 8 :

 78 Fig. 7 shows the response of the design process in scenario III. Fig. 7a indicates the average of the pressure quality and Fig. 7b the average of the maximum contact pressure along the simulation. It can be seen that the pressure quality is higher (up to 0.855) at lower values of v or higher values of τ * lim . However, such a combination almost eliminates the inhibition. As seen in Fig. 1, low values of υ and high values of τ * lim cause S to be almost a constant function. If S becomes a

Figure 9 :

 9 Figure 9: Evolution (from blue to red lines) of the contact pressure along the contact line in scenario IV with τ * lim = σ * lim = 0.5 and de circular initial geometry. υ = 3.6 in (a) and υ = 7.1 in (b). l refers to arc length and the red shapes refer to the geometry of the last iteration.

Figure 10 :

 10 Figure 10: Geometric characteristics of the last profile obtained with υ = 7.1.

Figure 11 :

 11 Figure 11: Evolution (from blue to red lines) of the contact pressure along the contact line in scenario IV with τ * lim = σ * lim = 0.5 and the rectangular initial geometry. υ = 3.6 in (a) and υ = 7.1 in (b). l refers to arc length and the red shapes refer to the geometry of the last iteration.

Figure 12 :

 12 Figure 12: Evolution (from blue to red lines) of the contact pressure along the contact line for the plane-stress case with υ = 7.1. l refers to arc length.

Figure 13 :

 13 Figure13: Results of the axisymmetric case with υ = 3.6 after the growth process. (a) shows the contact interface of the last geometry in the coordinate system l and n seen in Fig.4 and (b)shows the evolution (from blue to red) of the contact pressure (l refers to arc length). The red shape refers to the geometry of the last iteration.

  Fig. 14d indicates the distribution of the contact pressure along the design iterations.As we can see, the contact interface changes rapidly at the beginning and, once Q p and Q v are close to 1, the evolution slows down. For this shape: max(|p c |) = 1.749GPa, max(σ v ) = 0.999GPa, a c = 3.369mm, Q p = 0.96 and Q v = 1.001. Fig.14shows that the design process is able to satisfy

Figure 14 :

 14 Figure 14: Performance of the contact interface along the design procedure. (a) Evolution of maximum contact pressure and von Mises stress. (b) Evolution of the contact area. (c) Evolution of the contact pressure and von Mises stress qualities. (d) Evolution (from blue to red lines) of the contact pressure along the contact line for the design case. l refers to arc length and the red shape refers to the last iteration.

Table 3 :

 3 Evaluation of the pressure distribution quality present in other works based on the reported graphs. We used a plot digitiser tool to compute Q p .

	Reference	Q p Strategy
	Contact
	profile as
	design
	variable.

Acknowledgements

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


Pablo Arrighi
email: parrighi@imag.fr.

Gilles Dowek
email: gilles.dowek@inria.fr.

The principle of a finite density of information

The possibility to describe the laws of the Universe in a computational way seems to be correlated to a principle that the density of information is bounded. This principle, that is dual to that of a finite velocity of information, has already been investigated in Physics, and is correlated to the old idea that there is no way to know a magnitude with an infinite precision. It takes different forms in classical Physics and in quantum Physics.

Why cellular automata?

Stephen Wolfram has advocated in [Wol02] the idea that cellular automata might be a relevant formalism to study the laws of the Universe. This thesis can be seen as a consequence of three more fundamental principles that are implicit in the definition of the notion of a cellular automaton.

To define a cellular automaton, we must discretize space by partitioning it into an infinite set of identical cells, for instance cubes. We must also discretize time and observe the Universe at regular time steps. Then come three assumptions 1. that the state space of each cell is finite, 2. that the state space of a cell at a given time step depends only on the state of a finite number of neighbours at the previous time step,

3. that the state space is the same for all cells and the evolution acts the same everywhere and everywhen.

Assumption (3.) is a reformulation of a well-known principle of Physics: the homogeneity of space and time. So is assumption (2.), which is a reformulation of the bounded velocity of information. Assumption (1.), in contrast, seems to express the new idea, that the density of information also is bounded. Such a principle can be formulated as the fact that the amount of information that can be stored in a bounded region of space is bounded. It can also be formulated as the fact that the cardinality of the state space of a bounded region of space is bounded, since the amount of information storage is the logarithm of this cardinality.

There are other assumptions in the definition of a cellular automaton, such as the fact that space-time is absolute and passive: the neighbourhood relation on cells does not evolve as their states do. Weakening these assumptions is possible, leading for instance to causal networks [START_REF] Wolfram | A new kind of science[END_REF] or to causal graph dynamics [START_REF] Arrighi | Causal graph dynamics[END_REF].

There is a historical discrepancy between the idea of a finite velocity of information and that of a finite density of information. The first has had a clear status as a principle of Physics, since special relativity, and we even know the bound on the velocity. The second seems to be less established. Yet, it is not completely new, as all three principles (1.), (2.), and (3.), have been stated by Robin Gandy in [START_REF] Gandy | Church's thesis and principles for mechanisms[END_REF], and the principle of a bounded density of information has also been stated by Jacob Bekenstein-although the Bekenstein bound is not on the amount of information but on the quotient of the amount of information and the energy contained in a sphere [START_REF] Bekenstein | Universal upper bound to entropy-to-energy ratio for bounded systems[END_REF] .

Two dual principles

The principles of a bounded velocity and of a bounded density of information play a dual role.

For instance, the amount I of information that can be output by a wire of section S in a given time t is bounded by the amount of information contained in the cylinder of section S and length l, where l is the distance that information can travel in time t.

If information has a bounded velocity, then l is bounded and if moreover information has a finite density, the amount of information I is bounded. Thus we need both principles to prove that the information flow of a wire is

l G G G G G GG G G G G G G G G G G G G G G G G G G
bounded. If one of these principles were invalid, it would be easy to build a wire with unbounded flow.

A more evolved example is the impossibility to build a infinitely parallel computing machine, i.e. a machine formed with an infinite number of processors communicating with a common server. If information has a bounded density, then the size of each machine has a lower bound and some machines must be arbitrarily far from the server and if moreover information has a finite velocity, some machines are too far to communicate with the server in due time. Again, if one of these principles were invalidated it would be easy to build an infinitely parallel machine, either with processors of the same size communicating infinitely fast with the server, or with arbitrarily small ... ... machines, of size 1, 1/2, 1/4, 1/8, 1/16, . . . that would lie at a finite distance of the server.

... It seems that these two principles could even be summarized in one, as the fact that in a given time, a piece of information can only travel in a space which is itself populated by a finite amount of information.

Finally, the two principles are comparable, and both differ from the third: homogeneity of space and time. Indeed, the latter is pre-supposed by the very notion of a 'Physics law'. For instance, suppose that some law would be valid

G G G G G G G G G G
in Paris but not in Boston, then this would not be much of a law. The word 'law' itself carries a notion of universality.

But what about the other two? It turns out that a it is possible to do very respectable Physics without these principles. Yet, it is still the case a notion of applicability of a Physics laws underpins both the finite velocity of information and the bounded density of information principles. If a law broke the finite velocity of information principle, this would entail that, in order to determine the new state of a physical system after one second has past, one may potentially need to look at an unbounded number of neighbours. This is no doubt unpractical, unless far-away contributors can be neglected.

In the same way, if a law broke the finite density of information principle, in order to determine the new state of a physical system, one would need to use a description of the state of the neighbours that might potentially contain an unbounded amount of information. Again this law would hardly be applicable, unless this information can be rounded up in some way.

Newtonian Physics

It is well-know that Newtonian Physics contradicts the principle of a finite velocity of information. Whenever a bee flies half a centimeter, it instantaneouly modifies the gravitational field on a distant galaxy. It is only with General Relativity that we understood how to control the consequences of the flight of a bee.

Newtonian Physics not only contradicts the finite velocity principle, it also contradicts the principle of a finite density of information as the distance between the jaws of a caliper for instance is measured by a real number. So if we imagine an infinite amount of information expressed, for instance, as an infinite sequence of 0 and 1, we can see this infinite sequence as the digits of a real number and place the jaws of a caliper at this exact distance, recording this way the infinite amount of information in a finite volume.

Of course, we have all been taught that this idea is nonsense and that the distance between the jaws of a caliper is defined with a finite precision: three or four digits. It is surprising that this thesis that a physical magnitude is always defined with a finite precision has never been given the status of a Physics principle as have the homogeneity of space and time. Yet, this thesis may be a very early occurrence of this idea of a finite density of information.

The finite-density versus the superposition principle

The principle of a finite density of information seems to be challenged by quantum theory and its superposition principle. Indeed, regardless of whether a region of space is bounded or not, whenever its state space contains two states u and v, then it must also contain all the linear combinations λu + µv up to a renormalization factor. Hence it is infinite, and the finite density principle cannot be formulated as the fact that the set of states of a bounded region of space is finite. These complex amplitudes λ and µ can be compared with probabilities. But this comparison has its limits. Indeed, since probabilities are positive real numbers, they can only add up constructively to make and event more likely. In contrast amplitudes are complex numbers, then can also be substracted to one another, i.e. add up destructively, making an event either more likely of less likely; or sometimes even impossible. This is why when two different scenarios happen in a superposition, we cannot quite think of them as two non-interacting branches of what may happen in the Universe: these amplitudes cannot be ignored. This makes the superposition principle is a difficulty when extending the finite-density principle to quantum theory.

Yet, the infinity of the state space does not mean that the amount of possible outcomes, when measuring the state of this region, is itself infinite, because in quantum theory, we cannot measure the state of the system, but only one of its observables. Thus, an alternative formulation of the bounded density of information principle is that each projective measurement of a bounded region, at any given point in time, may only yield a finite number of possible outcomes. This requirement amounts to the fact that the state space of a bounded region of space is a finite-dimensional vector space. This constitutes a good quantum alternative to the formulation of the finite density of information principle, one which does not demand anything to be actually measured in any way.

The finite-density principle versus correlations

Moreover, yet another problem arises in the quantum setting, mamely that of 'correlations'; which in this context are also referred to as 'entanglement'. To understand what entanglement is, one must apply the superposition principle again, but this time to pairs of systems. For instance if systems A and B may be in state uu, meaning that both are in state u, or in state vv, meaning that both are in state v, then the state λ(uu) + µ(vv) is also allowed. This entangled state corresponds to a superposition of two global, and everywhere different scenarios. It is possible to get an intuition of the meaning of this entangled state by appealing to our understanding about usual probabilities. In probability theory, it is a commonplace that knowing the marginal probabilities over the individual systems does not entail knowing the joint probability distributions. For instance Alice and Bob may have a half-half chance of having a blue or red ball in each of their boxes, but it could be that each of their boxes have been prepared in such a way that both balls are of the same colour. The state λ(uu) + µ(vv) tells the same story not for probabilities, but for amplitudes. This is an issue, because the bounded velocity of information principle seems very weak in that respect: after all, all it says is that the state of each individual system is determined by that of the neighbours, but it does not say how the entanglement between the individual states is determined. If we retake the probability analogy; the bounded velocity of information principle tells you that the probabilities of Alice finding a red or a blue ball in her box are determined locally; but what about the probability that this ball is of the same color as Bob's?

Fortunately on this issue of whether the correlations can be determined locally, quantum theory is a bit friendlier than probability theory. Indeed, it has been shown that whereas the bounded velocity of information principle is way too weak in the context of probabilistic cellular automata [START_REF] Arrighi | Applying causality principles to the axiomatization of Probabilistic Cellular Automata[END_REF], surprisingly it turns out to be sufficient in the context of quantum cellular automata [START_REF] Arrighi | Unitarity plus causality implies localizability (full version)[END_REF] for the evolution to be described as a composition of local functions mapping neighbourhoods to neighbourhoods, that is functions from a finite-dimensional vector space to itself. Another way to phrase this is that in the quantum setting, any evolution that respects the bounded velocity of information principle necessarily takes the form of a circuit of gates being applied on systems and their neighbours. Thus correlations can be 'tamed' and the global evolution broken into functions from finite-dimensional vector spaces to themselves.

A computable quantum Universe

Yet, moving from a finite set to a vector space, even a finite-dimensional one, is a big leap and many advantages of describing the Universe as a cellular automaton might be lost by doing so. Indeed, the main advantage of having a finite state space S for cells was that the local evolution function of the cellular automaton, that is a function from S n to S, where n is the number of neighbours of a cell, was a function from a finite set to a finite set. Such a function can be described by a table giving the image of each element of the domain and it is always computable. This point was the key to the computability of cellular automata, and thus to the idea that the Universe is computable if it can be described as a cellular automaton. This idea of a computable Universe-the physical Church-Turing thesis-is one of the main goals of both Robin Gandy and Stephen Wolfram.

Fortunately, another feature of the quantum theory comes to repair the damages created by the superposition principle. According to the quantum theory, theses local functions of a quantum cellular automata cannot be any function from neighbourhoods to neighbourhoods, but must be linear. There are much less linear functions than functions from a finite-dimensional vector space to another. In particular any such linear function can be described by a matrix, and matrix multiplication boils down to additions and multiplications that are computable operations. Thus if scalars are chosen in an adequate manner as discussed below, all such functions are computable. Thus, infinity is tamed by finite-dimension and linearity and this formulation of the principle of a finite density of information is sufficient to prove the Universe computable.

Scalars

When saying that addition and multiplication of scalars are computable, we need to be a little careful, as scalars are complex numbers and computability is mostly defined for natural numbers.

Computability can be extended from natural numbers to real and complex numbers [START_REF] Weihrauch | Computable analysis: an introduction[END_REF], but then the picture is a somewhat different. Not only functions from complex numbers to complex numbers can be computable or not, but complex numbers themselves can be computable or not. Thus, although multiplication is computable, multiplication by a fixed scalar may be non-computable if the fixed scalar is non-computable.

Michael Nielsen has noticed that having a process building a superposed state λu + µv, where λ and µ are non computable numbers, can per se lead to non computability, as this state can be tomographed by repeated measurements and λ be extracted with increased precision from the state.

Thus, a physical system has two ways to encode an infinite amount of information: either by having an infinite dimensional state space, or by encoding an infinite amount of information in a scalar, exactly in the same way as an infinite amount of information can be encoded in a distance in Newtonian Physics.

To express the finite-density principle, it seems that it not sufficient to restrict to finite-dimensional vector spaces for the state spaces of cells, but we must also restrict to a smaller set of scalars, yielding an intersting problem for physicists about the fields that can be used for the scalars in quantum theory.

With this precise formulation of the finite density of information, we can formally prove that the physical Church-Turing thesis is a consequence of the three hypotheses we have mentioned: homogeneity of space and time, bounded velocity of information and bounded density of information [START_REF] Arrighi | The physical Church-Turing thesis and the principles of quantum theory[END_REF].

Some scientists, such as David Deutsch [START_REF] Deutsch | Quantum theory, the Church-Turing principle and the universal quantum computer[END_REF] propose to take the physical Church-Turing thesis as a principle of Physics. Some others, such as Michael Nielsen [START_REF] Nielsen | Computable functions, quantum measurements, and quantum dynamics[END_REF], propose to deduce this thesis from more fundamental principles of Physics. Like Robin Gandy, we have chosen to start from these three fundamental principles. An advantage of this approach is that these three principles are falsifiable provided explicit bounds are given. Thus, from a logical point of view, they are purely universal formulas, while the physical Church-Turing thesis itself has a much more complex logical form ∀∃∀, namely for all physical system m, there exists a computable function f such that forall initial state x, the evolution of m leads from x to f (x).

Conclusion: cellular automata and quantum cellular automata

It is one thing to say that the Universe is computable and thus can be simulated by any model of computation, e.g Cellular Automata. It is another to seek to closely model physical systems in a space-preserving manner. Pursuing this second, more demanding aim leads us to investigate how to model the state space of a bounded region of space. Supposing that this state space is finite is adequate in the classical setting, but must be refined in the quantum setting to a finite-dimensional vector space over a reasonable field. Indeed it seems that simulating a quantum systems by classical cellular automata not only leads to an exponential slowdown, but also fails to respect the geometry of the system. This justifies Feynman's proposal to model quantum systems with Quantum Cellular Automata. This quantum extension of the Cellular Automata model shows its liveliness and robustness.


Anson Ma

Francisco Chinestat

Malcolm Mackley

The Rheology and Microstructure of Carbon N anotube ( CNT) Suspensions

Keywords: Carbon Nanotube suspension, Processing, Rheo-optics, Shear-induced assembly PACS: 6L46,Fg, 6L48,+c, 8Ll6,qVc, 8Ll6,Fg

This conference paper is concerned with the rheology and microstructure of Carbon Nanotube suspensions, Carbon Nanotubes (CNTs) are cylinders of rolled graphene sheets of high aspect ratio and they belong to a relatively new class of fibrous material that can potentially be used for high-performance nano-composites and nano-devices, For these applications, understanding the rheology and microstructure of CNT suspensions is of key importance as this would allow for processing CNTs into macroscopically usable form with controlled properties, This paper reports the rheology and optical microstructure of surface-treated and untreated CNT suspensions, In terms of rheology, both types of CNT suspensions exhibited a shear-thinning characteristic in steady simple shear flows and their rheological responses were subsequently modelled using a Fokker-Planck (FP) based orientation model that incorporates aggregation effect in the case of untreated CNT suspension, In terms of unusual mesostructures, it was observed that by applying appropriate flow conditions to an untreated CNT suspension, a novel type of mesoscopic filament structure could be formed, where the CNT filaments aligned perpendicular to the flow direction and in the vorticity direction, The formation of these anisotropic structures has similarities to the hand-spinning of cotton into fibres and these flow-induced structures might have strategic benefits for some potential applications of CNT

INTRODUCTION

Carbon Nanotubes (CNTs) are rolled cylinders of graphene sheets having high aspect ratio, low density, high mechanical strength and high thermal conductivity [START_REF] Iijima | Helical microtubules of graphitic carbon[END_REF] [2], Because of these intrinsic properties, they have been researched in a number of application areas such as high-performance nanocomposites and nanodevices [START_REF] Calvert | Nanotube composites: A recipe for strength[END_REF] [START_REF] De Heer | A Carbon Nanotube field-emission electron source[END_REF], For most of these applications, understanding the rheology and microstructure of CNT suspensions is important in relation to processing this type of nanomaterials into some macroscopically usable forms with predictable properties, The objectives of this paper are to highlight several key experimental findings, both in terms of rheological and optical studies of CNT suspensions, and secondly, to discuss briefly the basis on which CNT suspensions can be modelled using a Fokker-Planck (FP) based diffusion-advection modeL Finally, flow-induced formation of some unusual CNT mesostructures is also described,

CARBON NANOTUBE SUSPENSIONS

At present, enhancing the bulk conductivity of a polymer film is identified as one of the most promising potential applications for CNTs [START_REF] Sandler | Development of a dispersion process for Carbon Nanotubes in an epoxy matrix and the resulting electrical properties[END_REF] [START_REF] Shaffer | Dispersion and packing of Carbon Nanotubes[END_REF] and for such an application, it is thought desirable to have a well-dispersed system of CNT. However, this can be a laborious and challenging task because in general, a mixing time of several hours is required in order to establish a stable rheology for CNT suspensions, and even then dispersion at the nano-scale is not guaranteed [START_REF] Huang | Dispersion rheology of Carbon Nanotubes in a polymer matrix[END_REF]. In this regard, introducing functional groups onto the surface of CNTs, also known as functionalization, has proven to be one possible way of enhancing the dispersion of CNT in a processing matrix [START_REF] Banerjee | Rational chemical strategies for Carbon Nanotube functionalization[END_REF]. The actual chemistry involved in surface-treating CNTs may differ, but functionalization generally works on the principle of increasing steric hindrance or introducing electrostatic repulsion between CNTs, thereby preventing effective 7r -7r stacking and subsequent aggregation of CNTs (see for example [START_REF] Dyke | Unbundled and highly functionalized Carbon Nanotubes from aqueous reactions[END_REF]). Figure 1 shows typical optical micrographs of surface-treated and untreated CNTs suspended in epoxy resin1 [START_REF] Ma | Rheological modelling of suspensions involving surfacetreated Carbon Nanotubes[END_REF]. It is clear from the micrographs that the untreated CNT suspension exhibited an optical microstructure (Fig. l(a)), whilst by contrast, there is no discernable microstructure for the treated CNT suspension (Fig. l(b)), indicating that the mixture was well-dispersed at the micron-level.

RHEOLOGY

The rheology of both treated and untreated CNTs suspended in epoxy resin were studied and modelled for simple rheometric flows. It was observed experimentally that both types of suspensions exhibited a shear-thinning characteristic when subject to steady simple shear and the apparent viscosity ( 1Ja) of the suspension decreased asymptotically to the matrix viscosity as shear rate increased (Fig. 2). Given the same weight concentration of CNTs, the untreated CNT suspension exhibited a significantly larger viscosity enhancement effect at low shear rates compared with the treated CNT suspension, as also shown in Fig. 2. . .. For treated CNT suspensions, the observed shear-thinning characteristic can be modelled using a simple orientation model that in the past has been used for short-fibre suspensions [START_REF] Petrie | The rheology of fibre suspensions[END_REF]. The model basically assumes that CNTs behave essentially as rigid rods and that they can orient and align in a simple shear flow. There are, however, also randomizing events due to Brownian motion and tube interaction (for semi-dilute or concentrated suspensions), which can be represented using a rotary diffusion coefficient (Dr) [START_REF] Folgar | Orientation behaviour offibers in concentrated suspensions[END_REF] and the time evolution of CNT orientation distribution can be described using a Fokker-Planck based diffusion-advection equation [START_REF] Hinch | The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles[END_REF]. Dr and an adjustable parameter NP that depends on the concentration and aspect ratio are the key parameters in the model and such model can successfully describe the shear-thinning characteristics of treated CNT suspensions as shown in Fig. 3. Besides simple shear flows, the model was also capable of describing some kinetic responses of CNT suspensions in smallamplitude oscillatory shear flows, thereby supporting the basis for the model itself [START_REF] Ma | Rheological modelling of suspensions involving surfacetreated Carbon Nanotubes[END_REF].

• • • • • • • • • • • • • • = = • • • • ~ 0 0 0 V V ~ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ~ ~ $ $ ""' (
However, in the case of untreated CNT suspensions where aggregates of CNTs were optically observed, the simple orientation model failed to describe the shear-thinning characteristic observed experimentally, and therefore a more sophisticated model is needed. An Aggregation/Orientation (AO) model has been proposed in this context by considering extra contribution from CNT aggregates [START_REF] Ma | Rheological modelling of suspensions involving untreated Carbon Nanotubes in shear flow[END_REF]. The model assumes that in an aggregating suspension of CNTs, CNTs can exist in different entanglement states and they can have different mobility to orient in a simple shear flow. The Fokker-Planck equation was then modified to accommodate for a hierarchy of microstructures. In this model, transformation between microstructures is possible depending on the shear rate applied to the system, and thus matching the phenomena of CNT aggregation and disaggregation as observed experimentally [15] [16].

UNUSUAL FLOW MICROSTRUCTURE

The flow microstructure of untreated CNT suspensions was followed using the Cambridge Shear System (Linkam Instruments Inc.), which essentially consists of a shear cell coupled with an optical microscope. Although shear can induce alignment of CNTs, it was also reported that low shear could induce aggregation of CNTs [15] [16]. In most cases, a homogeneous dispersion of CNT is desired, however there are other situations where a highly heterogeneous structure can be of strategic benefit. It was discovered that some unusual filament mesostructures could be formed by applying appropriate flow conditions to an untreated CNT suspension [START_REF] Ma | Experimental observation on the flow-induced assembly of Carbon Nanotube suspensions to form Mesoscopic Helical Filaments[END_REF]. Figure 4 shows the time evolution of optical texture as low shear was applied to some untreated CNTs (0.03%) suspended in a low-viscosity epoxy matrix2 • Initially, there were only isotropic aggregates of CNTs (Fig. 4(a)), but as they rotated within a simple shear flow, they drew in nearby CNT aggregates to form long Mesoscopic Helical Filaments (MHFs) that aligned perpendicular to the flow and in the vorticity direction. The way in which CNT aggregates are helically wound to form this kind of filament structure is intriguing

and not yet fully understood. In terms of applications, although the advantage of having a MHF structure of CNTs has not yet been explored, it is expected that the anisotropic nature of these mesostructures will produce preferential conductivity along the length of the filaments and thereby providing unusual anisotropic electrical properties to the system as a whole. 1 ) was applied to 0.03% untreated CNTs suspended in a low-viscosity epoxy matrix (gap size= 130 J-Lm) [START_REF] Ma | Experimental observation on the flow-induced assembly of Carbon Nanotube suspensions to form Mesoscopic Helical Filaments[END_REF].

CONCLUSIONS

In conclusion, this paper reported recent findings on the rheo-optical studies of Carbon Nanotube (CNT) suspensions which are of general relevance to the processing of CNTs.

It was found that for surface-treated CNT suspensions that exhibit no discernable optical microstructure, CNTs can essentially be modelled as rigid fibres that can orient in a shear flow. The rheology of flowing CNTs can be successfully described using a simple orientation model and essential information on the orientation distribution of CNTs can be obtained from a Fokker-Planck (FP) description. The FP based orientation model was capable of predicting both steady and transient responses of treated CNT suspensions, but failed in the case of untreated CNT suspensions where CNT aggregates were optically observed. This has led to the development of a more appropriate Aggregation/ Orientation (AO) model, in which a hierarchy of CNT aggregate structures was considered. Finally, in terms of optical microstructure, although low shear can induce aggregation and should be avoided in cases where a uniformly dispersed system of CNTs is desired, it can also be used strategically to produce some unusual Mesoscopic Helical Filament (MHF) structures that are highly anisotropic.

FIGURE 1 .

 1 FIGURE 1. Optical micrographs of 0.3% (a) surface-treated and (b) untreated CNT suspensions [1 0].

  0.3% Untreated CNT . • ... / . . . . . . . .. . . . . . . . .

FIGURE 2 .

 2 FIGURE 2. Apparent viscosity (1Ja) as a function of shear rate (i) for (a) surface-treated and (b) untreated CNT suspensions [1 0].

FIGURE 3 .

 3 FIGURE 3. Simple orientation model fitting of 0.2% and 0.3% treated CNT suspensions. Dr and Np are the key fitting parameters in the model. (1)a =apparent viscosity and y =shear rate.) [l 0]

FIGURE 4 .

 4 FIGURE 4. Time evolution of optical texture as low shear (y = 0.5 s-1 ) was applied to 0.03% untreated CNTs suspended in a low-viscosity epoxy matrix (gap size= 130 J-Lm)[START_REF] Ma | Experimental observation on the flow-induced assembly of Carbon Nanotube suspensions to form Mesoscopic Helical Filaments[END_REF].

The surface-treated CNTs were single-walled and were supplied by Nanocomposites Inc., whereas the untreated CNTs were multi-walled and were provided by the Department of Materials Science and Metallurgy, University of Cambridge. The epoxy resin used was Araldite LY556 from Huntsman Inc.

The low-viscosity epoxy has a viscosity of0.4 Pa.s and was supplied by Epoxies Inc. (UV 60-7155).

ACKNOWLEDGMENTS

We would like to thank Nanocomposites Inc. and Prof. A. H. Windle from the Department of Materials Science and Metallurgy at the University of Cambridge for providing the base Carbon Nanotubes. Anson Ma would also like to to acknowledge the Croucher Foundation Scholarship (Hong Kong) and the Overseas Research Students Awards Scheme (ORSAS) (United Kingdom) for providing financial support.