Multiscale model for ultrashort pulsed parallel laser structuring—part II. The macroscale model

Abstract. The increasing pixel density in displays demands high quality in the production of fine metal masks (FMMs). The production process of FMMs boils down to structuring tiny holes in thin metal sheets or foils. The manufacturing requirements of FMMs are high precision in terms of the hole geometry to let enough light escape from each diode and high productivity to produce the required amount. To achieve both objectives, high power ultrashort pulse (USP) lasers can be utilized. Because USP lasers fall short of the productivity requirements, they are combined with multibeam scanners. During production, the multibeam scanners deposit a lot of heat in the metal foil, which can ultimately yield temperature-induced distortions. To understand and finally avoid such distortions, a process simulation is sought. In a preceding study, the structuring of a single hole (the microscale) was investigated, but due to the large differences in the time and spatial scales involved, it was not feasible to simulate the production of the whole part (the macroscale). Within this treatise, a multiscale approach that takes into account the necessary information from the microscale to describe temperature-induced distortions on the macroscale is described. This approach targets laser ablation processes with pulse durations ranging from picoseconds up to nanoseconds provided the ablation is not melt-driven. First, a representative volume element (RVE) is generated from the results of the microscale model. Then, this RVE is utilized in the thermo-elastic structural mechanics simulation on the macroscale. The multiscale model is validated numerically against a hole-resolved computation, which shows good agreement. Naturally, the simulation is highly dependent on the microscale model, which in turn depends on the material properties. To handle material changes well, an experimental calibration has to be performed. This calibration is not part of this treatise, but it will be described in a future publication. In addition to the calibration process, the validation with experiments will be conducted in future research. Additionally, the authors envision the automation of the whole process, resulting in a first-time-right approach for the development of FMMs. Finally, the procedure might be extended to the requirements of other filtration purposes.


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. 1 A technology that allows for such quality is high power ultrashort pulsed lasers. 2 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. 2 To meet the productivity requirements, the laser can be combined with a multibeam scanner. 3 During the production of FMMs, it is important to avoid *Address all correspondence to Christian Heinigk, christian.heinigk@nld.rwth-aachen.de distortions or discoloration due to heat accumulation. 4 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, 5 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][7][8][9][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. 11 The unit cell, which is the computational domain in the microscale simulation, is the direct result of the simulation described in the prior publication. 5 Both the heat conduction and the thermo-elastic structural mechanics problems are discretized using a Bubnov-Galerkin method. [12][13][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][13][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.

Multiscale Model
In a prior publication, the ablation of a structure with a single beam was explained. 5 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.

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 < ζ ≔ l f l c ≪ 1. Variables referring to a particular length scale are denoted with a subscript x f for fine and x c for coarse. Furthermore, the notation for functions that account for information on all scales uses a preceding superscript ζ 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∕ζ lives on the fine scale. The computational domain is denoted as Ω with the boundary ∂Ω. To account for different physical interactions with the environment, it is necessary to split this boundary into a Neumann part Γ N and a Dirichlet part Γ D , such that ∂Ω ¼ Γ N ∪ Γ D and Γ N ∩ Γ D ¼ 0 (c.f. Fig. 1). Analogously, the domain and boundaries for the unit cell are defined, except the symbol Θ is used. Whenever a distinction between Dirichlet and Neumann Boundary has to be made, Γ D and Γ N , respectively, are used.
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, 5 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 ∂ x i f to emphasize vector field operators such as the divergence. Averaged quantities are denoted with an overline, e.g., x, vectors and matrices are highlighted in a boldface x and uppercase boldface X, respectively. Finally, the subscript f ði;x j Þ is used to denote the symmetrized gradient tensor entry, i.e. E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 4 5 8 More complex geometries can be created using a level set function ΦðxÞ∶R n → R. Such a level set function is the key result in the companion publication 5 wherein Φ 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 Φð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.

Two-Scale Heat Conduction
The heat conduction equation for a material with specific density ρ∶ ζ Ω → R, specific heat capacity c p ∶ ζ Ω → R, and thermal conductivity λ∶ ζ Ω → R heated by a heat source Q∶ ζ Ω → R reads  Fig. 2(a)]. The final unit cell after removal of marked elements is shown in Fig. 2 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.
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∶ ζ Ω → 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., ρ ≔ ρðyÞ and c p ≔ c p ðyÞ. Now, an ansatz function for the temperature is defined by means of an asymptotic expansion 15 using the scale ζ, which reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 6 4 ζ TðxÞ ≔ Tðx; yÞ ¼ ζ 00 Tðx; yÞ þ ζ 11 Tðx; yÞ þ Oðζ 2 Þþ · · · : (3) Note that terms of order 2 and above have been omitted as they do not contribute significantly to the temperature due to the condition ζ ≪ 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 ζ yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 4 8 5 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 yÞ − Oðζ 2 Þ− · · · q i ðxÞ ¼ ζ −1−1 q i ðx; yÞ þ ζ 00 q i ðx; yÞ þ ζ 11 q i ðx; yÞ þ Oðζ 2 Þþ · · · : 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 ζ yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 3 1 9 Here, the functions' arguments are dropped for brevity. Inserting the ansatz functions-Eq.
Note the truncation at OðζÞ. The terms are now rearranged and collected in orders of ζ. More rigorously, Eq. (7) is multiplied by ζ 2 first, then by ζ 1 , which, after taking the limit ζ → 0 þ , yields a system of three equations: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 1 4 3 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 8 6 Oðζ −1 Þ∶ Henceforth, it is assumed that the heat source ζ Q and the boundary conditions ζ q N and ζ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 6 2 which, after applying the product rule of the divergence operator and the divergence theorem, reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 5 9 3 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 5 0 0 From λ > 0 follows ∂ 0 T ∂y i ¼ 0, which implies an independence on the fine scale, i.e., 0 T ≔ 0 TðxÞ. Through Eq. (5), this also holds for the corresponding part of the heat flux E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 4 2 7 Now, Eq. (9) is considered. Inserting Eq. (14) results in E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 3 6 9 ∂ 0 q i ∂y i ¼ 0: (15) Again, the expression for the corresponding part of the heat flux from Eq. (5) can be inserted, which yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 3 0 0 Using the ansatz 1 Tðx; yÞ ¼ H j ðyÞ∂ x j 0 TðxÞ and factoring out ∂ x j 0 TðxÞ gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 2 4 0 with the temperature influence function H j ∶ ζ Ω → R ∈ C 0 ðΘÞ and the usual Dirac 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 ðyÞ, such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 1 4 4 With the current information, the temperature ansatz Eq. (3) is given as On the coarse scale, the temperature is defined to be the multiscale temperature averaged over the unit cell: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 6 7 4 c TðxÞ ¼ 1 volðΘÞ and thus, Inserting the corresponding relation of the heat flux ansatz Eq. (5) into Eq. (17) yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 6 ; 5 7 2 in which Λ ij is termed the heat flux influence function. Finally, the zeroth-order terms of ζ Eq. (10) are averaged over the unit cell, resulting in E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 6 ; 4 7 2 1 volðΘÞ Rearranging the first integral and applying the divergence theorem to the second gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 6 ; 4 1 4 wherein the boundary integral vanishes due to 1 ; t e m p : i n t r a l i n k -; e 0 2 4 ; 1 1 6 ; 3 2 9 Henceforth, averaged quantities are denoted with a horizontal line above them, e.g., ρc p ≔ 1 volðΘÞ ∫ Θ ρc p dy. With this notation in place, the coarse scale heat conduction reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 1 1 6 ; 2 5 6 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 ðyÞ, such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 1 1 6 ; 1 7 2 From the solution, the heat flux influence function Λ ij can be computed using Eq. (21). The coarse scale problem is then defined as Find 0 TðxÞ, such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 1 1 6 ; 7 2 3 The overall temperature function ζ T can then be constructed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 1 1 6 ; 6 1 8

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∶ ζ Ω → R n . Note that it does not depend on time. The body forces acting on the work piece are written as b∶ ζ Ω → R n . The notations for strains and stresses resulting from the loads read ϵ∶ ζ Ω → R n×n and σ∶ ζ Ω → R n×n , respectively. They are related through the stiffness tensor C ∈ R n×n×n×n . Using these definitions, the static equilibrium of stresses in index notation is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 1 1 6 ; 4 5 7 Herein, the total ϵ tot strain is the sum of the elastic ϵ el , plastic ϵ pl , and thermal ϵ th strains: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 0 ; 1 1 6 ; 3 4 4 The stiffness tensor obeys a symmetry condition E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 1 ; 1 1 6 ; 2 5 3 and is positive in the sense that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 2 ; 1 1 6 ; 2 0 7 ∃ c > 0∶ ζ C ijkl n ij n kl ≥ cn ij n ij ; ∀ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 3 ; 1 1 6 ; 1 4 9 ζ u i ðxÞ ≔ u i ðx; yÞ ¼ 0 u i ðx; yÞ þ ζ 1 u i ðx; yÞ þ ζ 22 u i ðx; yÞ þ Oðζ 3 Þþ · · · : (33) 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 ζ yields Inserting into the strain-displacement relation in Eq. (30) gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 5 ; 1 1 6 ; 6 2 2 Assuming the stiffness tensor is independent of the coarse scale, i.e., C ijkl ðx; yÞ ≔ C ijkl ð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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 6 ; 1 1 6 ; 5 0 1 Again, the derivative of the stress with respect to space is needed as it appears in Eq. (29). Thus, inserting the derivative, which reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 7 ; 1 1 6 ; 3 9 3 in the stress balance and sorting in ascending order of ζ yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 8 ; 1 1 6 ; 3 3 4 Analogously to the derivation of the two-scale heat conduction problem, the leading order terms of the stress balance are obtained by multiplication with ζ 2 and ζ, respectively, and taking the limit ζ → 0: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 9 ; 1 1 6 ; 2 5 1 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 0 ; 1 1 6 ; 1 9 2 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 1 ; 1 1 6 ; 1 5 5 From the symmetry and positiveness of the stiffness tensor C, it can be inferred that ∂ ∂y j C ijkl 1 2 This shows a direct dependence between the displacement of the zeroth and first scale. Hence, the separation of variables ansatz with the definition I klmn ≔ 1 2 ðδ km δ ln þ δ lm δ kn Þ. Because 0 u ðm;x n Þ is arbitrary, it is required that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 5 2 ; 1 1 6 ; 1 0 7 ∂ ∂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 ðyÞ ¼ 0; on ∂Θ vert , and 2. ∫ Θ mn H k ðyÞdy ¼ 0; in Θ.
Herein, ∂Θ vert denotes the vertices of the unit cell boundary ∂Θ. As always, each has advantages over the other. Although condition 1 is simpler to implement, condition 2 associates 0 u i ðxÞ with the average displacement c u i through Within this treatise, condition 2 has been applied. From this, the leading order strain is given as 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 with the stress influence function Σ ijmn ≔ C ijkl E klmn .
Finally, the highest order, i.e., Eq. (41), averaged over the unit cell reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 0 ; 1 1 6 ; 7 2 3 1 volðΘÞ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 1 ; 1 1 6 ; 6 2 2 1 volðΘÞ with the according definitions of c C ijmn and c σ ij . In summary, on the fine scale, the problem reads Find mn H i ðyÞ, such that Finally, the overall displacements can be recovered with

Weak Formulation
To solve the aforementioned systems numerically, a Galerkin method is applied in the space dimension. The time discretization is performed using the Euler method. First, the two-scale heat conduction is explained, and second, the thermo-elasticity problem is discussed.

Two-Scale Heat Conduction
Remembering the equation for the temperature influence function, i.e., Eq. (26), and choosing trial and test functions from the Hilbert space ϕ; φ ∈ H 1 0 ðΩÞ, containing functions vanishing on the boundary, the temperature influence function can be approximated within this space to read H j ðyÞ ≈Ĥ jk ϕ k ðyÞ; k ¼ 1; : : : ; K: (65)

Multiplying Eq. (26) by arbitrary test functions and integrating over the unit cell yields
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 6 ; 1 1 6 ; 6 9 9 Z Θ ∂ ∂y i λ δ ij þ ∂Ĥ jk ϕ k ∂y i φ l dy ¼ 0; ∀ φ l ; l ¼ 1; : : : ; L: (66) Constant terms can be factored out and brought to the right side. Using the product rule for the divergence operator and applying the divergence theorem gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 7 ; 1 1 6 ; 6 2 8 Note that the boundary integral vanishes due to vanishing test functions. With the definition of a system matrix A ∈ R L×K and a thermal load vector b j ∈ R L , such that E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 6 8 ; 1 1 6 ; 5 5 9 The problem can be written in matrix form, i.e.
On the macroscale, the discretization in space is performed analogously. Again, the same notation is used to refer to trial and test functions. Note, however, that the functions are defined on the macroscopic space this time, e.g., ϕ ¼ ϕðxÞ. Representing the temperature in the Hilbert space ∀ φ l ; l ¼ 1; : : : ; L: Henceforth, the arguments of test and trial functions are dropped for a terser notation. The second term can be rewritten using the divergence theorem and integrating the Neumann boundary condition from Eq. (27) and, therefore, it can be brought to the right side, too: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 3 ; 1 1 6 ; 1 7 3 Now, the time domain is discretized using an explicit Euler scheme. Henceforth, the time step is denoted with a superscript n and the time step with Δt. With this notation in place, the temperature coefficients at the next time stepT nþ1 k can be computed as With the definition of the system matrix A ∈ R L×K and the thermal load vector b n ∈ R L , i.e.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 5 ; 1 1 6 ; 6 8 7 the matrix form of Eq. (74) reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 6 ; 1 1 6 ; 6 0 4

Thermo-Elastic Deformations
In analogy to Sec. 3.1, the weak formulation is now derived for the thermo-elastic deformations.
For the fine-scale problem, the test and trial functions are chosen from the same Hilbert space Applying the product identity for the divergence operator followed by the divergence theorem and taking into account the test functions vanishing on the boundary gives E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 7 9 ; 1 1 6 ; 3 7 6 − Z Θ C ijkl ð mnĤq q ϕ ðk;y l Þ þ I klmn Þ ∂φ p ∂y j dy ¼ 0: The known terms can be written on the right side, which reads E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 0 ; 1 1 6 ; 3 1 9 To turn this expression into matrix form, the system matrix A ∈ R P×Q and load vector b ∈ R P are defined. In coefficient notation, they read E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 1 ; 1 1 6 ; 2 4 9 With this, the matrix form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 2 ; 1 1 6 ; 1 6 2 can be solved for the coefficient vector mnĤ ∈ R Q . The solution to the discretized fine-scale problem is used to define the homogenized coarsescale material properties. First, the displacement influence function can be reconstructed using E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 3 ; 1 1 6 ; 9 1 mn H i ðyÞ ≈ mnĤ q ϕ qi ðyÞ: Then, the strain influence function is given as Using, again, the product rule of the divergence operator and the divergence theorem, after inserting the traction defined on the Neumann boundary, gives Inserting the stress-strain and in turn the strain-displacement relations yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 8 9 ; 1 1 6 ; 3 9 5 Z Ω c C ijkl ð cû q q ϕ ðk;x l Þ þ αΔ c Tδ kl Þ Again, constant terms are brought to the right, i.e.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 0 ; 1 1 6 ; 3 4 0 Factoring out the coefficients yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 1 ; 1 1 6 ; 2 8 3 which can be brought into the matrix form E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 9 2 ; 1 1 6 ; 2 2 6

From a Single Beam to Multibeam Patches
To scale up production, the laser beam can be split into multiple beams, each structuring the same cavity in parallel. This scenario can also be accounted for in the described model. The trick is to replace the periodic unit cell domain to account for a multibeam patch. Therefore, instead of the unit cell shown in Fig. 1(a), the unit cell consists of multiple holes ablated from a cuboid (c.f. Fig. 3). Naturally, the volume heat source, which is included in the two-scale heat conduction problem in Sec. 2.2, has to account for the multiple beams as well. In addition, the two-scale formulation is the same.

Numerical Experiments
To validate the two-scale model, the structuring of a foil is simulated with the multiscale approach and compared with a hole-resolved simulation. The geometric setup is similar to the foil shown in Fig. 1(b). The foil's thickness is 2 × 10 −5 m. The holes are assumed to be rectangular (c.f. Fig. 1) Table 1 for the beam parameters).
To structure the hole foil, a hatch strategy is needed. The scanner is set to move in the y-direction first and then in the x-direction with a speed of 0.05 m∕s. For the hole-resolved simulation, each hole is discretized with a mesh of 2 × 2 × 1 elements in the x, y, and z directions, respectively. On the other hand, in the multiscale simulation, the discretization is patch resolved; therefore, each patch corresponds to 1 element.
Beneath the foil, which is made of INVAR36 (c.f. Ref. 16 for the material parameters), an air layer and a soda-lime float glass layer, with a thickness of 5 × 10 −6 m each, are added to dissipate heat and avoid thermally induced plastic distortions which cannot be handled by the described model. The material parameters for air and glass are shown in Tables 2 and 3, respectively.  It remains to define the initial conditions, the boundary conditions, and the sources. The load vectors for the microscale heat conduction problem read b 0 ¼ ð1;0; 0Þ T , b 1 ¼ ð0;1; 0Þ T , and b 2 ¼ ð0;0; 1Þ T . For the mechanical microscale problem, the fine scale stiffness tensor has to be defined. In Voigt notation, i.e., after compressing the symmetric coefficients, the stiffness tensor can be represented as a symmetric matrix, which in turn can be stored as a vector with the entries C ¼ ðC 11 ; C 22 ; C 33 ; C 12 ; C 23 ; C 13 Þ. The six load cases read C 6 ¼ ð1;0; 0;0; 0;0Þ, C 7 ¼ ð0;1; 0;0; 0;0Þ, C 8 ¼ ð0;0; 1;0; 0;0Þ, C 9 ¼ ð0;0; 0;1; 0;0Þ, C 10 ¼ ð0;0; 0;0; 1;0Þ, and C 11 ¼ ð0;0; 0;0; 0;1Þ. For the heat conduction problem, the initial condition is ambient temperature, i.e., T D ≔ T a ¼ 298.15 K. The Neumann boundary is chosen to be isolating on all boundaries except for the bottom one, on which a heat transition condition is set, i.e., q N ≔ 10ðT − T a Þ. Finally, the source is given by Q ≔ AP h ∕ðV el n eph Þ with the absorption coefficient A ¼ 0.4, the hole power P h ¼ 0.0163416 W, the volume of the hole V el n eph , and the number of elements per hole n eph . For the mechanical problem, there is neither an initial condition nor a source. Instead, two boundary conditions have to be provided. The foil is clamped uðxÞ ¼ ð0;0; 0Þ T at x ¼ ð−1000; −1000;0Þ T . Additionally, it is constrained in the x-direction at x ¼ ð−1000;1000; 0Þ T and in the y-direction at x ¼ ð1000; −1000;0Þ T . There are no traction forces applied to the surface; therefore, c t i ¼ 0.
To solve the time dependent system Eq. (82), the time step is chosen to yield a Courant-Friedrichs-Lewy number of C ¼ 1000 and a total of 400 time steps is simulated. Then, the system is solved using the GMRES 17 algorithm and the ILU0 18 preconditioner from the PETSc project. 19 The algorithm stops when either a tolerance of 1 × 10 −6 is reached or the number of iterations surpasses 100. The structural mechanics problem Eq. (92) is solved using the iterative Newton SOR method 20 with a relaxation factor of 1.5.

Results
The solutions to the microscale problem were computed once per load case. This information was then fed into the patch-resolved simulation. Therefore, there are no comparisons to the hole-   Fig. 4). Figure 5 shows the displacement influence function for the first three load cases. Again, it can be seen that the geometry and the stiffness tensor play the dominant role as the displacement influence function in the right most plot is not as affected as in the other two cases.
In Fig. 6, it can be seen that the temperature field is highest in the middle of the structured area. Even though the ablation area is rectangular, there is an almost radial symmetric temperature dissipation. Note, however, that the midpoint of this radial symmetric field is slightly shifted toward the right and the top of the structured area. This is an effect of the hatch strategy, which drives the laser from the lower left to the upper right corner. Because the heat cannot dissipate into the holes, the cooling rate is highest at the corners, high at the border of the pattern, and lowest at the midpoint. Hence, it is a radial symmetric temperature field.  The displacements in Fig. 7 are higher in the positive x or y direction, respectively. As the deformation is solely driven by the temperature differences, this is again an impact of the hatching strategy. The elastic deformations in the left or bottom edge are already reducing to the initial configuration, whereas the deformation at the upper and right edges have just occurred. The somewhat lower displacements at the corner are a direct result of the clamping boundary conditions, which act on the bottom side of the foil.
The von Mises stresses in Fig. 8 are lower at the corners and the edges of the foil and highest in the middle. This actually visualizes the clamping conditions and explains the reduced deformations toward the corners. Because the process is temperature based, it is obvious that the von Mises stress show a similar pattern as the temperature field.
Finally, it is observed that the results of both the hole-resolved and the patch-resolved simulations coincide. The total runtime for the hole-resolved simulation was 36 067 s. The patchresolved simulation took 421 s to compute the RVEs and 23 s to simulate the foil using the RVEs. In total, it is an improvement of over 98%, but because the RVEs only need to be computed once per material and hole shape, the improvement for further runs is even better at over 99%. Fig. 6 The temperature of the hole-resolved simulation is displayed in Fig. 6(a) and the multiscale simulation in Fig. 6(b).

Conclusion and Outlook
To develop the necessary process understanding to allow for a first-time-right production of FMMs, it is essential to have a fast simulation that aids in the search for better processing strategies. Therefore, a multiscale approach was developed within this treatise. This mathematical model consists of a two-scale heat conduction problem and a two-scale thermo-elasticity problem. In total, four tasks needed to be solved, two for the heat conduction and two for the structural mechanics. The discretization in space was performed using a Bubnov-Galerkin method, and the time discretization of the heat conduction problem was realized with the Euler method.
The implementation of the multiscale model is capable of reproducing the results simulated with a microscale model in a fraction of the time. The total achieved improvement in runtime is over 98% with an absolute runtime of a few seconds for the shown experiment.
The multiscale simulation enables scientists and researchers to explore and evaluate different scanning strategies with respect to the temperature-induced distortions in the work piece. Additionally, the processing strategies for filters or FMMs can be found faster compared with experiments, which might yield a quicker production ramp-up.
The current model does not account for plastic deformations. However, the yield stress is a good indicator for the occurrence of plastic deformations and can therefore be used to prevent situations in which the model fails to work. In future research, the authors plan to validate the model against physical experiments. In addition, the automation of the chain of calibrations and simulations is desired. Finally, the simulation might be used as a basis for an optimization algorithm, which selects the best scanning strategy for structuring thin foils.