# On the effective conductivity and the apparent viscosity of a thin rough polymer interface using PGD-based separated representations

Amine Ammar, Chady Ghnatios, Frank Delplace, Anaïs Barasinski, Jean-louis

Duval, Elías Cueto, Francisco Chinesta

► To cite this version:

Amine Ammar, Chady Ghnatios, Frank Delplace, Anaïs Barasinski, Jean-louis Duval, et al.. On the effective conductivity and the apparent viscosity of a thin rough polymer interface using PGD-based separated representations. International Journal for Numerical Methods in Engineering, 2020, 121, pp.5256-5274. 10.1002/nme.6448. hal-02902716

HAL Id: hal-02902716

https://univ-pau.hal.science/hal-02902716

Submitted on 26 Aug 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the effective conductivity and the apparent viscosity of a thin rough polymer interface using PGD-based separated representations

Amine Ammar<fnref n="1" /> | Chady Ghnatios<fnref n="2" /> | Frank Delplace<fnref n="3" /> | Anais Barasinski<fnref n="4" /> | Jean-Louis Duval<fnref n="3" /> | Elías Cueto<fnref n="5" /> | Francisco Chinesta<fnref n="6" />

**Summary** Composite manufacturing processes usually proceed from pre-impregnated preforms that are consolidated by simultaneously applying heat and pressure, so as to ensure a perfect contact compulsory for making molecular diffusion possible. However, in practice, the contact is rarely perfect. This results in a rough interface where air could remain entrapped, thus affecting the effective thermal conductivity. Moreover, the interfacial melted polymer is squeezed flowing in the rough gap created by the fibers located on the prepreg surfaces. Because of the typical dimensions of a com-posite prepeg, with thickness orders of magnitude smaller than its other in-plane dimensions, and its surface roughness having a characteristic size orders of mag-nitude smaller than the prepreg thickness, high-fidelity numerical simulations for elucidating the impact of surface and interface roughness remain today, despite the impressive advances in computational availabilities, unattainable. This work aims at elucidating roughness impact on heat conduction and the effective viscosity of the interfacial polymer squeeze flow by using an advanced numerical strategy able to reach resolutions never attained until now, a sort of numerical microscope able to attain the scale of the smallest geometrical detail.

**KEYWORDS:**

Surface and interface roughness; Thermal conductivity; Effective viscosity, Numerical simulation; Model Order Reduction; Space separated representations; PGD.

# 1 | INTRODUCTION

Many composite forming processes for elaborating structural parts are based on the consolidation of preimpregnated preforms in form of sheets or tapes, that are put in contact, heated and compressed to ensure its bonding.

However, the welding of two thermoplastic layers requires specific physical conditions: an almost perfect contact, also called intimate contact, and a temperature that has to be high enough during a time large enough to ensure the diffusion of macromolecules across the interface enabling the bonding, while limiting as such as possible thermal-induced material degradation.

The efficient numerical simulation of automated tape placement, including all the just referred mechanisms, was addressed by the authors in<fnref n="1" /> (the interested reader can also refer to the references therein concerning the process modeling) where the use

<footnote n="1">¹Arts et Métiers ParisTech - Angers, 2 boulevard du Ronceray, BP 93525, 49035 Angers cedex 01, France</footnote>

<footnote n="2">²Notre Dame University - Louaize, P.O. Box 72, Zouk Mikael, Zouk Mosbeh, Lebanon</footnote>

<footnote n="3">³ESI Group, Bâtiment Seville, 3 bis rue Saarinen, F-50468, Rungis, France</footnote>

<footnote n="4">⁴Composites Chair E2S UPPA Arkema-Canoe, IPREM, Helioparc, 2 Av. Pdt. Angot, 64000 Pau, France</footnote>

<footnote n="5">⁵Aragon Institute of Engineering Research, Universidad de Zaragoza. Maria de Luna, s.n. Zaragoza, Spain</footnote>

<footnote n="6">⁶ESI Chair @ Arts et Métiers ParisTech, PIMM Lab, 151 Boulevard de l'Hôpital, 75013 Paris, France</footnote>
of advanced numerical techniques (revisited and discussed later) was emphasized. The same numerical techniques were then considered for describing the compression of pre-impregnated laminates within a fully resolved 3D approach.²,³

As just mentioned, consolidation implies putting plies in contact while supplying heat and pressure. The former promotes molecular diffusion at the plies interface and both (heat and pressure) facilitate the intimate contact by squeezing surface asperities. Interfacial thermal properties of two rough surfaces in contact are of major relevance for evaluating the process performances and manufactured part properties. On the other hand evaluating the interfacial polymer flow in the rough interface gap allows evaluating the asperities squeeze inducing intimate contact.

Two widely considered descriptions of surfaces are based on the use of a fractal⁴ whose parameters can be extracted from the spectral power density of the surface profile,⁵ or by assuming that the surface is composed by a set of rectangular (parallelepipeds in 3D) elements representing the so-called surface asperities.⁶ More elaborated representations consists of approximating the real surface by using an adequate functional basis, being wavelets specially adapted because its inherent multi-resolution (multi-scale) properties.⁷

Since the main phenomenon which drives the adhesion of surfaces in contact is molecular diffusion, intimate contact becomes of major relevance. Many studies in the past considered the deformation of loaded surfaces, from the seminal work of Hertz⁸ to more recent studies ⁹,¹⁰,¹¹,¹²,¹³, where distribution of interstices, i.e., the difference of heights between both interfaces, the shape of the surface tips, etc. were taken into consideration.

High-resolution simulation requires specific discretization techniques, as the ones described later. By using them the analysis of thermal properties of fractal surfaces was addressed in⁵ and the subsequent consolidation (squeezing flow) occurring as response to the applied pressure was considered in¹⁴. Later in⁷ real surfaces were modeled and simulated, being described by using a wavelet-based multi-resolution analysis. However, these studies considered the fractal surface in contact with a planar one and consequently the conclusion were too dependent on that modeling hypotheses as was proved in¹⁵.

On the other hand, flows in narrow gaps were successfully simulated using advanced discretization techniques based on the use of in-plane-out-of-plane separated representations ²,³ allowing extremely fine representations in such degenerated domains (with the characteristic thickness dimension orders of magnitude smaller than the in-plane dimensions), situations where mesh-based discretizations fail to obtain a valid solution of the flow problem.

In practice contact is rarely perfect, a fact that results in a roughness interface where air could remain entrapped. Because of the typical dimensions of a composite prepeg, with its thickness orders of magnitude smaller than its other in-plane dimensions, and its surface roughness having a characteristic size orders of magnitude smaller than the prepreg thickness, high-fidelity numerical simulations for elucidating the impact of surface and interface roughness remains today, despite the impressive advances in computational availabilities, unreachable. This work aims at elucidating roughness impact on heat conduction and the effective viscosity of the interfacial polymer squeeze flow by using an advanced numerical strategy able to reach resolutions never attained until now as discussed later, a sort of numerical microscope able to attain the scale of the smallest geometrical detail.

In order to explore two kinds of configuration, the first consisting of harmonic interfaces (looking for elucidating the role of the gap amplitude, frequency and phase) and the second with the topography described from a random motion (fractal description) we will consider the former in the heat transfer analysis whereas the last will be considered when addressing the flow problem in the thin gap representing the thick interface.

After this short introduction, next section revisits the in-plane-out-of-plane representations for simulating heat transfer and fluid flows in narrow gaps. Then, Section 3 will focus on the separability issues as well as on the roughness representation. Section 4 formulates the problem in the reference and fully separable domain. Section 5 concerns the thermal problem and its associated effective thermal conductivity whereas Section 6 focusses on the flow problem and the associated apparent viscosity (or effective thickness gap). Finally Section 7 addresses some major conclusions.

## 2 | REVISITING HEAT TRANSFER AND FLUID FLOWS IN NARROW GAPS USING PGD-BASED SEPARATED REPRESENTATIONS

The main issue of solving problems in narrow gaps is the domain degeneracy, being one dimension (here the one related to the thickness) much smaller than the others. In this situation standard mesh-based discretization techniques fail, because the resulting meshes involve too many degrees of freedom.

Space-separated representations within the Proper Generalized Decomposition (PGD) framework seems a valuable route for alleviating such difficulties as discussed below. The interested reader can refer to the abundant references on the PGD.¹⁶
These separated representations were successfully applied for addressing the solution of problems defined in plate¹⁷ and shell¹⁸ geometries. The same approach was extensively considered in structural plate and shell models.¹⁹,²⁰

## 2.1 | Heat transfer in a multilayered plate

In what follows we are illustrating the construction of the Proper Generalized Decomposition of a generic model defined in a plate domain $\Omega = \Xi \times I$ with $\Xi \subset \mathbb{R}^2$ and $I = [0, H] \subset \mathbb{R}$. For the sake of simplicity we consider the model related to the steady state heat conduction equation

$$ \nabla \cdot (\mathbf{K}\nabla u) = 0, \qquad (1) $$

in a plate geometry that contains $P$ plies along the plate thickness. Each ply is characterized by its conductivity tensor $\mathbf{K}_i(x, y)$ which is assumed constant through the ply thickness. Moreover, without any loss of generality, we assume the same thickness $h$ for the different plies constituting the laminate. Thus, we can define a characteristic function representing the position of each ply $i = 1, \dots, P$

$$ \chi_i(z) = \begin{cases} 1 & z_i \le z \le z_{i+1}, \\ 0 & \text{otherwise} \end{cases}, \qquad (2) $$

where $z_i = (i - 1)h$ defines the location of the $i$-th ply in the laminate thickness. The laminate conductivity can be given in the following separated form

$$ \mathbf{K}(x, y, z) = \sum_{i=1}^{P} \mathbf{K}_i(\mathbf{x}) \cdot \chi_i(z), \qquad (3) $$

where $\mathbf{x}$ denotes the in-plane coordinates, i.e., $\mathbf{x} = (x, y) \in \Xi$.

The weak form of Eq. (1), with appropriate boundary conditions (that will be specified later) reads

$$ \int_{\Omega} \nabla u^* \cdot (\mathbf{K} \nabla u) \, d\Xi = 0, \qquad (4) $$

with the test function $u^*$ defined in an appropriate functional space. The solution $u(x, y, z)$ is searched under the so-called “in-plane-out-of-plane” separated form,

$$ u(\mathbf{x}, z) \approx \sum_{j=1}^{N} X_j(\mathbf{x}) \cdot Z_j(z). \qquad (5) $$

In what follows we are illustrating the construction of such a decomposition. For this purpose we assume that at enrichment step $n < N$ the solution $u^{n-1}(\mathbf{x}, z)$ is already known,

$$ u^{n-1}(\mathbf{x}, z) = \sum_{j=1}^{n-1} X_j(\mathbf{x}) \cdot Z_j(z), \qquad (6) $$

and that at the present step we look for the solution enrichment $R(\mathbf{x}) \cdot S(z)$:

$$ u^n(\mathbf{x}, z) = u^{n-1}(\mathbf{x}, z) + R(\mathbf{x}) \cdot S(z). \qquad (7) $$

The test function involved in the weak form is searched under the form

$$ u^*(\mathbf{x}, z) = R^*(\mathbf{x}) \cdot S(z) + R(\mathbf{x}) \cdot S^*(z). \qquad (8) $$

By introducing Eqs. (7) and (8) into Eq. (4) it results

$$ \int_{\Omega} \left( \left( \begin{pmatrix} \tilde{\nabla} R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} + \begin{pmatrix} \tilde{\nabla} R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \right) \cdot \begin{pmatrix} \mathbf{K} \\ R \cdot \frac{dS}{dz} \end{pmatrix} \right) d\Xi = - \int_{\Omega} \left( \left( \begin{pmatrix} \tilde{\nabla} R^* \cdot S \\ R^* \cdot \frac{dS}{dz} \end{pmatrix} + \begin{pmatrix} \tilde{\nabla} R \cdot S^* \\ R \cdot \frac{dS^*}{dz} \end{pmatrix} \right) \cdot \mathbf{Q}^n d\Xi, \qquad (9) $$

where $\tilde{\nabla}$ denotes the plane component of the gradient operator, i.e., $\tilde{\nabla}^T = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$ and $\mathbf{Q}^{n-1}$ denotes the flux

$$ \mathbf{Q}^{n-1} = \mathbf{K} \sum_{j=1}^{n-1} \begin{pmatrix} \tilde{\nabla} X_j(\mathbf{x}) \cdot Z_j(z) \\ X_j(\mathbf{x}) \cdot \frac{dZ_j(z)}{dz} \end{pmatrix}. \qquad (10) $$

Now, as the enrichment process is non-linear we propose to search the couple of functions $R(\mathbf{x})$ and $S(z)$ by applying the alternating direction fixed point algorithm. Thus, assuming $R(\mathbf{x})$ known, we compute $S(z)$, and then we update $R(\mathbf{x})$. The process continues until reaching convergence. The converged solutions allow defining the next term in the finite sum decomposition:
$R(\mathbf{x}) \to X_n(\mathbf{x})$ and $S(z) \to Z_n(z)$. The global enrichment procedure stops when the norm of the equation residual given y Eq. (1) becomes small enough, even if more advanced stopping criteria exist.¹⁶

## 2.2 | Squeeze flow in composite laminates

The in-plane-out-of-plane separated representation allows the solution of 3D flow models defined in plate geometries with a computational complexity characteristic of 2D simulations. In the present case the 3D velocity field reads

$$ \mathbf{v}(\mathbf{x}, z) = \begin{pmatrix} u(\mathbf{x}, z) \\ v(\mathbf{x}, z) \\ w(\mathbf{x}, z) \end{pmatrix} \approx \begin{pmatrix} \sum_{i=1}^{N} P_i^1(\mathbf{x}) \cdot T_i^1(z) \\ \sum_{i=1}^{N} P_i^2(\mathbf{x}) \cdot T_i^2(z) \\ \sum_{i=1}^{N} P_i^3(\mathbf{x}) \cdot T_i^3(z) \end{pmatrix}, \quad (11) $$

which leads to a separated representation of the strain rate that, when introduced into the flow problem weak form, allows the calculation of functions $P_i(x, y)$ by solving the corresponding 2D equations and functions $T_i(z)$ by solving the associated 1D equations, as described later.

Eq. (11) can be rewritten, using the Hadamard product, in the compact form

$$ \mathbf{v}(\mathbf{x}, z) \approx \sum_{i=1}^{N} \mathbf{P}_i(\mathbf{x}) \circ \mathbf{T}_i(z). \quad (12) $$

The velocity gradient $\nabla \mathbf{v}(\mathbf{x}, z)$ can be written as

$$ \nabla \mathbf{v} = \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\ \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{pmatrix} \approx \sum_{i=1}^{N} \begin{pmatrix} \frac{\partial P_i^1}{\partial x} & \frac{\partial P_i^1}{\partial y} & P_i^1 \\ \frac{\partial P_i^2}{\partial x} & \frac{\partial P_i^2}{\partial y} & P_i^2 \\ \frac{\partial P_i^3}{\partial x} & \frac{\partial P_i^3}{\partial y} & P_i^3 \end{pmatrix} \circ \begin{pmatrix} T_i^1 & T_i^1 & \frac{\partial T_i^1}{\partial z} \\ T_i^2 & T_i^2 & \frac{\partial T_i^2}{\partial z} \\ T_i^3 & T_i^3 & \frac{\partial T_i^3}{\partial z} \end{pmatrix} = \sum_{i=1}^{N} \mathbb{P}_i(\mathbf{x}) \circ \mathbb{T}_i(z). \quad (13) $$

### 2.2.1 | Stokes model

The Stokes flow model is defined in $\Omega = \Xi \times I$, $\Xi \subset \mathbb{R}^2$ and $I \subset \mathbb{R}$, and for an incompressible fluid, in absence of inertia and mass terms reduces to

$$ \left\{ \begin{array}{l} \nabla \cdot \boldsymbol{\sigma} = \mathbf{0} \\ \boldsymbol{\sigma} = -p\mathbf{I} + 2\eta\mathbf{D}, \\ \nabla \cdot \mathbf{v} = 0 \end{array} \right. , \quad (14) $$

where $\boldsymbol{\sigma}$ is the Cauchy's stress tensor, $\mathbf{I}$ the unit tensor, $\eta$ the fluid viscosity, $p$ the pressure (Lagrange multiplier associated with the incompressibility constraint) and the rate of strain tensor $\mathbf{D}$, symmetric component of the velocity gradient. A penalty formulation is here considered, even if stable mixed formulations compatible with the use of in-plane-out-of-plane separated representations were proposed by the authors.²¹ The penalty formulation can be written as

$$ \nabla \cdot \mathbf{v} + \lambda p = 0, \quad (15) $$

or, more explicitly

$$ p = - \frac{\nabla \cdot \mathbf{v}}{\lambda} = - \frac{\mathrm{Tr}(\mathbf{D})}{\lambda}, \quad (16) $$

where Tr() refers the trace operator.

The weak form of the penalized Stokes problem, for a test velocity $\mathbf{v}^*$ vanishing on the boundary in which the velocity is prescribed, and assuming null traction in the remaining part of the domain boundary, can be written as

$$ \int_{\Omega} \left\{ \frac{1}{\lambda} \mathrm{Tr}(\mathbf{D}^*) \mathrm{Tr}(\mathbf{D}) + 2\eta \mathbf{D}^* : \mathbf{D} \right\} d\mathbf{x} \, d\mathbf{z} = 0, \quad (17) $$
with $\mathbf{D}^*$ is the strain rate related to the test field $\mathbf{v}^*$. The two terms involved in the previous integral form can be expressed in the fully separated form³

$$2\eta \mathbf{D}^* : \mathbf{D} \approx \frac{\eta}{2} \sum_{j=1}^{N} \sum_{k=1}^{4} \left( A_{jk}^*(\mathbf{x}) : \mathbb{B}_{jk}(z) + A_{jk}(\mathbf{x}) : \mathbb{B}_{jk}^*(z) \right), \quad (18)$$

and

$$\frac{1}{\lambda} \mathrm{Tr}(\mathbf{D}^*) \cdot \mathrm{Tr}(\mathbf{D}) \approx \frac{1}{\lambda} \sum_{j=1}^{N} \left( \mathbb{F}_j^*(\mathbf{x}) : \mathbb{G}_j(z) + \mathbb{F}_j(\mathbf{x}) : \mathbb{G}_j^*(z) \right), \quad (19)$$

where the different matrices appearing in the expressions (18) and (19) were detailed in a former work of the authors.³.

The construction of the solution separated representation is performed incrementally, a term of the sum at each iteration. Thus, supposing that at iteration $n-1$, $n \ge 1$, the first $n-1$ terms of the velocity separated representation were already computed

$$\mathbf{v}^{n-1}(\mathbf{x}, z) = \sum_{i=1}^{n-1} \mathbf{P}_i(\mathbf{x}) \circ \mathbf{T}_i(z), \quad (20)$$

the updated velocity field $\mathbf{v}^n(\mathbf{x}, z)$ at iteration $n$ reads

$$\mathbf{v}^n(\mathbf{x}, z) = \sum_{i=1}^{n} \mathbf{P}_i(\mathbf{x}) \circ \mathbf{T}_i(z) = \mathbf{v}^{n-1}(\mathbf{x}, z) + \mathbf{P}_n(\mathbf{x}) \circ \mathbf{T}_n(z), \quad (21)$$

with the test function $\mathbf{v}^*(\mathbf{x}, z)$ expressed as

$$\mathbf{v}^* = \mathbf{P}^* \circ \mathbf{T}_n + \mathbf{P}_n \circ \mathbf{T}^*. \quad (22)$$

The fixed-point algorithm for solving the resulting nonlinear problem proceeds by calculating functions $\mathbf{T}_n$ by assuming known functions $\mathbf{P}_n$, and then updating the first. The iteration continues until reaching the fixed point.

### 2.2.2 | Carreau fluid

The Stokes model extended to rheothing fluids, like the Carreau model, reads

$$\left\{ \begin{array}{l} \nabla p = \nabla \cdot \mathbf{T} \\ \nabla \cdot \mathbf{v} = 0 \end{array} \right. \qquad (23)$$

where the extra-stress tensor reads

$$\mathbf{T} = 2\eta_{\text{eff}} (D_{\text{eq}}) \mathbf{D}, \quad (24)$$

with the equivalent strain rate $D_{\text{eq}}$ given by

$$D_{\text{eq}} = \sqrt{2(\mathbf{D} : \mathbf{D})}, \quad (25)$$

where “:” denotes the tensor product twice contracted, and the effective viscosity $\eta_{\text{eff}}$ given by the Carreau model

$$\eta_{\text{eff}} = \eta_{\infty} + (\eta_0 - \eta_{\infty}) (1 + (\kappa D_{\text{eq}})^2)^{\frac{m-1}{2}}, \quad (26)$$

with the value of the different parameters considered in our simulations: $\eta_{\infty} = 0$, $\eta_0 = 1$, $\kappa = 100$ and $m = 0.5$.

Again, the solution can be carried out by using a penalty formulation. Moreover, at each iteration of the non-linear solver we must evaluate the equivalent strain rate, then the associated local effective viscosity (from the Carreau model) and write it in a separated form for enhancing the efficiency of the separated representation solver. The simplest way for performing such decomposition

$$\eta_{\text{eff}} \approx \sum_{i=1}^{M} F_i(x, y) \cdot G_i(z), \quad (27)$$

consists of using a singular value decomposition –SVD–.

Thus, the nonlinearity is addressed by reconstructing the nonlinear term and then performing its separated representation by applying the SVD, from the solution (velocity field) available at the previous enrichment step.

## 3 | SEPARABILITY ISSUES AND ROUGHNESS REPRESENTATION

In fact the previous separated representations work perfectly as soon as the domain accepts a separated representation, i.e. $\Omega = \Xi \times (0, H)$, however rough surfaces, with the thickness varying with the coordinate $\mathbf{x}$, i.e. $H(\mathbf{x})$, prevent such a decomposition
of the domain and consequently the use of the solution separated representation. In what follows we distinguish three situations, where without loss of generality, we consider problems defined in 2D domains involving the in-plane coordinate $x \in \Xi = [0, L]$ and the one related to the thickness $y \in I = [0, H(x)]$. Note that from now on, the y-coordinate refers to the domain thickness.

## 3.1 | Non-rectangular domains

To map any non-planar domain $\Omega$ into a reference rectangular domain $\mathcal{R}$ (accepting a trivial separated description as previously discussed), depicted in Fig. 1, we define the mapping $\Omega \to \mathcal{R}: (x, y) \to (r, s)$:

$$ \begin{cases} x = r \\ y = s H(r), \end{cases} \qquad (28) $$

with $(r, s) \in \mathcal{R} = (0, L) \times (0, 1)$.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 Mapping $\Omega \to \mathcal{R}$</figcaption></figure>

The associated Jacobian $J$, that allows writing the problem in the reference domain $\mathcal{R}$ (fully separable), reads

$$ \left\{ \begin{aligned} \frac{\partial x}{\partial r} &= 1 \\ \frac{\partial x}{\partial y} &= 0 \\ \frac{\partial s}{\partial r} &= s H'(r) \\ \frac{\partial s}{\partial y} &= H(r) \end{aligned} \right. . \qquad (29) $$

## 3.2 | Rectangular domains containing a nonplanar interface

Now, we consider the configuration depicted in Fig. 2 in which the nonplanar interface $\Gamma$ defined by $y_I = h(x)$ separates the domain $\Omega = (0, L) \times (0, H)$, $H$ being constant, in general $H \ll L$, in two regions, the upper-domain $\Omega_u$ and the one located at the bottom (with respect to the interface) $\Omega_b$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Rectangular domain $L \times H$ containing a nonplanar interface defined by $y = h(x)$</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Mapping associated with the nonplanar thick interface $h_1(x) < y < h_2(x)$, where domains $\Omega_u$, $\Omega_\Gamma$ & $\Omega_b$ are mapped into the rectangular domains $R_u$, $R_\Gamma$ & $R_b$ respectively</figcaption></figure>

Both sub-domains $\Omega_u$ and $\Omega_b$ can be mapped into the two rectangular domains $R_u = (0, L) \times (1, 2)$ and $R_b = (0, L) \times (0, 1)$ respectively, with $\Omega = \Omega_u \cup \Omega_b$, $R = R_u \cup R_b$ and $\Gamma = \bar{\Omega}_u \cap \bar{\Omega}_b$. As it can be noticed, in the reference domain the interface is given by the horizontal line $s_I = 1$.

Both mappings are expressed by:

1. Mapping $\Omega_b$ into $R_b$. The first mapping with $\mathbf{r} = (r, s) \in R_b = (0, L) \times (0, 1)$ reads

$$ \begin{cases} x = r \\ y = s h(r), \end{cases} \qquad (30) $$

or, equivalently,

$$ \begin{cases} r = x \\ s = \frac{y}{h(x)}. \end{cases} \qquad (31) $$

The components of the Jacobian matrix $\mathbf{J}_b$ result

$$ \begin{cases} \frac{\partial x}{\partial r} = 1 \\ \frac{\partial x}{\partial y} = 0 \\ \frac{\partial s}{\partial r} = s h'(r), \\ \frac{\partial s}{\partial y} = h(r) \end{cases}, \qquad (32) $$

2. Mapping $\Omega_u$ into $R_u$. Equivalently, for $\Omega_u$ we define

$$ \begin{cases} x = r \\ y = (s - 1)(H - h(r)) + h(r), \end{cases} \qquad (33) $$

with $(r, s) \in R_b = (0, L) \times (1, 2)$, that defines its own Jacobian $\mathbf{J}_u$

## 3.3 | Rectangular domain containing a nonplanar thick interface

Sometimes, interfaces become thicker than a simple curve (in 2D) or surface (in 3D) separating two domains. When considering pre-impregnated composite sheets, the sheet surface consists of a thin polymer layer. When two sheets are put in contact, a thin rough polymer layer is created at the interface.

Using the notation previously introduced the thick interface (intermediate layer) $\Omega_\Gamma$ can be defined by its lower and upper boundaries, $h_1(x)$ and $h_2(x)$ respectively, with $h_1(x) \le h_2(x)$, $\forall x \in (0, L)$, $\min h_1(x) > 0$ and $\max h_2(x) < H$.

The mapping (sketched in Fig. 3) that enables transforming $\Omega = \Omega_b \cup \Omega_\Gamma \cup \Omega_u$ into the reference domain composed by layers $R_b = (0, L) \times (0, 1)$, $R_\Gamma = (0, L) \times (1, 2)$ and $R_u = (0, L) \times (2, 3)$, with $R_\Gamma$ the reference domain related to the thick interface $\Omega_\Gamma$ and $R_b$ and $R_u$ respectively the ones related to the plies located at the bottom and top ($\Omega_b$ and $\Omega_u$), reads:

1. Mapping $\Omega_b \to R_b$: $(x, 0 < y < h_1(x)) \to (r, s)$

$$ \begin{cases} x = r \\ y = s h_1(r), \end{cases} \qquad (34) $$
with $(r, s) \in \mathcal{R}_b = (0, L) \times (0, 1)$.

2. Mapping $\Omega_{\Gamma} \to \mathcal{R}_{\Gamma}: (x, h_1(x) \le y \le h_2(x)) \to (r, s)$

$$ \left\{ \begin{array}{l} x = r \\ y = (s - 1)(h_2(r) - h_1(r)) + h_1(r) \end{array} \right. , \qquad (35) $$

with $(r, s) \in \mathcal{R}_{\Gamma} = (0, L) \times (1, 2)$.

3. Mapping $\Omega_u \to \mathcal{R}_u: (x, h_2(x) < y < H) \to (r, s)$

$$ \left\{ \begin{array}{l} x = r \\ y = (s - 2)(H - h_2(r)) + h_2(r) \end{array} \right. , \qquad (36) $$

with $(r, s) \in \mathcal{R}_u = (0, L) \times (2, 3)$.

## 4 | IN-PLANE/OUT-OF-PLANE SEPARATED REPRESENTATION IN THE REFERENCE DOMAIN $\mathcal{R}$

In the previous section, domains with (i) nonplanar surfaces; (ii) nonplanar interfaces, and (iii) nonplanar thick interfaces, were mapped into a reference domain $\mathcal{R}$ where interfaces become planar and parallel to the domain surfaces.

Thus, the domain $\mathcal{R}$ becomes a multilayered domain, fully separable in the in-plane-out-of-plane sense previously discussed. Thus, it becomes natural, for a given problem (differential operator) defined in $\Omega$, transform it into the reference domain $\mathcal{R}$ and use the mapping Jacobian for transforming the differential operators.

Then, the unknown field approximation in the reference coordinates $u(r, s)$ reads

$$ u(r, s) \approx \sum_{i=1}^{N} R_i(r) \cdot S_i(s). \qquad (37) $$

It is important to note that this separated representation only needs solving one-dimensional problems for building-up the fully 2D solution.¹⁶ Thus, one could consider extremely rich discretizations of both dimensions, the plane and the thickness, without compromising the solution efficiency, as discussed in our former works previously referred.

For the sake of completeness, we develop below the procedure when addressing a scalar field, e.g. the temperature (the extension to a vectorial field, e.g. the flow velocity, being straightforward) in the case of a rectangular domain with a nonplanar interface. The unknown field $u(x, y)$ verifies the balance equation

$$ \nabla \cdot (K(\mathbf{x})\nabla u(\mathbf{x})) = 0, \quad \mathbf{x} \in \Omega, \qquad (38) $$

with the following boundary conditions: $u = u_g$ on $\Gamma_D$ ($\Gamma_D \equiv (y = 0) \cup (y = H)$) and $\frac{\partial u}{\partial x} = 0$ in $\Gamma_N$ ($\Gamma_N \equiv (x = 0) \cup (x = L)$). The thermal conductivity is assumed isotropic and constant in the upper and lower domains, expressed respectively by $K(\mathbf{x} \in \Omega_u) = K_u$ and $K(\mathbf{x} \in \Omega_b) = K_b$.

The problem weak form reads

$$ \int_{\Omega} K(\mathbf{x}) \nabla u^* \cdot \nabla u \, d\mathbf{x} = 0, \qquad (39) $$

with $u^*(\mathbf{x} \in \Gamma_D) = 0$.

As previously discussed, we define the two mappings for transforming $\Omega_u$ and $\Omega_b$ into respectively, $\mathcal{R}_u$ and $\mathcal{R}_b$. Then, we define the characteristic function of both subdomains, $\chi_u(\mathbf{r})$ and $\chi_b(\mathbf{r})$:

$$ \chi_u(\mathbf{r}) = \begin{cases} 1 & \text{if } s > 1 \\ 0 & \text{if } s < 1 \end{cases}, \qquad (40) $$

and

$$ \chi_b(\mathbf{r}) = \begin{cases} 1 & \text{if } s < 1 \\ 0 & \text{if } s > 1 \end{cases}, \qquad (41) $$
that allow expressing the conductivity and the differential operators expressed from $\nabla u = \mathbf{B}\nabla_r u$ (where $\nabla$ and $\nabla_r$ are respectively the gradient in the original and reference domain coordinates) as follows:

$$ 
\left\{ 
\begin{aligned}
K(\mathbf{r}) &= \chi_u(\mathbf{r})K_u + \chi_b(\mathbf{r})K_b, \\
\mathbf{B}(\mathbf{r}) &= \chi_u(\mathbf{r})\mathbf{B}_u(\mathbf{r}) + \chi_b(\mathbf{r})\mathbf{B}_b(\mathbf{r}), \\
\det(\mathbf{J}(\mathbf{r})) &= \chi_u(\mathbf{r})\det(\mathbf{J}_u(\mathbf{r})) + \chi_b(\mathbf{r})\det(\mathbf{J}_b(\mathbf{r})),
\end{aligned}
\right. 
\quad (42) $$

such that the weak form in $\mathcal{R} = \mathcal{R}_u \cup \mathcal{R}_b$ reads

$$ \int_{\mathcal{R}} K(\mathbf{r}) (\nabla_r u^*)^T \mathbf{B}^T(\mathbf{r})\mathbf{B}(\mathbf{r}) \nabla_r u \det(\mathbf{J}(\mathbf{r})) d\mathbf{r} = 0. \quad (43) $$

The most appealing property of the proposed mapping is the fact of enabling a separated representation of the Jacobian determinant as well as the terms involved in matrix B. This fact allows the direct use of standard algorithms deeply described in our former works.¹⁶ The fast variation of the different terms of the Jacobian when addressing fractal surfaces or interfaces is finely described by the possibility of using extremely rich discretization in both directions as discussed later, without a significant impact on the computational efficiency.

# 5 | THERMAL PROBLEM AND EFFECTIVE CONDUCTIVITY

## 5.1 | Problem statement

We consider the thermal problem stated in the previous section while enforcing a through-the-thickness thermal flux by pre-scribing $u(x, y = H) = U > 0$ and $u(x, y = 0) = 0$. We assume the existence of a nonplanar thick interface $\Omega_Γ \subset \Omega$, as previously discussed, with the conductivities of the upper ad lower regions $K_u = K_b = K$ and the one of the thick interface noted by $K_Γ \ll K$.

In order to address the issue associated with the very different characteristic dimensions ($L \gg H$) compromising mesh-based discretization, we could consider a separated representation in $\Omega$

$$ u(x, y) \approx \sum_{i=1}^{N} X_i(x) \cdot Y_i(y), \quad (44) $$

but such a choice faces the difficulty related to the presence of the nonplanar thick interface, that implies a hardly separable conductivity field. The conductivity (as any function $f(\mathbf{x})$, $\mathbf{x} \in \Omega$) can always be separated by invoking, for instance, the singular value decomposition –SVD–, i.e., by expressing the conductivity as

$$ K(x, y) \approx \sum_{k=1}^{K} F_k(x) \cdot G_k(y), \quad (45) $$

where the number of terms in that finite sum K informs on the function separability. In our recent work²² it was proved that when such a decomposition involved many terms, the convergence rate of the solution procedure slows down significantly.

Thus, in the case of a planar thick interface defined by $h_1(x) = h_1$ and $h_2(x) = h_2$, with $h_2 > h_1$, $0 < h_1$ and $h_2 < H$, a single term suffices for separating the conductivity, i.e. $K = 1$,

$$ K^p(x, y) = F^p(x) \cdot G^p(y), \quad (46) $$

where the superscript $(\cdot)^p$ refers to the planar interface configuration, $F^p(x) = 1$ and $G^p(y)$ is defined by

$$ G^p(y) = K - (K - K_Γ)\chi_Γ(y), \quad (47) $$

with the characteristic function of the interface expressed from

$$ \chi(y) = \begin{cases} 1 & h_1 < y < h_2 \\ 0 & y \le h_1 \\ 0 & y \ge h_2 \end{cases}. \quad (48) $$

However, when interfaces deviate from the planar configuration, the number of modes K involved in the thermal conductivity separated representation, Eq. (45), increases prohibitively²², and with it the number of operators involved in the weak form (39).
Thus, both the performance and the efficiency of the solver are seriously compromised. In order to circumvent this issue, we consider the mapping $\Omega \to \mathbb{R}$ previously described.

The functions involved in the finite sum decomposition (44) are approximated by using standard one-dimensional finite elements (piecewise linear) and computed using the standard procedures, deeply described in our former works.¹⁶

## 5.2 | Thick interface description

We consider a thick interface $\Omega_Γ$ in a thin domain $\Omega = (0, L) \times (0, H)$, $H \ll L$, emulating a gap in between the two prepreg layers, occupied by the air, with lower thermal conductivity, at the origin of the so-called contact thermal resistance. The thick interface is described by $h_1(x)$ and $h_2(x)$ expressed from

$$h_1(x) = \frac{H}{2} - (A_1 + \delta_1) - A_1 \cos(n_1 x - \phi_1), \quad (49)$$

and

$$h_2(x) = \frac{H}{2} + (A_2 + \delta_2) + A_2 \cos(n_2 x - \phi_2), \quad (50)$$

with $\delta_1 > 0$ and $\delta_2 > 0$ to avoid the layers interpenetration.

## 5.3 | Effective thermal flux

Thermal fluxes across both lateral boundaries are prevented, and the net heat flux is computed on the bottom surface from

$$Q = \int_{0}^{L} K(y=0^{+}) \left. \frac{\partial u}{\partial y} \right|_{y=0^{+}} dx, \quad (51)$$

where $0^+$ refers to the bottom surface neighborhood, within the material domain.

## 5.4 | Parametric analysis

We consider the case study with the domain length $L = 2\pi$, $H = 0.002$ and the conductivities $K_u = K_b = 1$ and $K_Γ = 0.01$ respectively (metric system units). The present analysis aims at evaluating the impact of a thick interface with much smaller thermal conductivity for different geometrical configurations. Following the PGD rationale the conductivity contrast could be introduced as an extra-parameter, however such a route is not retained in the present study. It is important to emphasize that the parametric study that follows does not concern the composites forming practice, it only aims at analyzing: (i) the effect of interface waviness described by the interface profile wavelength (or frequency) that motivated the use of $2\pi$ as in-plane domain length; (ii) the effect of the phase difference between both harmonic profiles (upper and lower) comprising in-between the thick interface; (iii) the effect of extremely important dimension differences –domain length to thickness, $L/H$ (of about four orders of magnitude) and domain to thickness interface thickness, $H/(\delta_1 + \delta_2)$ (of about two orders of magnitude); as well as (iv) significant thermal contrasts –here of two orders of magnitude–. Thus, in our numerical tests it was assumed without loss of generality, $\delta_1 = \delta_2 = 0.01\text{mm}$, $A_1 = A_2 = 0.02\text{mm}$, $n_1 = n_2 = n$ (for better evaluating –without loss of generality– the effect of having the roughness in phase or in phase opposition) and $\phi_1 = 0$. The temperature on the upper boundary $y = H$ was prescribed to a unit value, i.e. $u(x, y = H) = U = 1$, vanishing on the one at the bottom.

Fig. 4 shows for $n = 1$ the temperature field for $\phi_2 = 0$ and $\phi_2 = \pi$ and Fig. 5 reports the parametric net thermal flux $Q(n, \phi)$, with $\phi = \phi_2 - \phi_1 = \phi_2$. From this last figure it can be noticed that the solution is very affected by the phase angle—a constant thin gap ($\phi_2 = \pi$) enhances the heat transfer—as well as by the frequency.

# 6 | FLOW PROBLEM AND APPARENT VISCOSITY

## 6.1 | Flow problem

This section focuses on the effects of surface roughness on the flow for both linear (Newtonian) and non-Newtonian (Carreau) fluids. We consider fractal rough surfaces characterized by its amplitude and fractal dimension. The gap thickness $H$ is defined
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Air gap between two rough surfaces. Temperature field: (left) $\phi_2 = 0$ (in phase); (right) $\phi_2 = \pi$ (in phase opposition)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5 Air gap between two rough surfaces: parametric thermal flux with respect to the roughness frequency and phase angle</figcaption></figure>

from the average value of $H(x)$

$$H = \frac{1}{L} \int_{0}^{L} H(x) \, dx, \qquad (52)$$

and it remains unchanged. The other parameter directly related to the roughness amplitude is $e$, the minimum gap thickness, i.e; $e = \min H(x)$, $\forall x$. All these parameters are sketched in Fig. 6.

Fig. 7 depicts some realizations of the fractal surfaces for different values of its fractal dimension $D$ and $e$, with $L = 10$ and $H = 1$ (dimensionless values chosen to enforce nonlinear effects when considering non-Newtonian fluid flows), for representative a thin rough flow domain, constructed by using the procedure described in$^{23}$ exploiting the resemblance between surface topography and random motion.

In the case of a gap with planar flat surfaces, $e = H$ and $D = 1$. In that case, if we apply a horizontal velocity $V$ to one of the surfaces (e.g. the lower one), while keeping the other (the upper one) at rest, the shear rate will be $\dot{\gamma} = \frac{V}{H}$. Thus, in the case of a linear fluid with viscosity, the shear stress $\tau$ is linearly related to the shear rate from $\tau = \eta\dot{\gamma}$.
<figure><img src="image_7.png" /><figcaption>FIGURE 6 Roughness descriptors</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7 Fractal surfaces with different values of *D* and *e*</figcaption></figure>

The separated representation enables the use of hundreds of millions of degrees of freedom for performing the in-plane approximation and also hundreds of millions of degrees of freedom for defining the through-of-thickness approximation, fact that results in an equivalent 3D resolution involving more than 10¹⁶ degrees of freedom, in our knowledge never attained, and able to represent the flow details at the level of each surface asperity.

The flow problems reported in Section 2.2, both defined in rough thin gaps Ω, some of them depicted in Fig. 7, are transformed into the reference one $\mathcal{R}$, of rectangular shape, by applying the mapping described in Section 3. Then, the standard PGD solver making use of the in-plane-out-of-plane separated representation¹⁶ described in Section 4 applies in the reference rectangular domain.
The velocity field is prescribed on the bottom surface for representing a given horizontal velocity $V$ ($V = 1$ in our simulations) whereas on the top surface, on which the rough surface was mapped, velocity vanishes.

When considering the Carreau fluid the local viscosity must be locally reconstructed and then expressed in a separated form (affine decomposition) by using the so-called “PGD in approximation”¹⁶, that in the considered case is equivalent to the application of the singular value decomposition –SVD–.

The solutions obtained when using the PGD solver were compared with the ones resulting from a standard FEM discretization in the limit of $D$ and $e$ close to the unity (quite flat surface), the one that is attainable with a fully 2D discretization, because when the fractal dimension increases a mesh able to resolve all the flow details at the asperities level becomes too expensive.

After solving the flow problem, the shear rate can be computed at the level of the lower surface, and then by using the local viscosity ($\eta = 1$ in the Newtonian case and the model given in Section 2.2.2 for the Carreau fluid), the shear stress can be computed, that integrated over the domain length, leads to the applied force per unit of length $F$

$$F = \frac{1}{L} \int_{0}^{L} (2\eta(\mathbf{x})\mathbf{D}(\mathbf{x}) \cdot \mathbf{n})_{y=0^{+}} dx, \qquad (53)$$

with $\mathbf{n}^T = (0, 1)$.

Now, if $H$ and $V$ are assumed unchanged (only the roughness effects are being analyzed), and we proceed to solve different flow problems for different values of $e$ and $D$, each solution allows computing the just defined force, whose parametric dependence on $e$ and $D$ is indicated by $F(e, D)$.

## 6.2 | Effective thickness and apparent viscosity

In the linear case (Newtonian fluid) the knowledge of $F(e, D)$ allows defining either an apparent viscosity $\tilde{\eta}$ or an effective gap thickness $\tilde{H}$, according to

1. Apparent viscosity.

In this case from

$$F(e, D) = \tilde{\eta}(e, D) \frac{V}{H}, \qquad (54)$$

we obtain the parametric apparent viscosity

$$\tilde{\eta}(e, D) = F(e, D) \frac{H}{V}. \qquad (55)$$

If $\tilde{\eta} > \eta$, the roughness effect resists the flow, while $\tilde{\eta} \approx \eta$ implies no net effect of the roughness on the flow.

2. Effective thickness.

In this case by assuming that the fluid viscosity is $\eta$, we consider that it is the effective gap thickness that changes, i.e.

$$\tilde{H}(e, D) = \eta \frac{V}{F(e, D)}. \qquad (56)$$

One could naively expect that the effective gap thickness should always be higher than $e$.

In the non-Newtonian setting, only the approach based on the effective gap thickness seems physically meaningful. However, in this case Eq. (56) deserves additional comments. In fact, Eq. (56) becomes nonlinear because $\eta$ depends on the shear rate that depends at its turn on the effective thickness, i.e.

$$\tilde{H}(e, D) = \eta \left( \frac{V}{\tilde{H}(e, D)} \right) \frac{V}{F(e, D)}, \qquad (57)$$

where the viscosity, in the case of a Carreau fluid reads was given in Section Carr.

Note that in the previous equation the viscosity is assumed depending on the effective shear rate that becomes constant everywhere in the flow domain as soon as it takes place in a rectangular domain with a constant effective thickness.

## 6.3 | Parametric analysis

In this section we summarize the results obtained when keeping $H$ constant, $H = 1$, while changing $1.2 \le D \le 1.9$ and $0 < e \le 1$.
<figure><img src="image_9.png" /><figcaption>FIGURE 8 Force vs thickness: Newtonian fluid</figcaption></figure>

### 6.3.1 | Newtonian case

By considering a constant viscosity, the evolution of the force with respect to the effective thickness reads

$$F = \eta \frac{V}{H}, \qquad (58)$$

that is represented in a log-scale in Fig. 8.

Fig. 9 shows the effective thickness with respect to the roughness parameters *e* and *D*. As it can be noticed when *e* approaches *H* (*e* ≈ 1 in the present case), the effective thickness approaches the value of *H* (here *H* = 1), independently on the surface fractal dimension, that always has a smaller effect in this case.

For better appreciation the effects of the fractal dimension, Fig. 10 depicts for *e* = 0.6, the evolution on the effective thickness with respect to the fractal dimension *D*.

### 6.3.2 | Non-Newtonian case

Similar simulations where performed by considering the Carreau fluid. In the present case the evolution of the force with respect to the effective thickness is shown in Fig. 11.

Fig. 12 shows the effective thickness with respect to the roughness parameters *e* and *D*. Again it can be noticed that when *e* approaches *H* the effective thickness approaches the value of *H* (here *H* = 1), independently on the surface fractal dimension that again has a smaller effect.

Again, for better appreciating the effects of the fractal dimension, Fig. 13 depicts for *e* = 0.6 the evolution on the effective thickness with respect to the fractal dimension *D*.

### 6.3.3 | Discussion

The main output of the just described results is that for a given *e* (minimum value of *H*(*x*), ∀*x* ∈ (0, *L*)), the effective thickness can become smaller than *e* depending on the roughness fractal dimension, as Figs. 10 and 13 reveal. Moreover, this effect is more noticeable in the case of the non-Newtonian fluid (it becomes stronger when the fluid power index *m*–related to the Carreau law–decreases exacerbating the non-Newtonian behavior). In all cases, when *e* approaches the value of *e* = 1, roughness disappears and the effective thickness approaches the gap thickness (*H* = 1).
<figure><img src="image_10.png" /><figcaption>FIGURE 9 Effective thickness: Newtonian fluid</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10 Effective thickness ($e = 0.6$) versus fractal dimension $D$: Newtonian fluid</figcaption></figure>

## 7 | CONCLUSIONS

This paper proposes an efficient numerical procedure for evaluating the impact of rough surfaces and interfaces on the effective thermal conductivity and the apparent fluid viscosity (or the effective gap thickness).

The proposed procedure consists of an appropriate mapping able to transform non separable domains into a separable one, on which a high-resolution separated representation enables its solution while circumventing issues related to the domain degeneracy (thicknesses of usual prepegs composites remain much smaller than the in-plane characteristic dimensions).
<figure><img src="image_12.png" /><figcaption>FIGURE 11 Force vs thickness: Non-Newtonian fluid</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12 Effective thickness: Non-Newtonian fluid</figcaption></figure>

The main conclusions of the study being from one side the significant impact of roughness on the thermal conductivity as well as the relative positioning of two roughness surfaces describing a thick rough interface. On the other hand, by increasing the roughness, the effective gap thickness can become smaller than the minimum value of the gap thickness.

These results should allow, from one side, to better defining the process conditions and the optimal material surfaces, and from the other, to make advances towards engineered surfaces and interfaces.
<figure><img src="image_14.png" /><figcaption>FIGURE 13 Effective thickness ($e = 0.6$) versus fractal dimension $D$: Non-Newtonian fluid</figcaption></figure>

## CONFLICT OF INTERESTS STATEMENT

The authors declare that they have no conflict of interest.

## References

1. F. Chinesta, A. Leygue, B. Bognet, Ch. Ghnatios, F. Poulhaon, F. Bordeu, A. Barasinski, A. Poitou, S. Chatel, and S. Maison-Le-Poec. First steps towards an advanced simulation of composites manufacturing by automated tape placement. *International Journal of Material Forming*, 7(1), 81-92, 2014.

2. Ch. Ghnatios, F. Chinesta, Ch. Binetruy. The Squeeze Flow of Composite Laminates. *International Journal of Material Forming*, 8, 73-83, 2015.

3. Ch. Ghnatios, E. Abisset-Chavanne, Ch. Binetruy, F. Chinesta, S. Advani. 3D Modeling of Squeeze Flow of Multiaxial Laminates. *Journal of Non-Newtonian Fluid Mechanics*, 234, 188-200, 2016.

4. G. Yang, R. Pitchumani. A fractal cantor set based description of interlaminar contact evolution during thermoplastic composites processing. *Journal of Materials Science*, 36(19), 4661-4671, 2001.

5. A. Leon, A. Barasinski, E. Nadal, F. Chinesta. High-resolution thermal analysis at thermoplastic pre-impregnated composite interfaces. *Composite Interfaces*, 22(8), 767-777, 2015.

6. W. Lee, G. Springer. A model of the manufacturing process of thermoplastic matrix composites. *Journal of Composite Materials*, 21, 1057-1082, 1987.

7. A. Leon, C. Argerich, A. Barasinski, E. Soccard, F. Chinesta. Effects of material and process parameters on in-situ consolidation. *International Journal of Material Forming*, 2018, https://doi.org/10.1007/s12289-018-1430-7

8. H. Hertz. Ueber die berührung fester elastischer körper. *Journal für die Reine und Angewandte Mathematik*, 92, 156-171, 1881.

9. J. F. Archard. Elastic deformation and the laws of friction. *Proceedings of Royal Society London A: Mathematical, Physical and Engineering Sciences*, 243(1233), 190-205, 1957.
10. J. A. Greenwood, J. B. P. Williamson. Contact of nominally flat surfaces. Proceedings of Royal Society London A: Mathematical, Physical and Engineering Sciences, 295(1142), 300-319, 1966.

11. F. Robbe-Valloire, B. Paffoni, R. Progri. Load transmission by elastic, elasto-plastic or fully plastic deformation of rough interface asperities. Mechanics of Materials, 33(11), 617-633, 2001.

12. B. N. J. Persson. Contact mechanics for randomly rough surfaces. Surface Science Reports, 61(4), 201-227, 2006.

13. G. Zavarise, M. Borri-Brunetto, M. Paggi. On the resolution dependence of micromechanical contact models. Wear, 262(1-2), 42-54, 2007.

14. A. Leon, A. Barasinski, F. Chinesta. Microstructural analysis of pre-impregnated tapes consolidation. International Journal of Material Forming, 10(3), 369-378, 2017.

15. A. Saoudi, A. Leon, G. Gregoire, A. Barasinski, H. Djebaili, F. Chinesta. On the interfacial thermal properties of two rough surfaces in contact in preimpregnated composites consolidation. Surface Topography: Metrology and Properties, 5, 045010, 2017.

16. F. Chinesta, R. Keunings, A. Leygue. The Proper Generalized Decomposition for advanced numerical simulations. A primer. Springerbriefs, Springer, 2013.

17. B. Bognet, A. Leygue, F. Chinesta, A. Poitou, F. Bordeu. Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity. Computer Methods in Applied Mechanics and Engineering, 201, 1-12, 2012.

18. B. Bognet, A. Leygue, F. Chinesta. Separated representations of 3D elastic solutions in shell geometries. Advanced Modelling and Simulation in Engineering Sciences, 2014, 1:4, http://www.amses-journal.com/content/1/1/4

19. L. Gallimard, P. Vidal, O. Polit. Coupling finite element and reliability analysis through proper generalized decomposition model reduction. International Journal for Numerical Methods in Engineering, 95(13), 1079-1093, 2013.

20. P. Vidal, L. Gallimard, O. Polit. Composite beam finite element based on the Proper Generalized Decomposition. Computers & Structures, 102, 76-86, 2012.

21. R. Ibanez, E. Abisset-Chavanne, F. Chinesta, A. Huerta. 3D mixed formulation for simulating squeeze flows in multiaxial laminates. International Journal of Material Forming, 10/5, 653-669, 2017.

22. C. Ghnatios, E. Abisset-Chavanne, A. Ammar, E. Cueto, J.L. Duval, F. Chinesta. Advanced spatial separated representations. Computer Methods in Applied Mechanics and Engineering, 354, 802-819, 2019.

23. K. Komvopoulos. Effects of multi-scale roughness and frictional heating on solid body contact deformation. C. R. Mecanique, 336, 149-162, 2008, 354, 802-819, 2019.