# Advanced Multiscale Modelling of Composite Welding Processes

Francisco Chinesta, Amine Ammar, F. Lemarchand, P. Beauchene, F. Boust

► To cite this version:

Francisco Chinesta, Amine Ammar, F. Lemarchand, P. Beauchene, F. Boust. Advanced Multiscale Modelling of Composite Welding Processes. 9th International Conference on Numerical Methods in Industrial Froming Processes, Numiform' 07, 2007, Porto, Portugal. pp.1-4, 10.1063/1.2740826. hal-00289716

HAL Id: hal-00289716

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

Submitted on 14 Feb 2018

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 CC BY 4.0 - Attribution - International License
# Advanced Multiscale Modelling of Composite Welding Processes

F. Chinesta*, A. Ammar†, F. Lemarchand**, P. Beauchene** and F. Boust**

*LMSP, UMR CNRS-ENSAM-ESEM, 151 Boulevard de l'Hopital, F-75013 Paris, France

†Laboratoire de Rhéologie, INPG-UJF-CNRS, 1301 rue de la piscine, BP 53 Domaine universitaire F-38041 Grenoble Cedex 9, France

**ONERA-DMSC, 29 avenue de la Division Leclerc, F-92320 Châtillon, France

**Abstract.** This work focuses on the solution of thermal models usually encountered in composite plates welding processes that involve a certain number of numerical difficulties related to: (i) the very fine mesh required due to the small domain thickness with respect to the other characteristic dimensions as well as to the presence of a heat source moving on the domain surface; (ii) the long simulation times induced by the small thermal conductivity of polymers, the important thermal solicitations as well as to the movement of the heat source; and (iii) the necessity to define a homogenized thermal conductivity which can vary significantly from one point to other in the domain. In this paper we explore the applicability of a novel high resolution strategy able to solve fast and accurately thermal models involving thousands of millions degrees of freedom. This technique opens new perspectives in computational mechanics, and more specifically in material homogenization.

**Keywords:** Thermal model; Numerical modelling; Homogenization; Separated representation; Tensor product approximation basis

**PACS:**

## INTRODUCTION

Innovative industrial processes dedicated to carbon fibers reinforced composites (CFRP) intends to achieve in a single stage the preform making, through Automated Tape Placement, and its consolidation. In this regard, fibres tapes are placed by a robot meanwhile they are welded to the previous layers, in the case of thermoplastics, or polymerized, in the case of thermosets, by a local heating (with the help of laser, hot gas torch or electronic beam).

These new processes are cost-effective only if they are quick enough which implies very fast heating and cooling and therefore important thermal gradients. These gradients associated to differences between dilatation coefficients of fibres and matrix lead to substantial residual strain and stresses that are even enhanced, in the case of thermoplastics, by variations of specific volumes due to crystallization phenomena.

There are important drawbacks issues. Some of them can be partly remedied, as for example the macroscopic parts distortion by controlling the tape placements; some others as micro-cracks cannot.

These thermal models involved in the numerical modelling of such processes call for the solution of a certain number of numerical difficulties related to: (i) the very fine mesh required due to the small domain thickness with respect to the other characteristic dimensions as well as to the presence of a heat source moving on the domain surface; (ii) the long simulation times induced by the small thermal conductivity of polymers, the impor-

tant thermal solicitations as well as to the movement of the heat source; and (iii) the necessity to define a homogenized thermal conductivity which can vary significantly from one point to other in the domain.

The consequence of the first difficulty is the necessity of using very fine meshes, and therefore a large number of degrees of freedom, with a significant impact on the simulation algorithms efficiency, are involved. Moreover, when the transient problem must be solved incrementally, the extremely large number of time increments induces a second difficulty. For alleviating this drawback we proposed in [3] a parallel time discretization combined with an adaptive reduced modelling.

The other difficulty concerns the material heterogeneity at the microscopic scale. Even if we could in a first approximation neglect the dependence of thermal parameters (specific heat, thermal conductivity ...) on the temperature -linear modelling-, the macroscopic thermal conductivity will depend on the considered point in the macroscopic part, because it depends strongly on the microstructure details (fibers volume fraction, fiber orientations, ...). Due to the small size of fibers, very fine models are required in order to analyze the microscopic patterns, which implies the necessity of using high resolution homogenization techniques. As just discussed, the main drawback comes from the extremely large of degrees of freedom involved in a such microscopic analysis where the level of detail could be in the order of the size of a fiber. In this paper we focus in a novel numerical technique that could be potentially an appealing choice for treating this kind of numerical models.
Thermal Model In The Whole Domain

From now on, the thermal model is assumed defined in the 3D plate domain $\Omega = ]0, L_x[ \times ]0, L_y[ \times ]0, L_z[$ where $L_z \ll L_x$ and $L_z \ll L_y$, in the time interval $[0, t_{max}]$, on whose boundary $\partial\Omega$ different kind of boundary conditions are prescribed, the temperature on $\partial_T \Omega$ and the heat flux on $\partial_Q \Omega$, with $\partial\Omega = \partial_T \Omega \cup \partial_Q \Omega$. The components of $\mathbf{x} \in \Omega$ will be denoted by $(x, y, z)$. Without loss of generality, we consider $\partial_T \Omega = \{\mathbf{x} : z = L_z\}$ and $\partial_Q \Omega$ the complementary part of the domain boundary. The temperature $T(\mathbf{x}, t)$ is then prescribed on $\partial_T \Omega$ to be a given temperature $T_g(\mathbf{x}, t)$. In $\partial_Q \Omega$, and we assume a null heat flux across the boundary, that is $(\mathbf{K}(T, \mathbf{x}) \nabla T) \cdot \mathbf{n} = 0$, where $\mathbf{K}(T, \mathbf{x})$ denotes the homogenized thermal conductivity tensor and $\mathbf{n}$ the unit outwards vector defined on the domain boundary. Moreover, an initial condition must be imposed at the initial time: $T(\mathbf{x}, t=0) = T^0(\mathbf{x})$. The usual variational formulation of the conduction heat equation is expressed by:

Find $T(\mathbf{x}, t \in H^1(\Omega \times ]0, t_{max}))$ verifying the essential boundary conditions $T(\mathbf{x} \in \partial_T \Omega, t) = T_g(\mathbf{x}, t)$, such that

$$ \int_{\Omega} \rho C(\mathbf{x}, T) T^* \frac{\partial T}{\partial t} d\Omega = - \int_{\Omega} \nabla T^* \cdot (\mathbf{K}(\mathbf{x}, T) \nabla T) d\Omega + \\ + \int_{\partial_Q \Omega} T^* (\mathbf{K}(\mathbf{x}, T) \nabla T) \cdot \mathbf{n} d\Gamma + \\ + \int_{\Omega} T^* f(\mathbf{x}, T, t) d\Omega, \quad \forall T^* \in H_0^1(\Omega) \quad (1) $$

where $H^1$ and $H_0^1$ are the usual Sobolev functional spaces, and the specific heat and the homogenized thermal conductivity tensor depend on space (in non homogeneous media) and on the temperature. The source term $f(\mathbf{x}, T, t)$ accounts for possible couplings or phase changes phenomena.

HIGH RESOLUTION SOLVERS BASED ON SEPARATED REPRESENTATIONS

We first consider a representative volume element $\Omega_{rve}$ that for composite materials applications must contain numerous fibers.

Standard homogenization techniques require the solution in such representative volume of three steady state thermal models with three different prescribed boundary conditions. These problems can be written as:

$$ \int_{\Omega_{rve}} \nabla T^* \cdot (\mathbf{k}(\mathbf{x}) \nabla T) d\Omega = 0 \quad (2) $$

with prescribed boundary conditions $T(\mathbf{x} \in \partial\Omega_{rve}) = T_g^i$, where the superscript $i$ refers to each one of these three problems to be solved.

Due to the characteristic size of those fibers a very fine mesh will be required for solving accurately the thermal models related to Eq. (2). Thus, it is usual to define 3D models consisting of thousands of millions of nodes. Nowadays, the performance of computers allow to solve, in a reasonable time, problems rarely exceeding some millions of degrees of freedom (dof).

In order to circumvent this difficulty we propose to use a separated representation in combination with a tensor product reduced approximation basis, whose allows to a discrete model of size $N \times D$ ($N$ being the number of nodes in each space direction and $D$ the dimension of the space) instead the $N^D$ characteristic of grid-based discrete models. Thus, if we suppose $N \approx 10^5$ and $D = 3$, the size of the discrete model will be $3 \times 10^5$ instead $10^{15}$.

In the following paragraphs we summarize the main elements of this technique that we have successfully applied for solving the multidimensional problems usually encountered in the kinetic theory models of complex fluids [1].

We consider the steady conduction heat transfer problem in $\Omega = \Omega_{rve} = ]0, L]^3$, the temperature vanishing on the domain boundary $\partial_T \Omega \equiv \partial\Omega$, $T(\mathbf{x} \in \partial_T \Omega) = 0$ (the non homogeneous case will be addressed later)

$$ \nabla (\mathbf{k}(\mathbf{x}) \nabla T) = -g(\mathbf{x}) \quad (3) $$

where the temperature $T$ and the microscopic thermal conductivity tensor $\mathbf{k}$ depend on the space coordinates $\mathbf{x} = (x_1, x_2, x_3)$. The meaning of the source term $g(\mathbf{x})$ will be addressed later.

The problem solution can be written in the form:

$$ T(x_1, x_2, x_3) = \sum_{j=1}^{\infty} \alpha_j \prod_{k=1}^{3} F_{kj}(x_k) \quad (4) $$

where $F_{kj}$ is the $j$th basis function, with unit norm, which only depends on the $k$th coordinate.

It is well known that the solution of numerous problems can be accurately approximated using a finite (sometimes very reduced) number of approximation functions, i.e.:

$$ T(x_1, x_2, x_3) \approx \sum_{j=1}^{Q} \alpha_j \prod_{k=1}^{3} F_{kj}(x_k) \quad (5) $$

The previous expression implies the same number of approximation functions in each dimension, but each one of these functions could be expressed in a discrete form using different number of parameters (nodes in the 1D grids).

Now, an appropriate numerical procedure is needed for computing the coefficients $\alpha_j$ as well as the $Q$ approximations functions in each dimension.
The proposed numerical scheme consists of an iteration procedure that solves at each iteration n the following two steps:

1. **Projection of the solution in a discrete basis.**

If we assume the functions $F_{kj}(\forall j \in [1, ..., n]; \forall k \in [1, ..., 3])$ known (verifying the boundary conditions), the coefficients $\alpha_j$ can be computed by introducing the approximation of T into the Galerkin variational formulation associated with Eq. (3):

$$ \int_{\Omega} \nabla T^* \cdot (\mathbf{k}(\mathbf{x}) \nabla T) \, d\Omega = \int_{\Omega} T^* g(\mathbf{x}) \, d\Omega \quad (6) $$

Introducing the approximation of T and $T^*$:

$$ T(x_1, x_2, x_3) = \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{3} F_{kj}(x_k) \quad (7) $$

and

$$ T^*(x_1, x_2, x_3) = \sum_{j=1}^{n} \alpha_j^* \prod_{k=1}^{3} F_{kj}(x_k) \quad (8) $$

we have

$$ \begin{aligned} & \int_{\Omega} \nabla \left( \sum_{j=1}^{n} \alpha_j^* \prod_{k=1}^{3} F_{kj}(x_k) \right) \times \\ & \qquad \left( \mathbf{k}(\mathbf{x}) \nabla \left( \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{3} F_{kj}(x_k) \right) \right) \, d\Omega = \\ & = \int_{\Omega} \left( \sum_{j=1}^{n} \alpha_j^* \prod_{k=1}^{3} F_{kj}(x_k) \right) g(\mathbf{x}) \, d\Omega \quad (9) \end{aligned} $$

Now, we assume that $g(\mathbf{x})$ and $\mathbf{k}(\mathbf{x})$ can be written in the form

$$ \left\{ \begin{array}{l} g(x_1, x_2, x_3) = \sum_{h=1}^{m_f} \prod_{k=1}^{3} g_{kh}(x_k) \\ k_{ij}(x_1, x_2, x_3) = \sum_{h=1}^{m_k^j} \prod_{k=1}^{3} (k_{ij})_{kh}(x_k) \end{array} \right. \quad (10) $$

This kind of separated representation can be easily performed using singular value decomposition (in 2D) or using an alternating least squares technique in the general multidimensional case. For this purpose we only need to know the values of the scalar or tensor field) in a cloud of nodes.

Eq. (9) involves integrals of a product of 3 functions each one defined in a different dimension. Let $\prod_{k=1}^{3} h_k(x_k)$ be one of these functions to be integrated. The integral over $\Omega$ can be performed by integrating each function in its 1D definition interval and then multiplying the 3 computed integrals according to:

$$ \int_{\Omega} \prod_{k=1}^{3} h_k(x_k) \, d\Omega = \prod_{k=1}^{3} \int_{0}^{L} h_k(x_k) \, dx_k, \quad (11) $$

which makes possible a fast numerical 1D integration.

Now, due to the arbitrariness of the coefficients $\alpha_j^*$, Eq. (9) allows to compute the n-approximation coefficients $\alpha_j$, solving the resulting linear system of size $n \times n$. This problem is linear and moreover rarely exceeds the order of tens of degrees of freedom. Thus, even if the resulting coefficient matrix is densely populated, the time required for its solution is negligible with respect to the one required for performing the approximation basis enrichment (step 2).

2. **Enrichment of the approximation basis.**

From the coefficients $\alpha_j$ just computed the approximation basis can be enriched by adding the new functions $\prod_{k=1}^{3} F_{k(n+1)}(x_k)$. For this purpose we solve the non-linear Galerkin variational formulation related to Eq. (3):

$$ \int_{\Omega} \nabla T^* \cdot (\mathbf{k}(\mathbf{x}) \nabla T) \, d\Omega = \int_{\Omega} T^* g(\mathbf{x}) \, d\Omega, \quad (12) $$

using the approximation of T given by:

$$ T(x_1, x_2, x_3) = \sum_{j=1}^{n} \alpha_j \prod_{k=1}^{3} F_{kj}(x_k) + \prod_{k=1}^{3} R_k(x_k). \quad (13) $$

whose associated weighting functions result:

$$ T^*(\mathbf{x}) = R_1^*(x_1)R_2(x_2)R_3(x_3) + R_1(x_1)R_2^*(x_2)R_3(x_3) + R_1(x_1)R_2(x_2)R_3^*(x_3) \quad (14) $$

This leads to a non-linear variational problem, whose solution allows to compute the 3 functions $R_k(x_k)$. Functions $F_{k(n+1)}(x_k)$ are finally obtained by normalizing, after convergence of the non-linear problem, the functions $R_1, R_2, R_3$.

To solve this problem we introduce a discretization of those functions $R_k(x_k)$. Each one of these functions is approximated using a 1D finite element description. If we assume than $p_k$ nodes are used to construct the interpolation of function $R_k(x_k)$ in the interval $[0, L]$, then the size of the resulting discrete non-linear problem is $\sum_{k=1}^{k=3} p_k$. The price to pay for avoiding a whole mesh in the domain is the solution of a non-linear problem. However, the size of those non-linear problems remains moderate and no particular difficulties have been found in its solution.
The just proposed technique exhibits four order con-
vergence when the Lebesque norm is considered and the
nodes used as integration points. Its extension to tran-
sient problems is quite simple and direct [2], as well as its
application for solving non-linear models [4]. We come
back to these properties in the next section.

As previously indicated, usual homogenization tech-
niques require the solution of three problems with pre-
scribed essential boundary conditions $T(\mathbf{x} \in \partial_T \Omega) =$
$T_g(\mathbf{x}) = x_i$, with $i = 1, 2, 3$. These problems can be trans-
formed in order to ensure homogeneous boundary condi-
tions by applying the change of variable $\tilde{T}(\mathbf{x}) = T(\mathbf{x}) -
x_i$, which leads to:

$$
\begin{align}
\int_{\Omega} \nabla \tilde{T}^* \cdot (\mathbf{k}(\mathbf{x}) \nabla \tilde{T}) \, d\Omega &= \nonumber \\
&= - \int_{\Omega} \nabla \tilde{T}^* \cdot (\mathbf{k}(\mathbf{x}) \nabla (x_i)) \, d\Omega \tag{15}
\end{align}
$$

whose second member takes the role of function g(x) in
the previous discussion.

However, in order to extend this procedure to more
general boundary conditions, other possibility, under
consideration at present, lies in the use of Lagrange mul-
tipliers or penalization techniques. The first one could
seem more elegant but it introduces a new unknown field
(the Lagrange multiplier) that must be expressed in a sep-
arated representation and then introduced in the resulting
mixed variational formulation. In this case we must ad-
dress also the question related to the LBB stability con-
dition, which could imply a different approximations of
the different fields. To avoid, or at least postpone, the
analysis of this problematic, the most simple and natural
way to enforce general initial or boundary conditions is
by using a penalization technique, that in the case of Eq.
(2) assuming a vanishing source term, results:

$$
\int_{\Omega} \nabla T^* \cdot (\mathbf{k}(\mathbf{x}) \nabla T) \, d\Omega + \alpha \int_{\partial\Omega} T^*(T - T_g^i) \, d\Gamma \quad (16)
$$

where $\alpha$ represents the penalization coefficient. The first
results concerning this alternative as well its application
for accounting for complex domains is an ongoing work.

CONVERGENCE ANALYSIS

In order to conclude about the convergence rates of the
proposed strategy we are considering a steady problem
as well as a transient one whose large time solution
coincides with the solution of the steady problem.

The steady problem is defined by:

$$
\Delta T = -N \prod_{k=1}^{N} \sin(x_k) \quad (17)
$$

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Convergence analysis for the multidimensional heat equation</figcaption></figure>

in $\Omega = ] -\pi, \pi]^N$, with $T(\mathbf{x} \in \partial\Omega) = 0$, whose exact solution reads:

$$
T^{\text{ex}} = \prod_{k=1}^{N} \sin(x_k). \quad (18)
$$

This problem is solved in different multidimensional
spaces considering the separated representation and an
increasing number of nodes for the functional approxi-
mation in each dimension (for the sake of simplicity we
are considering that all the one-dimensional functional
approximations are defined from the same number of
nodes $n_n$).

The error *E* is quantified from the norm:

$$
E = \frac{\int_{\Omega} (T^{\text{num}} - T^{\text{ex}})^2 d\Omega}{|\Omega|}. \quad (19)
$$

Figure 1 depicts the error versus the size of the one-
dimensional space discretization, that in this example re-
sults $h = 2\pi/(n_n - 1)$, when different multidimensional
spaces are considered. A convergence rate of 4 (slope in
the log-log representation) is noticed, two times higher
than the quadratic one characteristic of piecewise linear
approximations.

Despite at present the code is not optimized, we have
reported in figure 2 the CPU time associated to each so-
lution. We can notice that even in 100 dimensions and
using 1000 degrees of freedom in each dimension (note
that the size of the resulting model is of 10⁵ when the
separated representation is used instead of 10³⁰⁰ required
in the framework of the finite element or any other grid
method) the CPU time remains in the order of 1000 sec-
onds. A slope of two is noticed that seems independent of
the space dimension. Thus, the computing time increases
in a power of two with respect to the number of nodes
considered for the one-dimensional functional approxi-
mations.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Computing time against nodal density</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Computing time against space dimension</figcaption></figure>

To check its dependence on the space dimension we depict in figure 3 the evolution of the CPU time versus the space dimension for two meshes consisting of $n_n = 401$ and $n_n = 1601$ nodes (in each direction). The same slope of two is noticed proving that the CPU time increases in a similar way with the space dimension and with the number of degrees of freedom considered for the one-dimensional functional approximations.

## Numerical Example

We consider the steady conduction heat transfer problem defined in a 2D square domain $\Omega = ]0, L[^2$ where the conductivity is assumed isotropic and where the temperature is prescribed on its boundary $\partial_T \Omega = \partial \Omega$ according to $T(x \in \partial_T \Omega) = T_g(x) = x_1$ (usual homogenization strategies solve thermal models with these kind of linear boundary conditions):

$$ \nabla \cdot (k(x) \nabla T) = 0 \qquad (20) $$

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Thermal conductivity distribution.</figcaption></figure>

where the temperature $T$ and the microscopic thermal conductivity tensor $k$ depend on the space coordinates $x = (x_1, x_2)$. The associated variational formulation results

$$ \int_{\Omega} \nabla T^* \cdot (k(x) \nabla T) \, d\Omega = 0 \qquad (21) $$

This problem can be transformed into a one with homogeneous boundary conditions by introducing

$$ \tilde{T}(\mathbf{x}) = T(\mathbf{x}) - x_1 \qquad (22) $$

which implies

$$ \int_{\Omega} \nabla \tilde{T}^* \cdot (k(\mathbf{x}) \nabla \tilde{T}) \, d\Omega = - \int_{\Omega} \nabla \tilde{T}^* \cdot (k(\mathbf{x}) \nabla (x_1)) \, d\Omega \qquad (23) $$

where

$$ (k(\mathbf{x}) \nabla (x_1))^T = (k(\mathbf{x}), 0) \qquad (24) $$

Figure 4 illustrates the conductivity at each point. Finally Fig. 5 depicts the computed temperature field $T(\mathbf{x})$ as well as the one related to the field $\tilde{T}$ depicted in figure 6, where $T_0$ is $T_g(\mathbf{x})$.

## CONCLUSION

This paper explores the ability of separated representations for performing high resolution computations. We have proved that the use of separated representation allows a major reduction of the CPU time. At present this procedure is limited to the linear case but its extension to more complex geometries with general boundary conditions as well as to the homogenization of non-linear behaviors constitute some works in progress.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Computed temperature field.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Computed temperature field.</figcaption></figure>

REFERENCES

1. Ammar A., Mokdad B., Chinesta F., Keunings R., A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, Journal of Non-Newtonian Fluid Mechanics, 139, 153-176, (2006).
2. A. Ammar, B. Mokdad, F. Chinesta, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids. Part 2: Transient simulation using space-time separated representations. Journal of Non-Newtonian Fluid Mechanics, In press.
3. F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, F. Boust, Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization. Computer Methods in Applied Mechanics and Engineering, Submitted.
4. B. Mokdad, E. Pruliere, A. Ammar, F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach. Applied Rheology, In press.