# Ultrasonic Welding of Thermoplastic Composites, Modeling of the Process Using Time Homogenization.

Arthur Lévy, Steven Le Corre, Arnaud Poitou, Eric Soccard

► To cite this version:

Arthur Lévy, Steven Le Corre, Arnaud Poitou, Eric Soccard. Ultrasonic Welding of Thermoplastic Composites, Modeling of the Process Using Time Homogenization. International Journal for Multiscale Computational Engineering, 2011, 9 (6), pp.53-72. 10.1615/IntJMultCompEng.v9.i1.50 . hal-01007405v2

HAL Id: hal-01007405

https://hal.science/hal-01007405v2

Submitted on 9 Apr 2021

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
# Ultrasonic Welding of Thermoplastic Composites, Modeling of the Process Using Time Homogenization.

Arthur LEVY¹, Steven LE CORRE¹, Arnaud POITOU¹, Eric SOCCARD².

March 4, 2009

1. GeM - Ecole Centrale Nantes, Université de Nantes - 1 rue de la Noë 44321 Nantes, France.

e-mail: [arthur.levy; steven.le-corre; arnaud.poitou]@ec-nantes.fr.

2. EADS IW - Technocampus, 1 rue de la Noé 44321 Nantes, France.

e-mail: eric.soccard@eads.net.

## Abstract

The process of ultrasonic welding allows to assemble thermoplastic composite parts. A high frequency vibration imposed to the processing zone induces self heating and melting of the polymer. The main feature of this process is the existence of phenomena that occur on two very different time scales: the vibration (about $10^{-5}$ s) and the flow of molten polymer (about 1 s). In order to simulate accurately these phenomena without the use of a very fine time discretization over the whole process, we apply a time homogenization technique. First, the thermo-mechanical problem is formulated using a Maxwell viscoelastic constitutive law and then, it is homogenized using asymptotic expansion. This leads to three coupled problems: a micro-chronological mechanical problem, a macro-chronological mechanical problem and a macro-chronological thermal problem. This coupled formulation is actually simpler because the macro-chronological problems do not depend on the micro time scale and its associated fast variations. Lastly, a uniform simple test case is proposed to compare the homogenized solution to a direct calculation. It shows that the method gives good results provided that the vibration is fast enough compared to the duration of the process. Moreover, the time saving appears to be highly reduced down to one thousand times less.

**Keywords** : Polymer, Welding, Time-homogenization, Asymptotic Expansion, Viscoelasticity.
Introduction

In the wide field of thermoplastic composite material, the ability of the matrix to
melt allowed new processes to appear. Beside the forming processes, this paper focuses
on assembling processes, among which welding appears as a revolution in the domain
of composites materials. Nevertheless, whereas metal welding implies conductive me-
dia, thermoplastic matrix reveals a very insulating behaviour. This property led to de-
velop new techniques that produce heat very locally, at the welding interface. We can
mention resistance welding or induction welding [Ageorges et al., 2000]. This study
deals with the ultrasonic welding, where a tool, the sonotrode, applies an ultrasonic
vibration to the processing zone at the interface, where little triangles called “energy
directors” are molded to create a kind of controlled rugosity. This induces self-heating
and melting of the energy directors which enable welding.

Although this process has been used for few decades (for polymer sealing for in-
stance), the innovation here consists in the “dynamic” welding. Instead of having a so
called “static” welding where the sonotrode doesn’t move on the plate, in our case, the
sonotrode slides along the directors direction and performs a weld line. The induced
three dimensional flow of the polymer then results in high quality welds.

Many studies [Wang et al., 2006, Benatar and Gutowski, 1989, Tolunay, 1983, Suresh et al., 2007] have been carried out on the modeling of the classical “static process” where a plane strain assumption can be stated for the flow, the simulation of the fully three dimensional dynamic process is much heavier. The main problem that arises when trying to perform this simulation comes from the existence of two time scales. The first one which will be denoted by “micro-chronological” is linked to the short ultrasonic period. The second one, called the“macro-chronological” time scale, is the characteristic duration of the polymer flow at the interface which is the process duration. Simulating each ultrasonic cycle would need very small time steps (about 20 per cycle), over a quite long process time (about 3 s). This would induce a lot of time steps (about a hundred thousand for a 20 kHz load) which would lead to huge calculation times especially in a three dimensional context. A time homogenization method is therefore needed to allow calculation time reductions.

To overcome the existence of two time scales in a given process, many authors [Wang-
et al., 2006, Benatar and Gutowski, 1989, Tolunay, 1983, Suresh et al., 2007] consid-
ered the two mechanical phenomena to be dissociated. Given no flow of polymer, the
ultrasonic vibration of the sonotrode only induces small strains at a high frequency in
the processing zone. Therefore the polymer is assumed to have a linear elastic be-
haviour. This allowed them to calculate a self-heating source independently of the
macro-chronological flow. Nevertheless our process exhibit simultaneously both me-
chanical phenomena.

The same difficulty appears when dealing with structures subjected to fatigue. In
the continuum damage domain, in order to simulate numerically such a multi time
scale process, authors like Van-Paepegem et al. [Van Paepegem et al., 2001] or Cojo-
caru and Karlsson [Cojocaru and Karlsson, 2006] propose the “cycle jump” method.
It consists in solving a micro-chronological problem over a full cycle once for each macro-chronological time step. The main difficulty lies in determining the macro-chronological time-step reasonably. One has to determine how many micro-chronological cycles can be jumped without affecting the global accuracy. In the present paper, considering the two time scales as well separated, we rather propose a time homogenization technique which may avoid the use of cycle jump.

The homogenization techniques first appeared in the field of heterogeneous media with a multi spatial scale problem. In early eighties Benssoussan Lions and Papanicolaou [Benssousan et al., 1978] followed by Sanchez-Palancia [Sanchez-Palencia, 1980] considered a material presenting a periodic structure at the representative volume element scale. The solution is therefore supposed to be periodic in the small scale variable. One can then define a scale factor and seek the solution as a spatial asymptotic expansion in the power of this factor. The initial multiscale problem can then be split into several eventually coupled quite standard sub-problems. Francfort and Suquet [Francfort and Suquet, 1986] confirmed and extended the convergence results for a thermo-visco-elastic material. More recently, we can mention Geindreau and Auriault [Geindreau and Auriault, 1999] for the metallic alloys, Boutin and Auriault [Boutin and Auriault, 1990] for the bituminous concrete or Le Corre et al. [Le Corre et al., 2005] who considered a fibre suspension behaviour in composite forming.

In contrast to spatial homogenization, a time homogenization method can be applied when there exist two separated time scales in a problem. Back to the fatigue problem, for instance, from an initial elastoviscoplastic problem, Guennouni [Guennouni, 1988] obtained two sub problems: an elastic behaviour that describes the micro-chronological effects, and a elastoviscoplastic behaviour for the homogenized macro-chronological behaviour. Since macro-chronological domain evolution and dependency of material parameters are not taken into account, the micro-chronological problem does not depend on the macro chronological time scale. Therefore, it can be solved independently once only. This enables those rather standard problems to be solved very simply. With the same mathematical tool, Oskay and Fish [Oskay and Fish, 2004] described the fatigue phenomena using the continuum damage evolution approach. They obtain two coupled micro and macro-chronological mechanics and a damage evolution law. This was validated on a test case by a numerical FEM implementation. To summarize, the micro-chronological resolution allows to determine the field of interest which is the damage field. Similarly, in our approach, the resolution of the micro-chronological ultrasonic loading allows to determine the heat source and therefore the temperature field.

Closer to our application, we can mention Boutin and Wong [Boutin and Wong, 1998], who treated a coupled visco-elastic and transient thermal problem in two steps. First they applied a space homogenization technique to the mechanical problem and obtained a macroscopic spatial problem. Then, after having discussed the cases where both spatial and temporal scale can be homogenized, they propose an homogenized thermal problem where the source is simply averaged over a micro-chronological period. Note that, they did not use any asymptotic expansion in time. For a fully simulta-
neous time and space homogenization, we can refer to Yu and Fish [Yu and Fish, 2002] who considered the spatial and temporal scale factors to be equal and therefore propose a “double scale asymptotic expansion”.

In the first part of this paper, after a short description of the initial thermo-mechanical problem, a dimensional analysis is carried out in relation with the industrial process orders of magnitude. We then propose a time homogenization technique using asymptotic expansion and obtain three coupled sub-problems (a micro and a macro chronological mechanical problems, and a macro chronological thermal problem). Even if in the scope of the ultrasonic welding process, the methodology proposed in this work aims to be rather general, so that it may be applied to any high frequency processing. In the second part, with the use of an analytical test case, we validate the proposed method by comparing results with a direct resolution.

# 1 Formulation of the Thermo - Mechanical Problem

Let us consider an energy director represented by a domain $\Omega$ as described in figure 1. It is subjected to an evolving small oscillatory displacement on the boundary $\Gamma_u$ (see figure 1). We aim to describe its thermo-mechanical evolution equations. Note that even though the presented figure is two dimensional, the analysis proposed is more general and can be led for a three dimensional problem.

## 1.1 Mechanical problem.

Because of the small dimensions of our domain, we can assume a quasi static description for the mechanics. Neglecting the body force as well, we get:

$$ \nabla \cdot \boldsymbol{\Sigma} = 0, \qquad (1) $$

where $\boldsymbol{\Sigma}$ is the Cauchy stress tensor.

<figure><img src="image_2.png" /><figcaption>Figure 1: Geometry of the problem.</figcaption></figure>
Considering the polymer to be incompressible we also assume:

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

where $\mathbf{u}$ is the displacement field. Therefore, the total stress $\mathbf{\Sigma}$ may be expressed as $\mathbf{\Sigma} = -p\mathbf{I} + \mathbf{\sigma}$ where $p$ is the hydrostatic pressure that can be seen as the Lagrange multiplier associated to incompressibility condition (2) and $\mathbf{\sigma}$ is the extra stress tensor which is determined by given constitutive equations for the polymer. In order to simplify the analysis, a Maxwell law is adopted as a first approach. In this first part, study is restricted to a small displacement framework where the constitutive equations can be written as:

$$ \lambda \dot{\mathbf{\sigma}} + \mathbf{\sigma} = 2\eta \dot{\mathbf{\epsilon}}, \qquad (3) $$

where $\lambda$ is the relaxation time, $\eta$ the viscosity, $\mathbf{\epsilon}$ the strain tensor and the dot expresses the time derivation.

The domain $\Omega$ is subjected to evolving oscillatory displacement boundary conditions on $\Gamma_u$ (see figure 1) which is the sum of a macro-chronological displacement boundary condition $\mathbf{u}_s$ linked to the squeezing of the energy director, and a micro-chronological harmonic boundary condition of amplitude $\mathbf{a}$, linked to the ultrasonic vibration:

$$ \mathbf{u} = \mathbf{u}_s(t) + \mathbf{a}\sin(\omega t) \quad (\text{on } \Gamma_u) \qquad (4) $$

The boundary $\Gamma_\sigma$ is assumed to be free:

$$ \mathbf{\Sigma} \cdot \mathbf{n} = \mathbf{0} \quad (\text{on } \Gamma_\sigma). \qquad (5) $$

**Initial conditions.** At the initial time, we consider the displacement to be zero and a relaxed configuration:

$$ \mathbf{u}(x, t = 0) = \mathbf{0}; \quad \mathbf{\sigma}(x, t = 0) = \mathbf{0}; \quad p(x, t = 0) = 0. \qquad (6) $$

## 1.2 Thermal problem.

$\Omega$ is supposed to be thermally insulated. Considering the elasticity to be fully entropic, which is a classical assumption when dealing with polymer at the rubbery state above the glassy temperature, the internal energy $e$ depends on temperature only and the energy balance can be written as:

$$ \left\{ \begin{array}{l} \rho \frac{de}{d\theta} \dot{\theta} = \rho c \dot{\theta} = \nabla \cdot (k \nabla \theta) + \mathbf{\sigma} : \dot{\mathbf{\epsilon}} \\ k \nabla \theta \cdot \mathbf{n} = 0 \quad \text{on } (\Gamma_u \cup \Gamma_\sigma) \end{array} \right. , \qquad (7) $$

where $\theta$ is the temperature, $k$ the thermal conductivity, $c$ the specific heat capacity and $\rho$ the density.
# 2 Dimensional Analysis

## 2.1 Dimensionless variables

In order to deal with comparable quantities, we have to turn the initial problem defined by equations (1) to (7) into a dimensionless problem. We proceed by using characteristic magnitudes of the process. Note that these values are specific to the ultrasonic welding.

* The characteristic length used in this work is $e = 1 \text{ mm}$ which corresponds to the initial height of the processing zone in our industrial application.

* The time $\lambda_0 = 1 \text{ s}$ which is approximately the duration of the process is used as characteristic time.

* We choose a characteristic stress $\sigma_c = \eta_0/\lambda_0$, where $\eta_0 = 10^7 \text{ Pa.s}$ that is an upper bound of the Newtonian viscosity of the studied polymer (PEEK) close to the glassy temperature $\theta = \theta_g = 143^\circ\text{C}$: $\sigma_c = 10^7 \text{ Pa}$

* Calling $\theta_m$ the melting temperature of the polymer and $\theta_r$ the room temperature, the temperature variation $\Delta\theta_c = \theta_m - \theta_r$ is chosen as the process characteristic temperature variation. In the case of the studied PEEK polymer, $\theta_m \sim 330^\circ\text{C}$ so that $\Delta\theta_c \sim 300^\circ\text{C}$.

We then introduce new dimensionless coordinates:

$$x = e \cdot x^*$$

$$t = \lambda_0 \cdot t^*$$

and finally dimensionless variables $\boldsymbol{\sigma}^*$, $\boldsymbol{u}^*$, $\boldsymbol{\varepsilon}^*$ and $\boldsymbol{\theta}^*$:

$$ \begin{aligned} \boldsymbol{\sigma}(x=ex^*, t=\lambda_0 t^*) &= \boldsymbol{\sigma}_c \cdot \boldsymbol{\sigma}^*(x^*, t^*) \\ \boldsymbol{u}(x=ex^*, t=\lambda_0 t^*) &= e \cdot \boldsymbol{u}^*(x^*, t^*) \\ \boldsymbol{\varepsilon}(x=ex^*, t=\lambda_0 t^*) &= \boldsymbol{\varepsilon}^*(x^*, t^*) \\ \boldsymbol{\theta}(x=ex^*, t=\lambda_0 t^*) &= \Delta\boldsymbol{\theta}_c \cdot \boldsymbol{\theta}^*(x^*, t^*) + \boldsymbol{\theta}_{amb} \end{aligned} \quad (8) $$

Note that the form taken for the dimensionless temperature allows to have it varying between 0 and 1 when it rises from room temperature to melting temperature.

In the same way, dimensionless space derivation operators are introduced:

$$ \begin{aligned} \nabla^*(\cdot) &\equiv e \cdot \nabla(\cdot) \\ \nabla^* \cdot (\cdot) &\equiv e \cdot \nabla \cdot (\cdot) \\ \Delta^*(\cdot) &\equiv \nabla^* \cdot (\nabla^*(\cdot)) = e^2 \Delta(\cdot) \end{aligned} \quad (9) $$
## 2.2 Dimensionless resulting problems

Dimensionless equations corresponding to the general thermo-mechanical problem (equations (1) to (7)) can be deduced<fnref n="1" />:

$$ egin{cases} 
\Lambda \frac{\partial \boldsymbol{\sigma}^*}{\partial t^*} + \boldsymbol{\sigma}^* = 2N \frac{\partial \boldsymbol{\epsilon}^*}{\partial t^*} \\ 
\boldsymbol{\nabla}^* \cdot (\boldsymbol{\sigma}^* - p^* \mathbf{I}) = \mathbf{0} \\ 
\boldsymbol{\nabla}^* \cdot \mathbf{u}^* = 0 \\ 
\mathbf{u}^* = \mathbf{U}_s(t^*) + \mathbf{R} \sin(\omega \lambda_0 t^*) & (\text{on } \Gamma_u) \\ 
(\boldsymbol{\sigma}^* - p^* \mathbf{I}) \cdot \mathbf{n} = \mathbf{0} & (\text{on } \Gamma_\sigma) 
\end{cases} 
\quad (\text{in } \Omega) 
\tag{10} $$

$$ egin{cases} 
\frac{\partial \theta^*}{\partial t^*} = A \Delta^* \theta^* + B \boldsymbol{\sigma}^* : \left( \frac{\partial \boldsymbol{\epsilon}^*}{\partial t^*} \right) & (\text{in } \Omega) \\ 
\boldsymbol{\nabla}^* \boldsymbol{\theta}^* \cdot \mathbf{n} = 0 & (\text{on } \Gamma_\sigma \cup \Gamma_u) 
\end{cases} 
\tag{11} $$

The obtained problem is defined with 6 dimensionless numbers that depend on the boundary conditions and the material parameters:

$$ egin{aligned} 
\mathbf{R} &= a/e \\ 
\mathbf{U}_s &= u_s/e \\ 
A &= (k\lambda_0) / (\rho c e^2) \\ 
B &= \sigma_c / (\rho c \Delta \theta_c) \\ 
\Lambda &= \lambda / \lambda_0 \\ 
N &= \eta / \eta_0 
\end{aligned} 
\tag{12} $$

## 2.3 Scale factor

Let us introduce the dimensionless scale factor which defines how well the two time scales are separated. It is the ratio of a characteristic micro-chronological time and a characteristic macro-chronological one. As already mentioned, the macro-chronological time is $\lambda_0$, the characteristic duration of the squeezing of the energy director. Now, we introduce the micro-chronological time, characteristic of the vibrating loading as $\tau = 1/\omega$, $\omega = 2\pi f$ being the loading pulsation of the sonotrode. The scale factor is finally:

$$ \xi = \frac{1}{\omega\lambda_0}. \tag{13} $$

Considering the wide range of multi time scales phenomena, we can identify several cases. In the domain of high cycle fatigue, for instance, the number of cycles can reach few millions. Therefore, the scale factor may be around $\xi \sim 10^{-9}$. If we focus on ageing phenomena in civil engineering life-time calculation, the number of cycle may be of few tens of thousands. In this case, $\xi \sim 10^{-4}$. In the case of an ultrasonic

<footnote n="1">¹the constitutive law has been divided by $\sigma_c$ and the thermal equation by $\rho c \Delta \theta_c / \lambda_0$</footnote>
process as the considered welding technique, the typical ultrasonic frequency is around $f \sim 20kHz$, so that the scale factor is about:

$$ \xi \sim 10^{-5}. \qquad (14) $$

Scales appear to be well separated in such processes and justifies a homogenization approach.

## 2.4 Evaluation of the dimensionless parameters

As it will be detailed in the following section, the dimensionless parameters (12) have to be evaluated with respect to powers of $\xi$. Different powers of $\xi$ (i.e. different orders of magnitude with respect to the scale separation) can lead to different models.

**Boundary conditions.** The amplitude of the vibrations in ultrasonic welding is of few percents of the energy director's dimension, so that:

$$ R \sim \xi^0 1, \qquad (15) $$

where **1** is a unit displacement vector giving the direction of the loading. It remains true even in the case of very high power ultrasonic processing where the amplitude of the vibration might exceed 50% of the energy director's height. On the contrary, to get an **R** of order of $\xi$, we would need a vibration amplitude of about $10^{-5}$ times the height of the processing zone. Obviously, such a small vibration would not have any effect on the material and is therefore never used in ultrasonic processing.

The macro-chronological displacement $U_s$ is clearly of the order of magnitude of the energy director's height. Therefore, the macro-chronological displacement stays in the vicinity of $\xi^0$ over the whole process:

$$ U_s \sim \xi^0 1. \qquad (16) $$

**Thermal parameters.** The diffusivity of a polymer being around $10^{-7} m^2 s^{-1}$, we get $A \sim 1$ which allows to fix:

$$ A \sim \xi^0. \qquad (17) $$

If we generalize to other materials, the diffusivity ranges between $4 \times 10^{-8}$ for a rubber material and $10^{-4}$ for a copper alloy, so that $0.5 < A < 10^3$. Only for metals with a high thermal diffusivity, $A$ grows to the order of $\xi^{-1}$. The present analysis will be restricted to the case where $A \sim \xi^0$ as in the polymers case.

The specific thermal capacity $\rho c$ of all engineering material are of the order of $10^6 W.m^{-3}K^{-1}$, so that $B \sim 0.03$ with the characteristic stress chosen, and:

$$ B \sim \xi^0. \qquad (18) $$
**Constitutive law.** During the process the temperature of the domain rises from the room temperature to the melting temperature. The material properties ($\lambda$ and $\eta$) being highly temperature dependent, they also vary a lot. Indeed, we can not assume, as Boutin and Wong [Boutin and Wong, 1998], that the temperature does not vary much around the reference temperature, and can not apply a Taylor expansion of the thermo-dependent parameters around this reference temperature. Nevertheless, having limited the study to very well separated time scales ($\xi \sim 10^{-5}$), even if $\lambda$ and $\eta$ vary over few decades, it is quite likely that they will not vary over more than 1 order of $\xi$. We may physically consider 2 case: (i) $\Lambda \sim \xi^0$ and $N \sim \xi^0$ and (ii) $\Lambda \sim \xi^{-1}$ and $N \sim \xi^{-1}$.

Case (i) is strictly valid above the glassy temperature, when the polymer is fully viscoelastic, because the characteristic relaxation time is between 10ms and 100s, and the viscosity is between $10^7$ Pa.s and $10^4$ Pa.s, so that:

$$ \begin{cases} \lambda \sim 1s & \Lambda \sim \xi^0 \\ \eta \sim 5.10^5 Pa.s & N \sim \xi^0 \end{cases} \qquad (19) $$

Case (ii) represents the evolution towards ambient temperature. It is discussed in the appendix where we show that it leads to a similar equation set.

# 3 Homogenization

## 3.1 Time scales

The very noticeable time scale separation justifies the use of a time homogenization approach to model the process [Guennouni, 1988]. To proceed, we first define two independent dimensionless time variables:

* The micro-chronological time scale $\tau^* = \omega t = \omega \lambda_0 t^*$, related to the oscillation, is the expression of the fast variations.

* The macro-chronological time scale $T^* = \frac{t}{\lambda_0} = t^*$, associated with the observation time, is the expression of the process evolution, and in this case, is similar to the dimensionless time already defined.

Since the two times scales are independent, the time derivation of a function $\phi$ that depends on both variables is expanded as follow:

$$ \dot{\phi}(T^*, \tau^*) = \frac{1}{\lambda_0} \frac{\partial \phi}{\partial T^*} + \frac{1}{\xi} \frac{1}{\lambda_0} \frac{\partial \phi}{\partial \tau^*} \qquad (20) $$

## 3.2 Asymptotic expansion

In order to separate the different orders of the physical phenomena in the problem, a variable $\phi$ of the problem is sought as an expansion in power of $\xi$:

$$ \phi(t^*) = \phi_0(T^*, \tau^*) + \phi_1(T^*, \tau^*)\xi + \phi_2(T^*, \tau^*)\xi^2 + \dots, \qquad (21) $$
where $\phi$ stands for $\boldsymbol{\sigma}^*$, $p^*$, $\boldsymbol{u}^*$, $\boldsymbol{\epsilon}^*$ or $\boldsymbol{\theta}^*$. The two time scales being independent, a classical assumption then consists in seeking the solution $\phi_i$ as periodic in $\boldsymbol{\tau}^*$, as in [Guennouni, 1988] or [Boutin and Wong, 1998]. Therefore the global evolution is described by the $T^*$ dependency whereas the periodic vibration is contained in the $\boldsymbol{\tau}^*$ dependency. Doing so, one expects to obtain simpler problems at different powers of $\xi$.

**Velocity formulation** We define the dimensionless velocity:

$$ \boldsymbol{v}^* = \lambda_0 \dot{\boldsymbol{u}}^*. \qquad (22) $$

Starting from the assumption that the expansion of the displacement fields begins at order 0, in our time derivation framework, the velocity expansion starts at order $-1$ with:

$$ \left\{ \begin{array}{l} \boldsymbol{v}_{-1}^* = \frac{\partial \boldsymbol{u}_0^*}{\partial \boldsymbol{\tau}^*} \\ \boldsymbol{v}_i^* = \frac{\partial \boldsymbol{u}_i^*}{\partial T^*} + \frac{\partial \boldsymbol{u}_{i+1}^*}{\partial \boldsymbol{\tau}^*} \quad \forall i \ge 0 \end{array} \right. . \qquad (23) $$

In a similar way, the strain rate tensor is expanded as:

$$ \left\{ \begin{array}{l} \boldsymbol{D}_{-1}^* = \boldsymbol{\nabla}_s^* \boldsymbol{v}_{-1}^* = \frac{\partial \boldsymbol{\epsilon}_0^*}{\partial \boldsymbol{\tau}^*} \\ \boldsymbol{D}_i^* = \boldsymbol{\nabla}_s^* \boldsymbol{v}_{-1}^* = \frac{\partial \boldsymbol{\epsilon}_i^*}{\partial T^*} + \frac{\partial \boldsymbol{\epsilon}_{i+1}^*}{\partial \boldsymbol{\tau}^*} \quad \forall i \ge 0 \end{array} \right. . \qquad (24) $$

Notice that in the following, for the sake of legibility, the stated notation will be dropped, keeping in mind that all variables are dimensionless

### 3.3 Identification of the mechanical problems at first orders.

Once injected in the previous problem, identification of different orders of $\xi$ leads to sets of equations. The equilibrium equations and incompressibility constraint can be identified trivially at every order $i$ of $\xi$. We get:

$$ \forall i \ge 0 \left\{ \begin{array}{l} \nabla \cdot (\boldsymbol{\sigma}_i - p_i \mathbf{I}) = \mathbf{0} \\ \nabla \cdot \boldsymbol{u}_i = 0 \\ (\boldsymbol{\sigma}_i - p_i \mathbf{I}) \cdot \boldsymbol{n} = \mathbf{0} \quad \text{on } \Gamma_{\sigma} \end{array} \right. \qquad (25) $$

which can be written in terms of velocity as:

$$ \forall i \ge -1 \left\{ \begin{array}{l} \nabla \cdot (\boldsymbol{\sigma}_i - p_i \mathbf{I}) = \mathbf{0} \\ \nabla \cdot \boldsymbol{v}_i = 0 \\ (\boldsymbol{\sigma}_i - p_i \mathbf{I}) \cdot \boldsymbol{n} = \mathbf{0} \quad \text{on } \Gamma_{\sigma} \end{array} \right. . \qquad (26) $$
On the contrary, the boundary condition on $\Gamma_u$, given in the system (10), appears at order 0:

$$ \left\{ \begin{array}{l} \boldsymbol{u}_i = 0 \quad \forall i > 0 \\ \boldsymbol{u}_0 = \boldsymbol{u} = \boldsymbol{U}_s(T) + \boldsymbol{R} \sin(\boldsymbol{\tau}) \end{array} \right. \quad (\text{on } \Gamma_u), \qquad (27) $$

which can be written in terms of velocity as:

$$ \left\{ \begin{array}{l} \boldsymbol{v}_{-1} = \frac{\partial \boldsymbol{u}_0}{\partial \boldsymbol{\tau}} = \boldsymbol{R} \cos(\boldsymbol{\tau}) \\ \boldsymbol{v}_0 = \frac{\partial \boldsymbol{u}_1}{\partial \boldsymbol{\tau}} + \frac{\partial \boldsymbol{u}_0}{\partial T} = \frac{d\boldsymbol{U}_s}{dT} = \boldsymbol{v}_s \quad (\text{on } \Gamma_u) \\ \boldsymbol{v}_i = \mathbf{0} \quad \forall i > 0 \end{array} \right. \qquad (28) $$

The expansion of the constitutive law (10) leads to the following set of constitutive equations:

• at order -1 of $\xi$:

$$ \Lambda \frac{\partial \boldsymbol{\sigma}_0}{\partial \boldsymbol{\tau}} = 2N \frac{\partial \boldsymbol{\varepsilon}_0}{\partial \boldsymbol{\tau}} = 2N \boldsymbol{D}_{-1}. \qquad (29) $$

• at order 0 of $\xi$:

$$ \Lambda \frac{\partial \boldsymbol{\sigma}_0}{\partial T} + \Lambda \frac{\partial \boldsymbol{\sigma}_1}{\partial \boldsymbol{\tau}} + \boldsymbol{\sigma}_0 = 2N \frac{\partial \boldsymbol{\varepsilon}_0}{\partial T} + 2N \frac{\partial \boldsymbol{\varepsilon}_1}{\partial \boldsymbol{\tau}} = 2N \boldsymbol{D}_0. \qquad (30) $$

We now introduce the time average operator $\langle \cdot \rangle$ defined as $\langle \cdot \rangle = \frac{1}{|\kappa|} \int_{\kappa}(\cdot) d\tau$, $\kappa$ being the ultrasonic period. The periodicity of functions $\phi_i$ in asymptotic expansion of type (21) enables to write the following general property:

$$ \left\langle \frac{\partial \phi_i}{\partial \tau} \right\rangle = 0 \qquad (31) $$

where $\phi$ stands for $\boldsymbol{\sigma}$ or $\boldsymbol{\varepsilon}$. Applying the averaging process to equation (30) leads to:

$$ \Lambda \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} + \langle \boldsymbol{\sigma}_0 \rangle = 2N \frac{\partial \langle \boldsymbol{\varepsilon}_0 \rangle}{\partial T} = 2N \langle \boldsymbol{D}_0 \rangle. \qquad (32) $$

Then, combining constitutive laws (29) and (32), equilibrium, incompressibility and boundary conditions (26) and (28), we can write two mechanical problems detailed hereunder.

The first one only deals with short time $\tau$. It is the micro-chronological mechanical problem:

$$ \left\{ \begin{array}{l} \Lambda \frac{\partial \boldsymbol{\sigma}_0}{\partial \boldsymbol{\tau}} = 2N \boldsymbol{D}_{-1} \\ \boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}_0 - p_0 \boldsymbol{I}) = \mathbf{0} \\ \boldsymbol{\nabla} \cdot \boldsymbol{v}_{-1} = 0 \\ \boldsymbol{v}_{-1} = \boldsymbol{R} \cos(\boldsymbol{\tau}) \quad (\text{on } \Gamma_u) \\ (\boldsymbol{\sigma}_0 - p_0 \boldsymbol{I}) \cdot \boldsymbol{n} = \mathbf{0} \quad (\text{on } \Gamma_\sigma) \end{array} \right. \qquad (33) $$
It is an hypo-elastic problem which is equivalent to an elastic problem in the small displacement framework. Indeed, for short time loading (i.e. fast displacement), the Maxwell model is mainly elastic. This problem describes a small amplitude oscillatory stress of order $\xi^0$ that is linked to a velocity of order $\xi^{-1}$. Integrating this micro-chronological problem in $\tau$ we get the system:

$$ 
\left\{ 
\begin{array}{l}
\boldsymbol{\sigma}_0 = \frac{2N}{\Lambda} \boldsymbol{\varepsilon} + \boldsymbol{\sigma}_0(T, \tau = 0) \\
\boldsymbol{\nabla} \cdot (\boldsymbol{\sigma}_0 - p_0 \mathbf{I}) = \mathbf{0} \\
\boldsymbol{\nabla} \cdot \mathbf{u} = 0 \\
\mathbf{u} = R \sin(\tau) \quad (\text{on } \Gamma_u) \\
(\boldsymbol{\sigma}_0 - p_0 \mathbf{I}) \cdot \mathbf{n} = \mathbf{0} \quad (\text{on } \Gamma_\sigma)
\end{array}
\right. 
\quad (\text{in } \Omega)
\quad (34) $$

in which the constitutive equation shows that $\sigma_0$ can be decomposed into two terms: a micro-chronological evolution $\frac{2\Lambda}{N}\boldsymbol{\varepsilon}$ and the integration constant $\sigma_0(T, \tau = 0)$ that describes the macro-chronological evolution. Since the boundary condition on $\Gamma_u$ reduces to a micro-chronological harmonic boundary condition and problem is linear, the micro chronological displacement $\mathbf{u}$ can be sight as harmonic: $\mathbf{u} = \tilde{\mathbf{u}}.\sin(\tau)$ and solved in terms of $\tilde{\mathbf{u}}$, which does not depend on $\tau$ any more. Then, $\boldsymbol{\varepsilon}$ being a sinus, averaging the first equation of the system gives $\langle \boldsymbol{\sigma}_0 \rangle = \boldsymbol{\sigma}_0(T, \tau = 0)$ which confirms that $\boldsymbol{\sigma}_0(T, \tau = 0)$ is the macro-chronological evolution of $\boldsymbol{\sigma}_0$.

$\langle \boldsymbol{\sigma}_0 \rangle$ is determined using the second problem, which deals with long time $T$. It is the macro-chronological mechanical homogenized problem obtained from equations (32), (26) and (28):

$$ 
\left\{ 
\begin{array}{l}
\Lambda \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} + \langle \boldsymbol{\sigma}_0 \rangle = 2N \langle \boldsymbol{D}_0 \rangle \\
\boldsymbol{\nabla} \cdot (\langle \boldsymbol{\sigma}_0 \rangle - \langle p_0 \rangle \mathbf{I}) = \mathbf{0} \\
\boldsymbol{\nabla} \cdot \langle \mathbf{v}_0 \rangle = 0 \\
\langle \mathbf{v}_0 \rangle = \mathbf{v}_s(T) \quad (\text{on } \Gamma_u) \\
(\langle \boldsymbol{\sigma}_0 \rangle - \langle p_0 \rangle \mathbf{I}) \cdot \mathbf{n} = \mathbf{0} \quad (\text{on } \Gamma_\sigma)
\end{array}
\right. 
\quad (\text{in } \Omega)
\quad (35) $$

with the initial condition:

$$ \langle \boldsymbol{\sigma}_0 \rangle (\mathbf{x}, T = 0) = \mathbf{0} \quad (36) $$

It describes the slow mechanical variations as a viscoelastic Maxwell flow. The stress average $\langle \boldsymbol{\sigma}_0 \rangle = \boldsymbol{\sigma}_0(T, \tau = 0)$ that allows to complete the problem (34) is also of order $\xi^0$ but is linked to a velocity of order $\xi^0$. The boundary condition on $\Gamma_u$ is the macro-chronological velocity condition only.

## 3.4 Identification of the thermal problem at first orders.

In the thermal problem, the boundary condition of insulation of the system (11) is trivially identified at every order of $\xi$ as:

$$ \forall i \quad \boldsymbol{\nabla} \theta_i \cdot \mathbf{n} = 0 \quad (\Gamma_u \cup \Gamma_\sigma). \quad (37) $$
The constitutive equation (11) can be identified at order $\xi^{-1}$ as:

$$ \frac{\partial \theta_0}{\partial \tau} = B \boldsymbol{\sigma}_0 : \mathbf{D}_{-1}. \qquad (38) $$

Using the first equation of the system (33) we can write:

$$ \langle \boldsymbol{\sigma}_0 : \mathbf{D}_{-1} \rangle = \left\langle \boldsymbol{\sigma}_0 : \frac{\partial \boldsymbol{\sigma}_0}{\partial \tau} \right\rangle = \frac{1}{2} \left\langle \frac{\partial \boldsymbol{\sigma}_0^2}{\partial \tau} \right\rangle = 0 \qquad (39) $$

because $\boldsymbol{\sigma}_0^2$ is $\tau$-periodic. Therefore the micro-chronological thermal equation (38) will not induce global temperature evolution, but will only describe a periodic fast variation. In order to describe a macro-chronological evolution, we need to identify higher orders. At order $\xi^0$ we get:

$$ \frac{\partial \theta_0}{\partial T} + \frac{\partial \theta_1}{\partial \tau} = A \Delta \theta_0 + B \boldsymbol{\sigma}_0 : \mathbf{D}_0 + B \boldsymbol{\sigma}_1 : \mathbf{D}_{-1}, \qquad (40) $$

which can be averaged as:

$$ \frac{\partial \langle \theta_0 \rangle}{\partial T} = A \Delta \langle \theta_0 \rangle + B \langle \boldsymbol{\sigma}_0 : \mathbf{D}_0 \rangle + B \langle \boldsymbol{\sigma}_1 : \mathbf{D}_{-1} \rangle. \qquad (41) $$

Back to the identification of $\mathbf{D}_0$ and $\mathbf{D}_{-1}$ given in equations (29) and (30), we can develop:

$$ \langle \boldsymbol{\sigma}_0 : \mathbf{D}_0 \rangle + \langle \boldsymbol{\sigma}_1 : \mathbf{D}_{-1} \rangle = \left\langle \boldsymbol{\sigma}_0 : \left( \frac{\Lambda}{2N} \frac{\partial \boldsymbol{\sigma}_0}{\partial T} + \frac{\Lambda}{2N} \frac{\partial \boldsymbol{\sigma}_1}{\partial \tau} + \frac{1}{2N} \boldsymbol{\sigma}_0 \right) \right\rangle + \left\langle \boldsymbol{\sigma}_1 : \frac{\Lambda}{2N} \frac{\partial \boldsymbol{\sigma}_0}{\partial \tau} \right\rangle. \qquad (42) $$

Periodicity of $\boldsymbol{\sigma}_i$ given by the assumption (21) then gives:

$$ \left\langle \boldsymbol{\sigma}_0 : \frac{\partial \boldsymbol{\sigma}_1}{\partial \tau} \right\rangle + \left\langle \boldsymbol{\sigma}_1 : \frac{\partial \boldsymbol{\sigma}_0}{\partial \tau} \right\rangle = \left\langle \frac{\partial (\boldsymbol{\sigma}_0 : \boldsymbol{\sigma}_1)}{\partial \tau} \right\rangle = 0 \qquad (43) $$

so that:

$$ \langle \boldsymbol{\sigma}_0 : \mathbf{D}_0 \rangle + \langle \boldsymbol{\sigma}_1 : \mathbf{D}_{-1} \rangle = \frac{1}{2N} \left\langle \boldsymbol{\sigma}_0 : \left( \Lambda \frac{\partial \boldsymbol{\sigma}_0}{\partial T} + \boldsymbol{\sigma}_0 \right) \right\rangle. \qquad (44) $$

The stress $\boldsymbol{\sigma}_0$ can then be decomposed according to the first equation of the system (34) in its macro and micro chronological parts:

$$ \langle \boldsymbol{\sigma}_0 : \mathbf{D}_0 \rangle + \langle \boldsymbol{\sigma}_1 : \mathbf{D}_{-1} \rangle = \frac{2N}{\Lambda^2} \left\langle \boldsymbol{\varepsilon} : \left( \Lambda \frac{\partial \boldsymbol{\varepsilon}}{\partial T} + \boldsymbol{\varepsilon} \right) \right\rangle + \frac{1}{2N} \langle \boldsymbol{\sigma}_0 \rangle : \left( \Lambda \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} + \langle \boldsymbol{\sigma}_0 \rangle \right). \qquad (45) $$
Finally we obtain the following macro-chronological thermal problem:

$$ egin{cases} 
\frac{\partial \langle \theta_0 \rangle}{\partial T} = A \Delta \langle \theta_0 \rangle + \frac{2BN}{\Lambda^2} \left\langle \boldsymbol{\varepsilon} : \left( \Lambda \frac{\partial \boldsymbol{\varepsilon}}{\partial T} + \boldsymbol{\varepsilon} \right) \right\rangle + \frac{B}{2N} \langle \boldsymbol{\sigma}_0 \rangle : \left( \Lambda \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} + \langle \boldsymbol{\sigma}_0 \rangle \right) \\ 
\nabla \theta_0 \cdot \boldsymbol{n} = 0 \quad (\Gamma_u \cup \Gamma_\sigma) \\ 
\langle \theta_0 \rangle (\boldsymbol{x}, T=0) = 0 
\end{cases} 
\quad (46) $$

## 3.5 Summary

The three problems obtained thanks to the homogenization technique are coupled, which brings up the following comments:

* In spite of being defined on the small time scale $\tau$, the micro-chronological mechanical problem depends on the large time scale $T$ through:

- the term $\boldsymbol{\sigma}_0(T, \tau = 0)$ of equation (34),

- the geometry that evolves as the director is squeezed,

- and the evolution of the mechanical parameters with the temperature.

Therefore, in the context of a time integration scheme, it needs to be solved over one cycle at each macro-chronological time step.

* The source term of the macro-chronological thermal problem does not depend on the small time $\tau$. Indeed, even if the micro chronological mechanical strain $\boldsymbol{\varepsilon}$ appears in the source term of the equation (46), we only need the average $\left\langle \boldsymbol{\varepsilon} : \left( \Lambda \frac{\partial \boldsymbol{\varepsilon}}{\partial T} + \boldsymbol{\varepsilon} \right) \right\rangle$ that can be post processed once $\boldsymbol{\varepsilon}$ is known over a whole micro-period.

* The macro-chronological mechanical problem does not depend directly on the small time scale $\tau$. Nevertheless it is weakly coupled to the thermal problem through the dependency of material parameters with the temperature.

Finally we are able to propose an integration scheme where each resolution may be performed over a realistic time discretization. It is illustrated in the flow diagram given in figure 2.

As a conclusion, it is worth underlining the advantage of the method proposed over more pragmatic ones like those proposed by [Wang et al., 2006, Benatar and Gutowski, 1989, Tolunay, 1983, Suresh et al., 2007]. It is systematic and rigorous and requires no other assumption than the scale separation to get the resulting formulation. In particular, the original source term of the thermal problem (46) would not have been obtained without the use of the asymptotic homogenization technique. Indeed, the identification at each order of $\xi$ allows to discriminate the needed macro chronological term from the
<figure><img src="image_3.png" /><figcaption>Figure 2: Integration scheme proposed to solve the three coupled problems.</figcaption></figure>
fluctuations. Moreover, the method clearly gives the order of magnitude of the fields involved in the problems. Namely, the stress is of order 0 in both mechanical problems, but is induced by velocity fields of order −1 in the micro chronological mechanical problem (33) and of order 0 in the macro chronological mechanical problem (35).

# 4 Validation test

In order to illustrate the presented method, we apply it to an uniform test case. Independently from the material consideration, we analyze the convergence of the mathematical forms obtained as $\xi$ becomes very small. Ignoring the spatial dimension, we focus on a homogeneous test which enables a direct and simple numerical solving of the full thermo-mechanical problem given by equations (10) and (11) by using very small time steps. This will allow a comparison between the direct and the proposed homogenized resolution.

Let us consider a rectangular domain $\Omega$ of height $h(t)$ and width $b(t)$ as shown in figure 3. The material is supposed to follow a Maxwell law (3) and we assume a plane strain kinematics. We impose a fluctuating squeezing boundary condition on the upper face $\Gamma_{sup}$ as described on figure 4, such as:

$$h(t) = h_0 \cdot \exp(-\frac{t}{t_s}) + a \cdot \sin(\omega t), \qquad (47)$$

where $t_s$ is a characteristic time of the squeezing.

The velocity boundary condition on $\Gamma_{sup}$ becomes:

$$v_d = \dot{h} = (v_s(t) + a\omega \cdot \cos(\omega t)) \cdot e_y, \qquad (48)$$

with:

$$v_s(t) = -\frac{h_0}{t_s} \exp(-\frac{t}{t_s}). \qquad (49)$$

Considering the domain $\Omega$ to be thermally insulated, we fulfill the conditions of the previous part and therefore can apply our homogenization technique.

## 4.1 Direct resolution.

Given such boundary conditions, kinematics is homogeneous and:

$$D = \frac{v_d}{h} (-e_x \otimes e_x + e_y \otimes e_y). \qquad (50)$$

Considering the Maxwell law (3), we get 6 scalar differential equations among which four obviously leads to;

$$\sigma_{ij} = \sigma_{33} = 0 \quad \forall i \neq j, \quad \forall t. \qquad (51)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Geometry of the test case and boundary conditions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Fluctuating boundary condition on the upper face.</figcaption></figure>
Moreover, $\sigma_{11}$ and $\sigma_{22}$ follow two opposite differential equations:

$$ 
\begin{aligned}
\lambda \dot{\sigma}_{11} + \sigma_{11} &= 2\eta \frac{v_d}{h} \\
\lambda \dot{\sigma}_{22} + \sigma_{22} &= -2\eta \frac{v_d}{h}
\end{aligned}
$$

with the same initial condition $\sigma_{11}(t=0) = \sigma_{22}(t=0) = 0$ so that $\sigma_{11} = -\sigma_{22}$.

Defining the scalar $\sigma = \sigma_{11}$, the extra-stress $\sigma$ can be written as:

$$ \sigma = \sigma(t) (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y), \qquad (52) $$

with:

$$ \lambda \dot{\sigma} + \sigma = 2\eta \frac{v_d(t)}{h(t)}; \qquad \sigma(t=0) = 0. \qquad (53) $$

We restrict the study to an idealized test case where material parameters are supposed to be thermally independent. This allows to solve the mechanics and the thermal problems independently. The temperature field follows the partial differential equation (7). But the stress and strain fields being homogeneous in $\Omega$, the thermal dissipation $\sigma: D$ is also homogeneous, and because of the adiabatic boundary conditions, the temperature field itself is homogeneous in the domain. The temperature finally follows the ordinary differential equation:

$$ \rho c \dot{\theta} = \boldsymbol{\sigma} : \boldsymbol{D} = 2\sigma \frac{v_d}{h}. \qquad (54) $$

## 4.2 Use of the homogenization technique.

**Micro chronological problem:** Rewriting the micro-chronological problem (34) in terms of this test notation, we get:

$$ 
\left\{ 
\begin{array}{l}
\boldsymbol{\sigma}_0 = \frac{2N}{\lambda} \boldsymbol{\varepsilon} + \boldsymbol{\sigma}_0(T, \tau=0) \\
\nabla \cdot (\boldsymbol{\sigma}_0 - p_0 \mathbf{I}) = \mathbf{0} \\
\nabla \cdot \mathbf{u} = 0 \\
\mathbf{u} = a \cdot \sin(\tau) \mathbf{e}_y \quad (\text{on } \Gamma_{\text{sup}}) \\
(\boldsymbol{\sigma}_0 - p_0 \mathbf{I}) \cdot \mathbf{n} = \mathbf{0} \quad (\text{on } \Gamma_{\text{lat}})
\end{array}
\right.
\quad (\text{on } \Omega) \qquad (55) $$

Considering the kinematics of the problem, $\boldsymbol{\varepsilon}$ can be written as:

$$ \boldsymbol{\varepsilon} = \frac{\tilde{a}}{h} (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y) \qquad (56) $$

where $\tilde{a} = a \cdot \sin(\tau)$. The extra-stress $\sigma_0$ is therefore determined as soon as $\sigma_0(T, \tau=0)$ is known.
**Macro chronological problem:** As for the direct resolution, considering the kinematic of the macro-chronological squeezing, we can write:

$$ \langle \mathbf{D}_0 \rangle = \frac{v_s}{h} (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y). \quad (57) $$

The constitutive law of the macro-chronological problem (35) becomes:

$$ \Lambda \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} + \langle \boldsymbol{\sigma}_0 \rangle = 2N \frac{v_s}{h} (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y), \quad (58) $$

and the extra-stress $\langle \boldsymbol{\sigma}_0 \rangle$ can be searched in the direction of $\langle \mathbf{D}_0 \rangle$, as in the direct resolution:

$$ \langle \boldsymbol{\sigma}_0 \rangle = \boldsymbol{\sigma}_M (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y), \quad (59) $$

where $\boldsymbol{\sigma}_M$ follows the ordinary differential equation:

$$ \Lambda \frac{\partial \boldsymbol{\sigma}_M}{\partial T} + \boldsymbol{\sigma}_M = 2N \frac{v_s}{h} \quad \boldsymbol{\sigma}_M(T=0) = 0. \quad (60) $$

**Thermal resolution:** Let us now consider the thermal problem (46). Using equation (55), we remind that $\boldsymbol{\sigma}_0$ is the sum of a micro-chronological fluctuating stress and a macro-chronological stress:

$$ \boldsymbol{\sigma}_0 = \frac{2N}{\Lambda} \boldsymbol{\varepsilon} + \boldsymbol{\sigma}_0(T, \tau=0), \quad (61) $$

which implies that $\langle \boldsymbol{\sigma}_0 \rangle = \boldsymbol{\sigma}_0(T, \tau=0)$ because the micro-chronological strain dependence in $\tau$ is periodic. Therefore,

$$ \boldsymbol{\sigma}_0 = (\boldsymbol{\sigma}_{m}(T, \tau) + \boldsymbol{\sigma}_{M}(T)) (-\mathbf{e}_x \otimes \mathbf{e}_x + \mathbf{e}_y \otimes \mathbf{e}_y) $$

where

$$ \boldsymbol{\sigma}_m = \frac{2N}{\Lambda} \frac{\tilde{a}(\tau)}{h(T)}. \quad (62) $$

This helps us to develop the average term needed to calculate the heat generation:

$$ \langle \boldsymbol{\sigma}_0 : \left( \Lambda \frac{\partial \boldsymbol{\sigma}_0}{\partial T} + \boldsymbol{\sigma}_0 \right) \rangle = 2 \langle \boldsymbol{\sigma}_m^2 \left( -\Lambda \frac{\dot{h}}{h} + 1 \right) \rangle + 2 \boldsymbol{\sigma}_M \left( \Lambda \frac{\partial \boldsymbol{\sigma}_M}{\partial T} + \boldsymbol{\sigma}_M \right) \quad (63) $$

As in the direct resolution, the temperature field is homogeneous, $\Delta\theta = 0$ and the partial differential equation (46) becomes an ordinary differential equation:

$$ \frac{\partial \theta_0}{\partial T} = \frac{B}{N} \langle \boldsymbol{\sigma}_m^2 \left( 1 - \Lambda \frac{\dot{h}}{h} \right) \rangle + \frac{B}{N} \boldsymbol{\sigma}_M \left( \Lambda \frac{\partial \boldsymbol{\sigma}_M}{\partial T} + \boldsymbol{\sigma}_M \right) \quad (64) $$
<table><caption>Table 1: Parameters used in the test case.</caption><tr><td>η</td><td>3.10<sup>4</sup> Pa.s</td></tr><tr><td>λ</td><td>0.1 s</td></tr><tr><td>ρc</td><td>2.610<sup>6</sup> WK<sup>-1</sup>m<sup>-3</sup></td></tr><tr><td>h<sub>0</sub></td><td>300 μm</td></tr><tr><td>t<sub>s</sub></td><td>4 s</td></tr><tr><td>a</td><td>40 μm</td></tr></table>

<table><tr><td>Case A</td><td>Case B</td><td>Case C</td></tr><tr><td>f = 1 Hz</td><td>f = 100 Hz</td><td>f = 10000 Hz</td></tr><tr><td>ω = 6.3 rad.s<sup>-1</sup></td><td>ω = 630 rad.s<sup>-1</sup></td><td>ω = 6.310<sup>4</sup> rad.s<sup>-1</sup></td></tr><tr><td>ξ ~ 0.16</td><td>ξ ~ 1.610<sup>-3</sup></td><td>ξ ~ 1.610<sup>-5</sup></td></tr></table>

<figure><img src="image_6.png" /><figcaption>(a) f = 1 Hz, ξ = 0.16: low scale separation</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) f = 100 Hz, ξ = 0.0016: good scale separation</figcaption></figure>

<caption>Figure 5: Comparison of stresses calculated by the direct and homogeneous method.</caption>

## 4.3 Results

Since this case is homogeneous, solving directly the initial fluctuating problem is possible even for quite high values of the pulsation $ω$. This was done using Mat lab's Runge-Kutta (4,5) algorithm. Considering our industrial material (PEEK thermo-plastic) around 200°C, we adapted DMA modulus values found in the literature [Li, 1999, Goyal et al., 2006, Benatar and Gutowski, 1989] and got order of magnitudes of $λ$ and $η$. The calculation of the direct and homogeneous cases were performed on three different test cases. Case A and B were done over 3 seconds using two different pulsations $ω$, and case C using frequency of our process, which is very high. Therefore the resolution of this case was only performed over 0.5 s. The parameters used are summarized in table 1.

**High value of ξ:** The functions $σ$ of equation (53), representing the stresses computed directly, and the function $σ_m + σ_M$, representing the stresses calculated using the
<figure><img src="image_8.png" /><figcaption>Figure 6: Stress and strain on first periods computed by the direct method.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of temperatures calculated directly and with the presented method</figcaption></figure>
homogenization method, are compared on figure 5.

Analyzing figure 5(a) we notice that for a high value of ξ (here, 0.16), when the two time scales are not well separated, the homogenized stress does not fit the one calculated directly. This can partly be explained by the fact that the micro-chronological problem is calculated as elastic, whereas the direct problem involves quite long time fluctuation (ω is around $10\,\text{rad}\cdot\text{s}^{-1}$ compared to the characteristic Maxwell time of 0.1s). Therefore, the impact of the fluctuation on the Maxwell material are not elastic only. Indeed, as shown on figure 6(a), the fluctuating stress is not in phase with the strain.

Concerning the temperature, the directly calculated one was compared to the homogenized macro chronological temperature calculated from equation (64). Figure 7(a) illustrates their differences. With a final temperature relative error around 150%, our homogenization method is clearly not relevant when the scale factor is to big.

**Low values of ξ:** Increasing the pulsation allows to better separate scales. This is illustrated on case B, where the fluctuation frequency is increased by a factor of 100 (see table 2).

The time scale factor then drops to about $10^{-3}$ and the homogenization technique exhibits good results. The two stresses are represented on Figure 5(b). Because of the high frequency, only their envelopes, and their macro-chronological part (periodic average, for the direct stress $\sigma$) were represented. The two curves match almost perfectly. The relevance of the homogenized solution can be explained by the nearly elastic response of the material to the fluctuation. This is visible on the stress and strain representation over the first periods, on figure 6(b), which are totally in phase.

Furthermore, figure 7(b) shows the perfect matching of directly calculated temperature averaged over each micro-chronological period, and macro-chronological homogenized temperature. The relative difference between the final homogenized and directly calculated temperature drops to about 2%. When $\xi$ is much lower, the homogenization method therefore works very well.

**Very high frequency:** Because of the sharp conditions imposed by the very high frequency (10 kHz), the direct problem was solved on $[0, 0.5s]$ only. Figure 7(c) shows a good matching between direct and homogenized temperature evolutions. The final temperature relative error is now less than 1%. In such case, which is close to more complex industrial situations, the time integration on longer times is already difficult and an accurate solution is very difficult to obtain by the direct calculation. The homogenized solution even seems more relevant and eventually more precise if applied to more complex geometry.

**Time saving:** Concerning the calculation times, in case B, the present method is more than ten times faster than the direct method. This is of course explained by the conditions imposed to the Runge Kutta solver, which are much smoother than the one
<table><caption>Table 2: Comparison of the efficiency of the direct and presented method.</caption><thead><tr><th colspan="2">CPU time (s)</th><th>Case A</th><th>Case B</th><th>Case C</th></tr></thead><tbody><tr><td colspan="2">time scale factor ξ</td><td>0.16</td><td>1.610<sup>-3</sup></td><td>1.610<sup>-5</sup></td></tr><tr><td rowspan="2">CPU time (s)</td><td>direct</td><td>0.39</td><td>10.51</td><td>219</td></tr><tr><td>homogenized</td><td>0.17</td><td>0.16</td><td>0.12</td></tr><tr><td colspan="2">Speed up</td><td>2.3</td><td>65</td><td>1820</td></tr><tr><td colspan="2">Final temperature relative error</td><td>160%</td><td>1.8%</td><td>0.8%</td></tr></tbody></table>

imposed initially in the direct problem. The convergence is therefore obtained rapidly, since the time discretization is coarser. In case C, the time saving increases to a factor of 1800 (see table 2). Indeed, the homogenized solution calculation does not depend on the frequency whereas the direct resolution is calculated using a time discretization that is directly linked to the frequency.

**Self heating:** In addition to validating the efficiency of the presented method with regard to the calculation time saving, this uniform test case allows to put forward a conclusion concerning the process. Observing figure 7(c), one can evaluate the temperature increase around 30°C per seconds. Nevertheless, even in the static case, where the process lasts no more than few seconds, in the real process, the temperature increase is much higher (around 300°C/s). Thus it allows the fusion of the director and the flow of molten polymer. This illustrates the importance of the energy director shape. Indeed the heating, initiated at the tip of the director, would not be possible in the case of rectangular directors and bulk heating. In other words, energy directors have an important role of heat concentrator and the heating and melting are local phenomena, at least in the initial stage.

# 5 Conclusion

Unlike most studies where vibration and flow are considered as independent, we here consider a thermo-mechanical problem including micro and macro chronological phenomena. After having analyzed characteristic dimensions of the ultrasonic problem trying to stay rather general, we applied an asymptotic expansion method based on the introduction of two independent time scales. The method leads to three coupled problems depending on both time scales. It is worth noticing that those two time scales can be treated as independent and be discretized separately. Instead of imposing very small time steps on the large time scale in order to capture the fluctuations, these theoretical results show that we can solve the micro-chronological problem on one vibration period only, once for each macro-chronological time step. In a global simulation of the process, the total number of time iterations can therefore be reasonable.

In order to validate the present method and illustrate its efficiency, it was compared to the solution of the direct problem on a uniform test case. It was shown to give
good results even for quite high values of $\xi$. Furthermore, even in this homogeneous test, the resolution of the direct problem at 10kHz was very long because of the sharp conditions induced by such a high frequency. The present method is therefore naturally more suited in this case. Concerning the time saving, the efficiency is convincing: at 10kHz, the calculation time was shown to be divided by one thousand. One has to highlight the interest of such a time saving if applied to a 3 dimensional nonlinear problem.

Furthermore, the test case brings a first conclusion about the process: the shape of the energy director is essential to the process. It induces highly heterogeneous fields at the tip of the energy director which allow the process to initiate. With a simple rectangular cross section, the bulk heating would not be important enough. To further investigate this last point, a finite element simulation tool handling geometrical evolutions is currently under development.

Finally, thanks to a rigorous mathematical framework, this paper reinforces the approaches of the literature consisting in modelling straightaway the process with three coupled problems. Nevertheless, the small displacement theory may be considered as rather unrealistic in the case of our process, and the present work should be extended to finite transformations. In this case new assumptions will undoubtedly be required.

## A Homogenization with low temperature conditions

In this case, the order of magnitude of the Maxwell relaxation time and Maxwell viscosity of the initial problem (10) are changed. Let us consider that the orders of magnitude (19) become:

$$ \begin{cases} \lambda \sim 10^5 s & \Lambda \sim \xi^{-1} \\ \eta \sim 10^{11} Pa.s & N \sim \xi^{-1} \end{cases}, \qquad (65) $$

which would be the case when the temperature drops to ambient temperature.

Defining $\Lambda_0 = \Lambda\xi \sim \xi^0$ and $N_0 = N\xi \sim \xi^0$, the dimensionless constitutive mechanical equation (10) then becomes:

$$ \Lambda_0 \xi^{-1} \frac{\partial \boldsymbol{\sigma}}{\partial T} + \Lambda_0 \xi^{-2} \frac{\partial \boldsymbol{\sigma}}{\partial \tau} + \boldsymbol{\sigma} = 2N_0 \xi^{-1} \frac{\partial \boldsymbol{\varepsilon}}{\partial T} + 2N_0 \xi^{-2} \frac{\partial \boldsymbol{\varepsilon}}{\partial \tau}, \qquad (66) $$

which identification leads to the two following mechanical problems.

**Micro-chronological problem** Identification of order -2 in $\xi$ gives

$$ \Lambda_0 \frac{\partial \boldsymbol{\sigma}_0}{\partial \tau} = 2N_0 \frac{\partial \boldsymbol{\varepsilon}_0}{\partial \tau} \qquad (67) $$

which is exactly the same as the previous case given by equations (33)
**Macro-chronological problem** Identification of order -1 in $\xi$ gives

$$ \Lambda_0 \frac{\partial \boldsymbol{\sigma}_0}{\partial T} + \Lambda_0 \frac{\partial \boldsymbol{\sigma}_1}{\partial \tau} = 2N_0 \frac{\partial \boldsymbol{\varepsilon}_0}{\partial T} + 2N_0 \frac{\partial \boldsymbol{\varepsilon}_1}{\partial \tau} = 2N_0 \boldsymbol{D}_0 \quad (68) $$

which after averaging gives the macro-chronological constitutive equation:

$$ \Lambda_0 \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} = 2N_0 \frac{\partial \langle \boldsymbol{\varepsilon}_0 \rangle}{\partial T}. \quad (69) $$

At low temperature, the macro-chronological viscoelastic problem (35) turns into an hypo-elastic problem:

$$ \left\{ \begin{array}{l} \Lambda_0 \frac{\partial \langle \boldsymbol{\sigma}_0 \rangle}{\partial T} = 2N_0 \langle \boldsymbol{D}_0 \rangle \\ \boldsymbol{\nabla} \cdot (\langle \boldsymbol{\sigma}_0 \rangle - \langle p_0 \rangle \boldsymbol{I}) = \boldsymbol{0} \\ \boldsymbol{\nabla} \cdot \langle \boldsymbol{v}_0 \rangle = 0 \\ \langle \boldsymbol{v}_0 \rangle = \boldsymbol{v}_d(T) \quad (\text{on } \Gamma_u) \\ (\langle \boldsymbol{\sigma}_0 \rangle - \langle p_0 \rangle \boldsymbol{I}) \cdot \boldsymbol{n} = \boldsymbol{0} \quad (\text{on } \Gamma_\sigma) \end{array} \right. \quad (\text{on } \Omega) \quad (70) $$

**Thermal problem:** The thermal analysis remains the same, the temperature is of order 0 in $\xi$, and the macro-chronological evolution is still given by averaging equation (40):

$$ \frac{\partial \langle \theta_0 \rangle}{\partial T} = A \Delta \langle \theta_0 \rangle + B \langle \boldsymbol{\sigma}_0 : \boldsymbol{D}_0 \rangle + B \langle \boldsymbol{\sigma}_1 : \boldsymbol{D}_{-1} \rangle. \quad (71) $$

Nevertheless, in this case, $\boldsymbol{D}_0$ is given by equation (68). Applying the same method as in section 3.4, we finally get a macro chronological thermal problem of the form:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \langle \theta_0 \rangle}{\partial t} = A \Delta \langle \theta_0 \rangle + \frac{B}{2N} \left\langle \boldsymbol{\sigma}_0 : \left( \Lambda_0 \frac{\partial \boldsymbol{\sigma}_0}{\partial tT} \right) \right\rangle . \\ \displaystyle \nabla \theta_i \cdot \boldsymbol{n} = 0 \quad (\Gamma_u \cup \Gamma_\sigma) \end{array} \right. \quad (72) $$

The three obtained problems are therefore very similar to those of the first case discussed in main text. The viscoelastic contributions simply disappears. A general treatment over the whole temperature range with temperature dependent parameters can be envisaged with the first case formulation, which appears to be the richest one.

## References

[Ageorges et al., 2000] Ageorges, C., Ye, L. et Hou, M., 2000. Advances in fusion bonding techniques for joining thermoplastic matrix composites: a review. *The University of Sydney*, The University of Sydney.

[Benatar and Gutowski, 1989] Benatar, Avraham et Gutowski, Timothy G., 1989. Ultrasonic welding of peek graphite apc-2 composites. *Polymer Engineering and Science*, 29(23):1705.
[Benssousan et al., 1978] Benssousan, Alain, Lions, Jacques Louis et Papanicoulau, George. *Asymptotic analysis for periodic structures*. , volume 5. North Holland, 1978.

[Boutin and Auriault, 1990] Boutin, C. et Auriault, JL, 1990. Dynamic behaviour of porous media saturated by a viscoelastic fluid. application to bituminous concretes. *International journal of engineering science*, 28(11):1157–1181.

[Boutin and Wong, 1998] Boutin, Claude et Wong, Henri, 1998. Study of thermosensitive heterogeneous media via space-time homogenisation. *European Journal of Mechanics - A/Solids*, 17(6):939–968.

[Cojocaru and Karlsson, 2006] Cojocaru, D. et Karlsson, A. M., 2006. A simple numerical method of cycle jumps for cyclically loaded structures. *International Journal of Fatigue*, 28(12):1677–1689.

[Francfort and Suquet, 1986] Francfort, G.A. et Suquet, P.M., 1986. Homogenization and mechanical dissipation in thermoviscoelasticity. *Archive for Rational Mechanics and Analysis*, 96(3):265–293.

[Geindreau and Auriault, 1999] Geindreau, C. et Auriault, J.L., 1999. Investigation of the viscoplastic behaviour of alloys in the semi-solid state by homogenization. *Mechanics of Materials*, 31(8):535–551.

[Goyal et al., 2006] Goyal, R.K., Tiwari, A. N., Mulik, U. P. et Negi, Y. S., 2006. Dynamic mechanical properties of al2o3/poly(ether ether ketone) composites. *Journal of Applied Polymer Science*, 104(1):568–575.

[Guennouni, 1988] Guennouni, Tajeddine, 1988. Sur une méthode de calcul de structures soumises à des chargements cycliques: l'homogénéisation en temps. *Modélisation mathématique et analyse numérique*, 22(3):417–455.

[Le Corre et al., 2005] Le Corre, S., Dumont, P., Orgeas, L. et Favier, D., 2005. Une approche multi-échelle de la rhéologie des suspensions concentrées de fibres. *Revue des Composites et des Matériaux Avancés*, 15(3):355.

[Li, 1999] Li, Rongzhi, 1999. Time-temperature superposition method for glass transition temperature of plastic materials. *Materials Science and Engineering A*, 278(1-2):36–45.

[Oskay and Fish, 2004] Oskay, Caglar et Fish, Jacob, 2004. Fatigue life prediction using 2-scale temporal asymptotic homogenization. *International Journal for Numerical Methods in Engineering*, (61):329–359.

[Sanchez-Palencia, 1980] Sanchez-Palencia, E., 1980. Non-homogeneous media and vibration theory. *Lecture Notes in Physics*, 127.
[Suresh et al., 2007] Suresh, KS, Rani, M.R., Prakasan, K. et Rudramoorthy, R. -, 2007. Modeling of temperature distribution in ultrasonic welding of thermoplastics for various joint designs. *Journal of Materials Processing Tech.*, 186(1-3):138–146.

[Tolunay, 1983] Tolunay, MN, 1983. Heating and bonding mechanisms in ultrasonic welding of thermoplastics. *Polymer Engineering and Science*, 23(13):726–733.

[Van Paepegem et al., 2001] Van Paepegem, W., Degrieck, J. et De Baets, P., 2001. Finite element approach for modelling fatigue damage in fibre-reinforced composite materials. *Composites Part B*, 32(7):575–588.

[Wang et al., 2006] Wang, X., Yan, J., Li, R. et Yang, S., 2006. Fem investigation of the temperature field of energy director during ultrasonic welding of peek composites. *Journal of Thermoplastic Composite Materials*, 19(5):593.

[Yu and Fish, 2002] Yu, Qing et Fish, Jacob, 2002. Multiscale asymptotic homogenization for multiphysics problems with multiple spatial and temporal scales: a coupled thermo-viscoelastic example problem. *International Journal of Solids and Structures*, 39(26):6429–6452.