1.

Introduction

In today’s consumer electronics, displays with ever-increasing pixel density are required. Within the production of such displays, fine metal masks (FMMs) with the highest quality have to be fabricated. Conventionally, the holes are etched into the metal foil; however, this process is reaching its limit as the defects accumulate for pixel densities around 600 PPI.¹ A technology that allows for such quality is high power ultrashort pulsed lasers.² The main reason for resorting to USP lasers is their ability to insert a lot of energy into the material without too much heating, which would destroy such thin foils.² To meet the productivity requirements, the laser can be combined with a multibeam scanner.³ During the production of FMMs, it is important to avoid distortions or discoloration due to heat accumulation.⁴ Therefore, process parameters and processing strategies that avoid such problems are sought.

Within this treatise, an extension to a simulation that computes the structuring of a single borehole, which has been explained in a prior publication,⁵ is described. This time-consuming computation of a single borehole shape is used as the basis for a multiscale model with which the simulation of larger workpieces becomes feasible. The multiscale model consists of two scales.^6–10 First, the heat conduction and deformations are computed for a single borehole, the microscale. Second, the information on the microscale is collected in a representative volume element (RVE) with which the heat conduction and deformations of a whole workpiece, the macroscale, can be computed. The microscale model assumes pulse durations from femtoseconds up to nanoseconds as long as the ablation is not melt-driven.¹¹

The unit cell, which is the computational domain in the microscale simulation, is the direct result of the simulation described in the prior publication.⁵ Both the heat conduction and the thermo-elastic structural mechanics problems are discretized using a Bubnov–Galerkin method.^12–14 Averaging the results over the unit cell domain yields the material properties of the RVE.

On the macroscale, the computational domain is an FMM. It consists of two regions: one for a hole and one for the solid material. It is essential to use the RVE in the region that accounts for a hole and the ordinary material properties of the solid everywhere else. Again, a Bubnov–Galerkin scheme is employed to compute the temperature field and the distortions of the FMM.^12–14

This treatise presents the multiscale approach using a simplified hole geometry. The homogenization of both models, the heat-conduction and the thermo-elasticity models, is performed using an asymptotic expansion in the scale variable as described in Ref. 6. To account for multibeam scanners, the two-scale model can be adapted easily. This novel approach is capable of simulating the temperature field and elastic distortions occurring during the structuring of thin metal foils. Plastic strains and melt formation are not considered. The multiscale simulation computes accurate results in a fraction of the time compared with a fully-resolved simulation.

2.

Multiscale Model

In a prior publication, the ablation of a structure with a single beam was explained.⁵ This treatise extends the results to a multiscale model with which a periodic ablation pattern can be simulated for larger domains in a reasonable time.

The multiscale approach follows closely the works of Fish [8, Chapter I.3] or [6, Chapter 2.2]. To further reduce the computational demands, a residual-free method that allows for computing the microstructure off-line and in advance of the macroscale simulation was developed.

2.1.

Geometric Set-Up and Notation

Hereafter, it is assumed that two length scales that differ in size enough to be separated exist, i.e., $0<\zeta ≔\frac{{\ell }_{\mathrm{f}}}{{\ell }_{\mathrm{c}}}\ll 1$. Variables referring to a particular length scale are denoted with a subscript $x_{\mathrm{f}}$ for fine and $x_{\mathrm{c}}$ for coarse. Furthermore, the notation for functions that account for information on all scales uses a preceding superscript ${}^{\zeta }f$ and functions that are valid on the $I$’th scale use ${}^{I}f$. If not indicated otherwise, the variable $x$ lives on the coarse scale, whereas the variable $y≔x/\zeta$ lives on the fine scale. The computational domain is denoted as $\mathrm{\Omega }$ with the boundary $\partial \mathrm{\Omega }$. To account for different physical interactions with the environment, it is necessary to split this boundary into a Neumann part ${\mathrm{\Gamma }}_N$ and a Dirichlet part ${\mathrm{\Gamma }}_D$, such that $\partial \mathrm{\Omega }={\mathrm{\Gamma }}_N\cup {\mathrm{\Gamma }}_D$ and ${\mathrm{\Gamma }}_N\cap {\mathrm{\Gamma }}_D=0$ (c.f. Fig. 1). Analogously, the domain and boundaries for the unit cell are defined, except the symbol $\mathrm{\Theta }$ is used. Whenever a distinction between Dirichlet and Neumann Boundary has to be made, ${\mathrm{\Gamma }}_D$ and ${\mathrm{\Gamma }}_N$, respectively, are used.

Fig. 1

(a) An artificial unit cell geometry is shown. In reality, the hole shape depends on the ablation strategy and the laser beam properties. (b) Displays a foil created using the unit cell in (a), although the hole sizes are emphasized.

The reason the artificial unit cell geometry [c.f. Fig. 1(a)] is used throughout this treatise is twofold. First, the computation of a realistic hole is already described in a prior publication,⁵ and the geometries can be easily exchanged. Second, on the large scale, the difference between the artificial unit cell and a realistic unit cell is negligible.

In addition to the geometric notation, the Einstein summation convention is employed, i.e., repeated indices are summed. The derivative of a function $f(x_i)$ with respect to $x_i$ is sometimes abbreviated with ${\partial }_{x_i}f$ to emphasize vector field operators such as the divergence. Averaged quantities are denoted with an overline, e.g., $\overline{x}$, vectors and matrices are highlighted in a boldface $\mathbf{x}$ and uppercase boldface $\mathbf{X}$, respectively. Finally, the subscript $f_{(i,x_j)}$ is used to denote the symmetrized gradient tensor entry, i.e.

Eq. (1)

$$f_{(i,x_j)}≔\frac{1}{2}(\frac{\partial f_i}{\partial x_j}+\frac{\partial f_j}{\partial x_i}).$$

More complex geometries can be created using a level set function $\mathrm{\Phi }(\mathbf{x}):{\mathbb{R}}^n\to \mathbb{R}$. Such a level set function is the key result in the companion publication⁵ wherein $\mathrm{\Phi }$ is the solution that describes the ablation surface of a laser structuring process. A unit cell domain can be defined using this function on voxels by removing all elements that lie completely inside the ablated domain, i.e., for which $\mathrm{\Phi }(\mathbf{x})<0$ holds. This can be seen in Fig. 2. The second result of the simulation described in Ref. 5 is the heat distribution after the ablation process is done. This heat distribution can be used as a volume source in Sec. 2.2.

Fig. 2

A level set function can be used to mark elements that lie inside an ablated domain [c.f. Fig. 2(a)]. The final unit cell after removal of marked elements is shown in Fig. 2(b).

2.2.

Two-Scale Heat Conduction

The heat conduction equation for a material with specific density $\rho :{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}$, specific heat capacity $c_p:{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}$, and thermal conductivity $\lambda :{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}$ heated by a heat source $Q:{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}$ reads

Eq. (2)

$$\begin{gathered} \rho c_p\frac{\partial {}^{\zeta }T}{\partial t}+\frac{\partial {}^{\zeta }q_i}{\partial x_i}={}^{\zeta }Q\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }q_i=-\lambda \frac{\partial {}^{\zeta }T}{\partial x_i}\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }q_i{}^{\zeta }n_i={}^{\zeta }q_N\text{on}{}^{\zeta }\mathrm{\Gamma }_N \\ {}^{\zeta }T={}^{\zeta }T_D\text{on}{}^{\zeta }\mathrm{\Gamma }_D. \end{gathered}$$

Herein, $n_i$, refers to the $i$’th entry of the outward pointing unit normal $n$, e.g. $n={(n_1,n_2,n_3)}^T$ in three dimensions. The solution to this problem is the temperature field $T:{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}$. In general, the functions depend on both the coarse and fine scales, but for the density and the specific heat capacity a dependence on the fine scale variable only is assumed, i.e., $\rho ≔\rho (\mathbf{y})$ and $c_p≔c_p(\mathbf{y})$.

Now, an ansatz function for the temperature is defined by means of an asymptotic expansion¹⁵ using the scale $\zeta$, which reads

Eq. (3)

$${}^{\zeta }T(\mathbf{x})≔T(\mathbf{x},\mathbf{y})={\zeta }^0{}^{0}T(\mathbf{x},\mathbf{y})+{\zeta }^1{}^{1}T(\mathbf{x},\mathbf{y})+O({\zeta }^2)+\cdots .$$

Note that terms of order 2 and above have been omitted as they do not contribute significantly to the temperature due to the condition $\zeta \ll 1$. Fourier’s law in Eq. (2) relates the temperature gradient to the heat flux; therefore, the spatial derivative of Eq. (3) is taken. Rearranging the terms in accordance to the order of $\zeta$ yields

Eq. (4)

$$\begin{gathered} \frac{{\partial }^{\zeta }T}{\partial x_i}(\mathbf{x})=\frac{\partial T}{\partial x_i}(\mathbf{x},\mathbf{y})={\zeta }^{-10}T(\mathbf{x},\mathbf{y})+{\zeta }^0(\frac{{\partial }^{0}T}{\partial x_i}(\mathbf{x},\mathbf{y})+\frac{{\partial }^{1}T}{\partial y_i}(\mathbf{x},\mathbf{y}))+{\zeta }^1\frac{{\partial }^{1}T}{\partial x_i}(\mathbf{x},\mathbf{y}) \\ +O({\zeta }^2)+\cdots . \end{gathered}$$

Again, higher-order terms are omitted in the two-scale formulation. Multiplying this equation with the negative thermal conductivity and a comparison of the coefficients results in

Eq. (5)

$$\begin{gathered} -\lambda \frac{\partial T}{\partial x_i}(\mathbf{x},\mathbf{y})=-{\zeta }^{-1}\lambda {}^{0}T(\mathbf{x},\mathbf{y})-{\zeta }^0\lambda (\frac{\partial {}^{0}T}{\partial x_i}(\mathbf{x},\mathbf{y})+\frac{\partial {}^{1}T}{\partial y_i}(\mathbf{x},\mathbf{y}))-{\zeta }^1\lambda \frac{{\partial }^{1}T}{\partial x_i}(\mathbf{x},\mathbf{y})-O({\zeta }^2)-\cdots \\ {q}_i(\mathbf{x})={\zeta }^{-1}{{}^{-1}q}_i(\mathbf{x},\mathbf{y})+{\zeta }^0{{}^{0}q}_i(\mathbf{x},\mathbf{y})+{\zeta }^1{{}^{1}q}_i(\mathbf{x},\mathbf{y})+O({\zeta }^2)+\cdots . \end{gathered}$$

The divergence of the heat flux is needed, too. Taking the derivative with respect to the $i$’th space component and sorting in ascending orders of $\zeta$ yields

Eq. (6)

$$\frac{\partial q_i}{\partial x_i}={\zeta }^{-2}\frac{{\partial }^{-1}{q}_i}{\partial y_i}+{\zeta }^{-1}(\frac{{\partial }^{-1}{q}_i}{\partial x_i}+\frac{{\partial }^0{q}_i}{\partial y_i})+{\zeta }^0(\frac{{\partial }^0{q}_i}{\partial x_i}+\frac{{\partial }^1{q}_i}{\partial y_i})+{\zeta }^1(\frac{{\partial }^1{q}_i}{\partial x_i}+\frac{{\partial }^2{q}_i}{\partial y_i})+O({\zeta }^2)+\cdots .$$

Here, the functions’ arguments are dropped for brevity. Inserting the ansatz functions—Eq. (3) for the temperature and Eq. (5) for the heat flux—in the energy balance in Eq. (2) gives

Eq. (7)

$$\rho c_p\frac{\partial ({}^{0}T+\zeta {}^{1}T)}{\partial t}+{\zeta }^{-2}\frac{{\partial }^{-1}{q}_i}{\partial y_i}+{\zeta }^{-1}(\frac{{\partial }^{-1}{q}_i}{\partial x_i}+\frac{{\partial }^0{q}_i}{\partial y_i})+{\zeta }^0(\frac{{\partial }^0{q}_i}{\partial x_i}+\frac{{\partial }^1{q}_i}{\partial y_i})+O(\zeta )+\cdots ={}^{\zeta }Q.$$

Note the truncation at $O(\zeta )$. The terms are now rearranged and collected in orders of $\zeta$. More rigorously, Eq. (7) is multiplied by ${\zeta }^2$ first, then by ${\zeta }^1$, which, after taking the limit $\zeta \to 0^{+}$, yields a system of three equations:

Eq. (8)

$$O({\zeta }^{-2}):\frac{{\partial }^{-1}{q}_i}{\partial y_i}=0,$$

Eq. (9)

$$O({\zeta }^{-1}):\frac{{\partial }^{-1}{q}_i}{\partial x_i}+\frac{{\partial }^0{q}_i}{\partial y_i}=0,$$

Eq. (10)

$$O({\zeta }^0):\rho c_p\frac{\partial ({}^{0}T+\zeta {}^{1}T)}{\partial t}+\frac{{\partial }^0{q}_i}{\partial x_i}+\frac{{\partial }^1{q}_i}{\partial y_i}={}^{\zeta }Q.$$

Henceforth, it is assumed that the heat source ${}^{\zeta }Q$ and the boundary conditions ${}^{\zeta }q_N$ and ${}^{\zeta }T_D$ in Eq. (2) depend on the coarse scale only. The system still needs closure.

Multiplying Eq. (8) with ${}^{0}T$ and integrating over the unit cell gives

Eq. (11)

$${\int }_{\mathrm{\Theta }}{}^{0}T\frac{{\partial }^{-1}{q}_i}{\partial y_i}\mathrm{d}y=0.$$

which, after applying the product rule of the divergence operator and the divergence theorem, reads

Eq. (12)

$${\int }_{\partial \mathrm{\Theta }}{}^0{T}^{-1}{q}_i{n}_i\mathrm{d}y-{\int }_{\mathrm{\Theta }}\frac{{\partial }^{0}T}{\partial y_i}{}^{-1}{q}_i\mathrm{d}y=0.$$

Because periodic boundary conditions are employed on the unit cell, the temperature and heat flux on opposite points of the boundary are equal, but the normal points in exactly opposite directions, and therefore, the boundary integral vanishes. Inserting the corresponding coefficient from the ansatz in Fourier’s law Eq. (5) gives

Eq. (13)

$${\int }_{\mathrm{\Theta }}\frac{{\partial }^{0}T}{\partial y_i}\lambda \frac{{\partial }^{0}T}{\partial y_i}\mathrm{d}y=0.$$

From $\lambda >0$ follows $\frac{{\partial }^{0}T}{\partial y_i}=0$, which implies an independence on the fine scale, i.e., ${}^{0}T≔{}^{0}T(\mathbf{x})$.

Through Eq. (5), this also holds for the corresponding part of the heat flux

Eq. (14)

$${}^{-1}q_i=-\lambda \frac{\partial {}^{0}T}{\partial y_i}=0.$$

Now, Eq. (9) is considered. Inserting Eq. (14) results in

Eq. (15)

$$\frac{{\partial }^0{q}_i}{\partial y_i}=0.$$

Again, the expression for the corresponding part of the heat flux from Eq. (5) can be inserted, which yields

Eq. (16)

$$\frac{\partial }{\partial y_i}\left(\lambda \right(\frac{{\partial }^{0}T}{\partial x_i}+\frac{{\partial }^{1}T}{\partial y_i}\left)\right)=0.$$

Using the ansatz ${}^{1}T(\mathbf{x},\mathbf{y})=H_j(\mathbf{y}){\partial }_{x_j}{}^{0}T(\mathbf{x})$ and factoring out ${\partial }_{x_j}{}^{0}T(\mathbf{x})$ gives

Eq. (17)

$$\frac{\partial }{\partial y_i}\left(\lambda \right({\delta }_{ij}+\frac{\partial H_j(\mathbf{y})}{\partial y_i}\left)\right)=0,$$

with the temperature influence function $H_j:{}^{\zeta }\mathrm{\Omega }\to \mathbb{R}\in C^0(\mathrm{\Theta })$ and the usual Dirac delta ${\delta }_{ij}$. Boundary conditions depend on the coarse scale only, and the fact that ${}^{0}T$ already accounts for all information on the boundary, the fine scale problem reads

Find $H_j(\mathbf{y})$, such that

Eq. (18)

$$\begin{gathered} \frac{\partial }{\partial y_i}\left(\lambda \right({\delta }_{ij}+\frac{\partial H_j(\mathbf{y})}{\partial y_i}\left)\right)=0,\text{in}\mathrm{\Theta } \\ {H}_j(\mathbf{y})=0\text{on}\partial \mathrm{\Theta }. \end{gathered}$$

With the current information, the temperature ansatz Eq. (3) is given as

Eq. (19)

$${}^{\zeta }T(\mathbf{x})=T(\mathbf{x},\mathbf{y})\approx {}^{0}T(\mathbf{x})+{\zeta }^1{H}_j(\mathbf{y})\frac{{\partial }^{0}T(\mathbf{x})}{\partial x_j}.$$

On the coarse scale, the temperature is defined to be the multiscale temperature averaged over the unit cell:

Eq. (20)

$$\begin{gathered} {}^{\mathrm{c}}T(\mathbf{x})=\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }T(\mathbf{x})\mathrm{d}y\approx \frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{0}T(\mathbf{x})+{\zeta }^1{H}_j(\mathbf{y})\frac{{\partial }^{0}T(\mathbf{x})}{\partial x_j}\mathrm{d}y \\ =\frac{{}^{0}T(\mathbf{x})}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\mathrm{d}y+\frac{{\zeta }^1{\partial }_{x_j}{}^{0}T(\mathbf{x})}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{H}_j(\mathbf{y})\mathrm{d}y, \end{gathered}$$

and thus, ${}^{\mathrm{c}}T={}^{0}T$, if and only if ${\int }_{\mathrm{\Theta }}{H}_j(\mathbf{y})\mathrm{d}y=0$.

Inserting the corresponding relation of the heat flux ansatz Eq. (5) into Eq. (17) yields

Eq. (21)

$$\begin{gathered} {}^{0}q_i=-{\mathrm{\Lambda }}_{ij}\frac{{\partial }^{0}T}{\partial x_j} \\ {\mathrm{\Lambda }}_{ij}≔\lambda ({\delta }_{ij}+\frac{\partial H_j}{\partial y_i}), \end{gathered}$$

in which ${\mathrm{\Lambda }}_{ij}$ is termed the heat flux influence function.

Finally, the zeroth-order terms of $\zeta$ Eq. (10) are averaged over the unit cell, resulting in

Eq. (22)

$$\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\rho c_p\frac{{\partial }^{0}T}{\partial t}+\frac{{\partial }^0{q}_i}{\partial x_i}\mathrm{d}y+\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\frac{{\partial }^1{q}_i}{\partial y_i}\mathrm{d}y-\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }Q\mathrm{d}y=0.$$

Rearranging the first integral and applying the divergence theorem to the second gives

Eq. (23)

$$\frac{{\partial }^{0}T}{\partial t}\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\rho c_p\mathrm{d}y+\frac{\partial }{\partial x_i}\left(\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{0}q_i\mathrm{d}y\right)+\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\partial \mathrm{\Theta }}{}^{1}q_i{n}_i\mathrm{d}y-\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }Q\mathrm{d}y=0,$$

wherein the boundary integral vanishes due to ${}^{1}q_i=-\lambda {\partial }_{x_i}$ ${}^{1}T=-\lambda H_j{\partial }_{x_j}$ ${}^{0}T$ and $H_j=0$ for $y\in \partial \mathrm{\Theta }$, resulting in

Eq. (24)

$$\frac{{\partial }^{0}T}{\partial t}\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\rho c_p\mathrm{d}y+\frac{\partial }{\partial x_i}\left(\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{0}q_i\mathrm{d}y\right)-\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }Q\mathrm{d}y=0.$$

Henceforth, averaged quantities are denoted with a horizontal line above them, e.g., $\overline{\rho c_p}≔\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\rho c_p\mathrm{d}y$. With this notation in place, the coarse scale heat conduction reads

Eq. (25)

$$\frac{{\partial }^{0}T}{\partial t}\overline{\rho c_p}+\frac{\partial \overline{{{}^{0}q}_i}}{\partial x_i}-\overline{{}^{\zeta }Q}=0.$$

In summary, two problems have to be solved: one on the fine scale and one on the coarse scale. The fine scale problem reads

Find $H_j(\mathbf{y})$, such that

Eq. (26)

$$\begin{gathered} \frac{\partial }{\partial y_i}\left(\lambda \right({\delta }_{ij}+\frac{\partial H_j}{\partial y_i}\left)\right)=0,\text{in}\mathrm{\Theta } \\ {H}_j=0,\text{on}\partial \mathrm{\Theta }. \end{gathered}$$

From the solution, the heat flux influence function ${\mathrm{\Lambda }}_{ij}$ can be computed using Eq. (21). The coarse scale problem is then defined as

Find ${}^{0}T(\mathbf{x})$, such that

Eq. (27)

$$\begin{gathered} \overline{\rho c_p}\frac{{\partial }^{0}T}{\partial t}-\frac{\partial }{\partial x_i}(\overline{{\mathrm{\Lambda }}_{ij}}\frac{{\partial }^{0}T}{\partial x_j})=\overline{Q},\text{in}\mathrm{\Omega }\times [0,t_{\max }] \\ {}^{0}T={\overline{T}}_0,\text{in}\mathrm{\Omega }\times \{0\} \\ \overline{{\mathrm{\Lambda }}_{ij}}\frac{{\partial }^{0}T}{\partial x_j}{n}_i=q_{\mathrm{N}},\text{on}\partial \mathrm{\Omega }\begin{array}{c}.\end{array} \end{gathered}$$

The overall temperature function ${}^{\zeta }T$ can then be constructed as

Eq. (28)

$${}^{\zeta }T(\mathbf{x},\mathbf{y})\approx {}^{0}T+\zeta H_j(\mathbf{y})\frac{{\partial }^{0}T(\mathbf{x})}{\partial x_j}.$$

2.3.

Thermo-Elastic Deformations

To define the multiscale structural-mechanics equations, which describe the deformation of a solid under volume and thermal loading conditions, a few functions have to be introduced. For an $n$-dimensional domain, the displacement is denoted $u:{}^{\zeta }\mathrm{\Omega }\to {\mathbb{R}}^n$. Note that it does not depend on time. The body forces acting on the work piece are written as $b:{}^{\zeta }\mathrm{\Omega }\to {\mathbb{R}}^n$. The notations for strains and stresses resulting from the loads read $ϵ:{}^{\zeta }\mathrm{\Omega }\to {\mathbb{R}}^{n\times n}$ and $\sigma :{}^{\zeta }\mathrm{\Omega }\to {\mathbb{R}}^{n\times n}$, respectively. They are related through the stiffness tensor $C\in {\mathbb{R}}^{n\times n\times n\times n}$. Using these definitions, the static equilibrium of stresses in index notation is given as

Eq. (29)

$$\begin{gathered} \frac{\partial {}^{\zeta }\sigma _{ij}}{\partial x_j}+{}^{\zeta }b_i=0,\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }\sigma _{ij}={}^{\zeta }C_{ijkl}{}^{\zeta }ϵ_{kl}^{\mathrm{el}},\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }\sigma _{ij}{n}_j={}^{\zeta }\overline{t}_i,\text{on}{}^{\zeta }\mathrm{\Gamma }_N \\ {}^{\zeta }u_i={}^{\zeta }\overline{u}_i,\text{on}{}^{\zeta }\mathrm{\Gamma }_D. \end{gathered}$$

Herein, the total ${ϵ}^{\mathrm{tot}}$ strain is the sum of the elastic ${ϵ}^{\mathrm{el}}$, plastic ${ϵ}^{\mathrm{pl}}$, and thermal ${ϵ}^{\mathrm{th}}$ strains:

Eq. (30)

$$\begin{gathered} {}^{\zeta }ϵ_{kl}^{\mathrm{tot}}={}^{\zeta }ϵ_{kl}^{\mathrm{el}}+{}^{\zeta }ϵ_{kl}^{\mathrm{th}},\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }ϵ_{kl}^{\mathrm{tot}}=\frac{1}{2}(\frac{\partial {}^{\zeta }u_k}{\partial x_l}+\frac{\partial {}^{\zeta }u_l}{\partial x_k}),\text{in}{}^{\zeta }\mathrm{\Omega } \\ {}^{\zeta }ϵ_{kl}^{\mathrm{th}}=-\alpha \mathrm{\Delta }{}^{\zeta }T{\delta }_{kl},\text{in}{}^{\zeta }\mathrm{\Omega }. \end{gathered}$$

The stiffness tensor obeys a symmetry condition

Eq. (31)

$${}^{\zeta }C_{ijkl}={}^{\zeta }C_{jikl}={}^{\zeta }C_{ijlk}={}^{\zeta }C_{klij},$$

and is positive in the sense that

Eq. (32)

$$\exists c>0:{}^{\zeta }C_{ijkl}{n}_{ij}{n}_{kl}\ge c{n}_{ij}{n}_{ij},\forall n_{ij}=n_{ji}.$$

Analogously to the two-scale heat conduction in Sec. 2.2, the asymptotic ansatz for the displacement is defined to be

Eq. (33)

$${}^{\zeta }u_i(\mathbf{x})≔u_i(\mathbf{x},\mathbf{y})={}^{0}u_i(\mathbf{x},\mathbf{y})+{\zeta }^1{u}_i(\mathbf{x},\mathbf{y})+{\zeta }^2{{}^{2}u}_i(\mathbf{x},\mathbf{y})+O({\zeta }^3)+\cdots .$$

From Eq. (29), it is obvious that the derivative is required, too. Taking the derivative with respect to $x_j$, applying the chain rule, and sorting in ascending order of $\zeta$ yields

Eq. (34)

$$\begin{gathered} \frac{\partial u_i}{\partial x_j}(\mathbf{x},\mathbf{y})={\zeta }^{-1}\frac{{\partial }^0{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y})+{\zeta }^0(\frac{{\partial }^0{u}_i}{\partial x_j}(\mathbf{x},\mathbf{y})+\frac{{\partial }^1{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y}))+{\zeta }^1(\frac{{\partial }^1{u}_i}{\partial x_j}(\mathbf{x},\mathbf{y})+\frac{{\partial }^2{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y})) \\ +{\zeta }^2(\frac{{\partial }^2{u}_i}{\partial x_j}(\mathbf{x},\mathbf{y})+\frac{{\partial }^3{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y}))+O({\zeta }^3)+\cdots \\ ={\zeta }^{-1}\frac{{\partial }^0{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y})+\sum _{s=0}^2{\zeta }^s(\frac{{\partial }^s{u}_i}{\partial x_j}(\mathbf{x},\mathbf{y})+\frac{{\partial }^{s+1}{u}_i}{\partial y_j}(\mathbf{x},\mathbf{y}))+O({\zeta }^3)+\cdots . \end{gathered}$$

Inserting into the strain–displacement relation in Eq. (30) gives

Eq. (35)

$$\begin{gathered} {ϵ}_{ij}^{\mathrm{el}}(\mathbf{x},\mathbf{y})={\zeta }^{-1}\frac{1}{2}(\frac{{\partial }^0{u}_i}{\partial y_j}+\frac{{\partial }^0{u}_j}{\partial y_i})+\sum _{s=0}^2{\zeta }^s\frac{1}{2}(\frac{{\partial }^s{u}_i}{\partial x_j}+\frac{{\partial }^s{u}_i}{\partial x_j}+\frac{{\partial }^{s+1}{u}_i}{\partial y_j}+\frac{{\partial }^{s+1}{u}_i}{\partial y_j})+O({\zeta }^3)+\cdots \\ ={\zeta }^{-1-1}{ϵ}_{ij}^{\mathrm{el}}+\sum _{s=0}^2{\zeta }^s{}^sϵ_{ij}^{\mathrm{el}}+O({\zeta }^3)+\cdots . \end{gathered}$$

Assuming the stiffness tensor is independent of the coarse scale, i.e., $C_{ijkl}(\mathbf{x},\mathbf{y})≔C_{ijkl}(\mathbf{y})$, and inserting the ansatz for the elastic strain Eq. (35) into the stress–strain relation Eq. (29) results in the two-scale ansatz for the stress

Eq. (36)

$$\begin{gathered} {\sigma }_{ij}(\mathbf{x},\mathbf{y})=C_{ijkl}(\mathbf{y})({\zeta }^{-1-1}{\epsilon }_{ij}^{\mathrm{el}}+\sum _{s=0}^2{\zeta }^s{}^s\epsilon _{ij}^{\mathrm{el}})+O({\zeta }^3)+\cdots \\ ={\zeta }^{-1-1}{\sigma }_{ij}+\sum _{s=0}^2{\zeta }^s{}^s\sigma _{ij}+O({\zeta }^3)+\cdots . \end{gathered}$$

Again, the derivative of the stress with respect to space is needed as it appears in Eq. (29). Thus, inserting the derivative, which reads

Eq. (37)

$$\frac{\partial {\sigma }_{ij}}{\partial x_j}(\mathbf{x},\mathbf{y})={\zeta }^{-2}\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial y_j}+{\zeta }^{-1}(\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^0{\sigma }_{ij}}{\partial y_j})+{\zeta }^0(\frac{{\partial }^0{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^1{\sigma }_{ij}}{\partial y_j})+O(\zeta ),$$

in the stress balance and sorting in ascending order of $\zeta$ yields

Eq. (38)

$${\zeta }^{-2}\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial y_j}+{\zeta }^{-1}(\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^0{\sigma }_{ij}}{\partial y_j})+{\zeta }^0(\frac{{\partial }^0{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^1{\sigma }_{ij}}{\partial y_j}+{}^{\zeta }b_i)+O(\zeta )=0.$$

Analogously to the derivation of the two-scale heat conduction problem, the leading order terms of the stress balance are obtained by multiplication with ${\zeta }^2$ and $\zeta$, respectively, and taking the limit $\zeta \to 0$:

Eq. (39)

$$O({\zeta }^{-2}):\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial y_j}=0,$$

Eq. (40)

$$O({\zeta }^{-1}):\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^0{\sigma }_{ij}}{\partial y_j}=0,$$

Eq. (41)

$$O({\zeta }^0):\frac{{\partial }^0{\sigma }_{ij}}{\partial x_j}+\frac{{\partial }^1{\sigma }_{ij}}{\partial y_j}+{}^{\zeta }b_i=0.$$

It is, hereafter, assumed that the body force depends on the fine scale only ${}^{\zeta }b_i(\mathbf{x})≔b_i(\mathbf{y})$, and the boundary conditions in Eq. (29) are coarse scale functions ${}^{\zeta }\overline{t}_i(\mathbf{x})≔{\overline{t}}_i(\mathbf{x})$ and ${}^{\zeta }\overline{u}_i(\mathbf{x})≔{\overline{u}}_i(\mathbf{x})$. For closure, the moments are computed for each order. Therefore, multiplying Eq. (39) by ${}^{0}u_i$ and integrating over the unit cell gives

Eq. (42)

$${\int }_{\mathrm{\Theta }}{}^{0}u_i\frac{{\partial }^{-1}{\sigma }_{ij}}{\partial y_j}\mathrm{d}y=0.$$

Using the product rule of divergence and the divergence theorem, as well as noting that the resulting boundary integral vanishes due to periodicity of both integrands and opposing normals, results in

Eq. (43)

$${\int }_{\mathrm{\Theta }}\frac{{\partial }^0{u}_i}{\partial y_j}{}^{-1}{\sigma }_{ij}\mathrm{d}y=0.$$

Inserting the definitions of ${{}^{-1}\sigma }_{ij}$ and ${{}^{-1}ϵ}_{ij}^{\mathrm{el}}$ from Eqs. (36) and (35), respectively, yields

Eq. (44)

$${\int }_{\mathrm{\Theta }}\frac{{\partial }^0{u}_i}{\partial y_j}{C}_{ijkl}\frac{1}{2}(\frac{{\partial }^0{u}_k}{\partial x_l}+\frac{{\partial }^0{u}_l}{\partial x_k})\mathrm{d}y=0.$$

From the symmetry and positiveness of the stiffness tensor $C$, it can be inferred that

Eq. (45)

$$\frac{{\partial }^0{u}_i}{\partial x_j}=0,$$

and, with Eqs. (36) and (35),

Eq. (46)

$${}^{-1}{\sigma }_{ij}=0.$$

Inserting Eq. (46) into Eq. (40) gives

Eq. (47)

$$\frac{\partial {{}^0\sigma }_{ij}}{\partial y_j}=0,$$

which, with the definition of ${}^0\sigma _{ij}$ from Eq. (36) and ${}^0{ϵ}_{\mathrm{ij}}$ from Eq. (35), results in

Eq. (48)

$$\frac{\partial }{\partial y_j}\left(C_{ijkl}\right(\frac{1}{2}(\frac{{\partial }^0{u}_k}{\partial x_l}+\frac{{\partial }^0{u}_l}{\partial x_k})+\frac{1}{2}(\frac{{\partial }^1{u}_k}{\partial x_l}+\frac{{\partial }^1{u}_l}{\partial x_k})\left)\right)=0.$$

This shows a direct dependence between the displacement of the zeroth and first scale. Hence, the separation of variables ansatz

Eq. (49)

$${{}^{1}u}_i≔{}^{mn}H_i(\mathbf{y})\frac{1}{2}(\frac{{\partial }^0{u}_m}{\partial x_n}+\frac{{\partial }^0{u}_n}{\partial x_m}),$$

with the first-order displacement influence function $H$, is employed. It is assumed that the influence function is symmetric ${}^{mn}H_i={}^{mn}H_i$, locally periodic, and continuous $H\in C^0(\mathrm{\Omega })$. Inserting the ansatz in Eq. (48) yields

Eq. (50)

$$\frac{\partial }{\partial y_j}(C_{ijkl}({}^{0}u_{(k,x_l)}+{}^{mn}H_{(k,y_l)}^0{u}_{(m,x_n)}))=0,$$

and after factoring out ${}^{0}u_{(m,x_n)}$

Eq. (51)

$$\frac{\partial }{\partial y_j}(C_{ijkl}(I_{klmn}+{}^{mn}H_{(k,y_l)})){}^{0}u_{(m,x_n)}=0,$$

with the definition $I_{klmn}≔\frac{1}{2}({\delta }_{km}{\delta }_{ln}+{\delta }_{lm}{\delta }_{kn})$. Because ${{}^{0}u}_{(m,x_n)}$ is arbitrary, it is required that

Eq. (52)

$$\frac{\partial }{\partial y_j}(C_{ijkl}(I_{klmn}+{}^{mn}H_{(k,y_l)}))=0.$$

To have a well-defined problem in the sense of Hadamard, an additional condition is needed. In the literature,^6–8 two conditions are reported:

• 1. ${}^{mn}H_k(\mathbf{y})=0,\text{on}\partial {\mathrm{\Theta }}_{\mathrm{vert}}$, and

• 2. ${\int }_{\mathrm{\Theta }}{}^{mn}H_k(\mathbf{y})\mathrm{d}y=0,\text{in}\mathrm{\Theta }$.

Herein, $\partial {\mathrm{\Theta }}_{\mathrm{vert}}$ denotes the vertices of the unit cell boundary $\partial \mathrm{\Theta }$. As always, each has advantages over the other. Although condition 1 is simpler to implement, condition 2 associates ${}^{0}u_i(\mathbf{x})$ with the average displacement ${}^{\mathrm{c}}{u}_i$ through

Eq. (53)

$${}^{\mathrm{c}}u(x):=\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }u(x,y)\mathrm{d}y.$$

With the fine-scale displacement ansatz, the overall displacement function reads

Eq. (54)

$$u_i(\mathbf{x},\mathbf{y})={}^{0}u_i(\mathbf{x})+{\zeta }^{mn}{H}_i(\mathbf{y}){}^{0}u_{(m,x_n)}(\mathbf{x})+O({\zeta }^2),$$

which can be inserted into Eq. (53)

Eq. (55)

$$\begin{gathered} {}^{\mathrm{c}}u_i(\mathbf{x})≔\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{0}u_i(\mathbf{x})+{\zeta }^{mn}{H}_i(\mathbf{y}){{}^{0}u}_{(m,x_n)}(\mathbf{x})+O({\zeta }^2)\mathrm{d}y \\ ={}^{0}u_{(m,x_n)}(\mathbf{x})\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}1\mathrm{d}y \\ ={}^{0}u(\mathbf{x}). \end{gathered}$$

Within this treatise, condition 2 has been applied. From this, the leading order strain is given as

Eq. (56)

$$\begin{gathered} {}^0ϵ_{ij}^{\mathit{tot}}={}^{0}u_{(i,x_j)}+{}^{1}u_{(i,x_j)} \\ =(I_{ijkl}+{}^{mn}H_{(k,y_l)}(y)){}^{0}u_{(m,x_n)}=E_{mnkl}{}^{0}u_{(m,x_n)} \\ =:{}^fϵ_{\mathrm{kl}}^{\mathit{tot}}(x,y), \end{gathered}$$

with the strain influence function $E_{mnkl}≔I_{ijkl}+{}^{mn}H_{(k,y_l)}(y)$ and the fine-scale total strain ${}^fϵ^{\mathrm{tot}}$. Averaging the fine-scale total strain over the unit cell yields the coarse-scale total strain

Eq. (57)

$$\begin{gathered} {}^cϵ_{mn}^{\mathit{tot}}(x)≔\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^fϵ_{mn}^{\mathit{tot}}(x,y)\mathrm{d}y \\ ={}^{0}u_{(m,x_n)}. \end{gathered}$$

Herein, the integral vanishes due to the divergence identity, the divergence theorem, and the local periodicity of the displacement influence function. With this, the coarse scale total strain is given as

Eq. (58)

$${}^fϵ_{kl}^{\mathit{tot}}(x,y)≔E_{mnkl}(y){}^cϵ_{\mathrm{mn}}^{\mathit{tot}}(x).$$

The leading stress is computed according to

Eq. (59)

$${}^0\sigma _{ij}(x,y)≔{\mathrm{\Sigma }}_{ijmn}(y){}^cϵ_{mn}^{\mathrm{tot}}(x)=:{}^f\sigma _{ij}(x,y),$$

with the stress influence function ${\mathrm{\Sigma }}_{ijmn}≔C_{ijkl}{E}_{klmn}$.

Finally, the highest order, i.e., Eq. (41), averaged over the unit cell reads

Eq. (60)

$$\begin{gathered} \frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\frac{{\partial }^0{\sigma }_{ij}}{\partial y_j}+\frac{{\partial }^1{u}_{ij}}{\partial y_j}+{}^{\zeta }b_i\mathrm{d}y=0 \\ \frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}\frac{{\partial }^0{\sigma }_{ij}}{\partial y_j}+{}^{\zeta }b_i\mathrm{d}y=0. \end{gathered}$$

Herein, the second term vanishes after applying the divergence theorem and due to local periodicity. Inserting the stress–strain and strain–displacement relations for the respective scale yields

Eq. (61)

$$\begin{gathered} \frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{C}_{ijkl}{E}_{klmn}(\mathbf{y})\mathrm{d}y\frac{\partial }{\partial x_j}{{}^{0}u}_{(m,x_n)}(\mathbf{x}))+\frac{1}{\mathrm{vol}(\mathrm{\Theta })}{\int }_{\mathrm{\Theta }}{}^{\zeta }b_i\mathrm{d}y=0 \\ {{}^{\mathrm{c}}C}_{ijmn}\frac{\partial }{\partial x_j}({}^cϵ_{mn}^{\mathrm{tot}})+{}^{\mathrm{c}}b_i=0 \\ \frac{{\partial }^{\mathrm{c}}{\sigma }_{ij}}{\partial x_j}+{}^{\mathrm{c}}b_i=0, \end{gathered}$$

with the according definitions of ${{}^{\mathrm{c}}C}_{ijmn}$ and ${{}^{\mathrm{c}}\sigma }_{ij}$.

In summary, on the fine scale, the problem reads

Find ${}^{mn}H_i(\mathbf{y})$, such that

Eq. (62)

$$\begin{gathered} \frac{\partial }{\partial y_j}(C_{ijkl}({}^{mn}H_{(k,y_l)}+I_{klmn}))=0,\text{in}\mathrm{\Theta } \\ {\int }_{\mathrm{\Theta }}{}^{mn}H_i(y)\mathrm{d}y=0,\text{in}\mathrm{\Theta } \\ {}^{mn}H_i(y)={}^{mn}H_i(y+l),\text{on}\partial \mathrm{\Theta }, \end{gathered}$$

where $l$ is chosen, such that $y+l$ is the point on the boundary opposite of $y$. The problem on the coarse scale is

Find ${{}^{\mathrm{c}}u}_i(\mathbf{x})$, such that

Eq. (63)

$$\begin{gathered} \frac{\partial {}^c\sigma _{ij}}{\partial x_j}+{}^{c}b_i=0,\text{in}\mathrm{\Omega } \\ {}^c\sigma _{ij}={}^{c}C_{ijmn}{}^cϵ_{\mathrm{mn}}^{\mathrm{tot}},\text{in}\mathrm{\Omega } \\ {}^{c}u_i={}^c\overline{u}_i,\text{on}{\mathrm{\Gamma }}_D \\ {}^c\sigma _{ij}{n}_j={}^c\overline{t}_i,\text{on}{\mathrm{\Gamma }}_N. \end{gathered}$$

Finally, the overall displacements can be recovered with

Eq. (64)

$$u_i(\mathbf{x},\mathbf{y})={}^{0}u_i(\mathbf{x})+{\zeta }^{mn}{H}_k(\mathbf{y}){}^{0}u_{(m,x_n)}(\mathbf{x}).$$