ASYMPTOTIC EXPANSIONS AND EFFECTIVE BOUNDARY CONDITIONS: A SHORT REVIEW FOR SMOOTH AND NONSMOOTH GEOMETRIES WITH THIN LAYERS

Problems involving materials with thin layers arise in various application fields. We present here the use of asymptotic expansions for linear elliptic problems to derive and justify so-called approximate or effective boundary conditions. We first recall the known results of the literature, and then discuss the optimality of the error estimates in the smooth case. For non-smooth geometries, the results of [18, 57] are commented and adapted to a model problem, and two improvements of the approximate model are proposed to increase its numerical performance. Résumé. Les problèmes faisant intervenir des matériaux avec couches minces sont fréquemment rencontrés dans divers domaines d’application. On présente ici la dérivation et la justification de conditions aux limites approchées, à l’aide de techniques de développements asymptotiques. Après avoir rappelé les résultats connus dans la littérature, on discute l’optimalité des estimations d’erreurs dans le cas d’une géométrie régulière. Dans le cas non régulier, les résultats obtenus dans [18,57] sont commentés et adaptés au cas d’étude. Deux améliorations de la condition approchée sont proposées pour augmenter le taux de convergence.


Introduction
Problems involving materials with thin structures arise in various application fields.This is the case for the analysis of mechanical properties of thin rods, beams, plates, or shells for which reduced models are derived, allowing to deal with lower dimensional geometries without thickness.We will focus in this paper on problems where the thin structure lies around another material, or inside materials, and has significantly different properties.We have in mind a large variety of applications.In mechanical engineering, the study of the properties of composite materials is a critical issue, see [32,33], with the reinforcement by thin fibers or layers.Situations where two materials are glued together enter this scope as well, see e.g.[34].For electromagnetism, thin dielectric layers appear in many situations, see for example [38] for the eddy current problem in the context of copper deposites on tubes, or [61] for the skin effect problem, which has strong connections with thin layers.Biological tissues often involve thin parts, see [48] for a mathematical and numerical study of the electromagnetic field around and inside a biological cell, or [21] for the description of the diffusion Magnetic Resonance Imaging signal in biological tissues.Thin films are good examples as well, and various models have to be considered depending on their nature and size, see [53] and the references therein for falling films, or [56] for an electrochemical situation.We can mention large scale applications in geophysics, where the earth crust may be considered as a thin layer, see [16,51].This list is not exhaustive and many other applications could be cited.
The common issue raised by such problems is the following.If we completely omit the thin layer in the (analytical or numerical) study, the obtained solution may significantly differ from the expected one.On the other hand, incorporating directly the thin layer generally prevents from analytical results, and leads to serious difficulties in the numerical simulations.Indeed, the discretization of the domain needs a local refinement at the scale of the layer, and due to the number of degrees of freedom, the computation can become cumbersome, especially for three-dimensional problems.An alternative solution consists in replacing the initial transmission problem with an approximate model, where the thin layer no longer appears, but is replaced with a suitable approximate boundary condition -or approximate interface condition, depending whether the layer lies inside or around the medium -also referred as to effective boundary condition (or impedance condition in electromagnetics), see Figure 1.
This strategy has generated a large amount of mathematical studies for both the derivation, the justification, and the analysis of such approximate boundary conditions.During the nineties, many geometrical situations have been investigated, from the simplest problem of a layer with uniform thickness around a material [11,31], or non-uniform thickness layers [5], to the more general case of periodically oscillating layers [3,44] (wall laws for flows over rough surfaces) or [6] (scattering by thin periodic coatings).Likewise, different models have been considered: stationary Laplace-Dirichlet and Helmholtz problems in [11,31], harmonic Maxwell equations in [5,31], time dependent Maxwell problem in [39], Stokes system in [3,44].Let us finally mention some works which directly consider the thin layer problem, and develop adapted numerical strategies [20,36], and a nice paper on the problem of optimizing the thickness function [2].
The subject still generates active research activities.Let us mention [1,10,15] for mechanical applications, [12,48,54] for electromagnetic problems, [59] for the Stokes system, [46] for the heat equation, [17,29] for general purposes.On the other hand, the case of a random thickness has been investigated, see [9] in the context of rough surfaces, and [23] for a practical application of approximate boundary conditions to compute moments of solutions of boundary value problems inside random domains.Let us mention the works [4,19,49] on polarization tensor for thin inclusions of rough layers.
ε We restrict here ourselves to layers of uniform thickness.For a review on problems with rapidly oscillating layers, see [50].See also [28] for an example mixing homogenization and matched asymptotic expansions.

Approximate condition
For the derivation of approximate boundary conditions, several techniques have been used, most of them based on asymptotic expansions of the solution in the domain with thin layer.Obtaining error estimates is a more or less difficult task, depending on the chosen tools.In section 1, we detail the approach of [31] and [11], and discuss the obtention of estimates.The detailed analysis of a complete asymptotic expansion is then given, allowing the investigation of the optimality of the error estimates.This description is done for a very simple problem in order to make the presentation the clearest as possible.
The second section §2 concerns the investigation of approximate boundary conditions in the case of nonsmoooth geometries.Indeed, impedance boundary conditions are very common in electrical engineering, and are often used for geometries with corners or edges.Users are aware of their limitations, see [7,42,52,55,58].In [57], the question of the performance of the impedance condition in the presence of corners has been rigorously investigated.A positive answer is given for the convergence of the approximate problem to the transmission one, as the thickness goes to 0. But it is shown that the performance of the impedance condition is weakened by the presence of corners (the analysis is done for the two-dimensional Laplace equation).To our knowledge, no mathematical work has been done since, which describes a method to overcome this difficulty and recover the convergence rate of the smooth case.

Approximate boundary conditions for smooth geometries
In order to present and compare the different approaches, we concentrate on the following toy problem where the domain Ω ε is the unit square with an ε-layer on its top side, and f is some L 2 (Ω − ) given datum.The border domains are defined in Figure 2. Problem (1) is associated with the following variational formulation: find with the variational space V = {ϕ ∈ H 1 (Ω ε ) ; ϕ = 0 on Γ − ∪ Γ ε + }.This variational problem obviously enters the framework of the Lax-Milgram lemma, and thanks to the Poincaré inequality, we can easily obtain the following a priori estimate with a constant C > 0 independent of the parameter ε, due to the fact that the domain Ω ε is uniformly bounded with respect to ε.

A variational approach
In [31], Engquist and Nédélec propose an elementary construction of approximate boundary conditions for Problem (1) (and generalizations for Helmholtz and harmonic Maxwell equations).The basic idea consists in assuming that the solution u ε has a polynomial behavior in the transverse variable y (see Figure 2) -this fact will be discussed later on.More precisely, if we assume that u ε + (s, y) ∈ P 1 (y), then it reads Using the transmission conditions Exploiting now the external Dirichlet condition on Γ ε + , we obtain This suggests to introduce v ε solution to the problem which is a well-posed boundary value problem, where the effect of the thin layer of Problem (1) has been replaced with a Robin boundary condition on the interface Γ.
A similar derivation can be performed, assuming that u ε + is a polynomial of order 2 in the y-variable.Actually, the same approximate condition is recovered.In the more general case, where the geometry of the layer is curved, the condition reads This approach is very simple, but does not give any way of estimating the difference between the solution of the original problem u ε − , and the solution of the approximate problem v ε .To overcome this drawback, the authors propose a derivation based on the variational formulation (2).Considering test fonctions v which have a polynomial dependence on the transverse variable y leads to the same approximate boundary conditions on the interface Γ as above, and is suitable for obtaining error estimates.Precisely, the following is proved in [31] (for a flat boundary, i.e. c(s) = 0) In [11], Bendali and Lemrabet propose the derivation of an asymptotic expansion of u ε with respect to the thickness ε of the layer.This leads to the derivation of approximate models, together with error estimates.The key idea is to use a dilation in the layer: introducing the scaled variable Y = y/ε, and this variational formulation allows to identify the terms u n and U n .The derivation of approximate boundary conditions follows from the particular form of the first terms in the expansion (see §1.2 for more details on this point), leading to the same boundary conditions as in [31].
Again, error estimates are given, coinciding with (5).The same technique has been used in [31] and [11], which can be briefly summarized as follows.The variational formulation (2) of the original problem reads a(u ε , ϕ) = (ϕ) for ϕ ∈ V , whereas the approximate problem (4) can be written as and the approximation of V by V 1 , leads to (5).
However, numerical tests suggest that the error estimate ( 5) is suboptimal.Indeed, a convergence rate in ε 3 is observed, rather than the expected ε 5/2 , see Fig. 3.In the next section, we show how to obtain this optimal error estimate.

Exploiting an asymptotic expansion: optimal estimates
We present here the construction of a complete asymptotic expansion for Problem (1), and the derivation of optimal error estimates for the approximate model (4).Similarly to [11], we postulate the Ansätze 1 We do not seek any convergence for the series with respect to n.The point of view will consist in truncating the series at a fixed order, and let ε tend to 0. This corresponds to the so-called convergence in the sense of asymptotic expansions.
Inserting this power series into Problem (1), we get a cascade of uncoupled problems: for Y ∈ (0, 1), with the convention u −1 = 0 and U −1 = U −2 = 0, and where δ n,0 denotes the Kronecker symbol.Note that the problem solved by U n is one-dimensional, and the tangential variable s is only a parameter.Besides, the Neumann condition on vertical sides ∂ ν U n = 0 on Γ L ∪ Γ R will be fulfilled, even if this is not clear at this stage.
The question of the solvability for Problems ( 6) is rather simple: since f ∈ L2 (Ω − ), the condition s → U n (s, 0) ∈ H 1/2 (Γ) ensures existence and uniqueness of u n in the H 1 -variational sense.This condition is fulfilled as soon as u n−1 belongs to H 2 (Ω − ) (note that no smoothing effect in the tangential direction can be expected from the resolution of the one-dimensional boundary problem defining U n ).This regularity "consumption" requires that u n ∈ H ∞ (Ω − ) for any n, if we wish to build the asymptotic expansion at any order.A sufficient condition is f ∈ C ∞ (Ω − ) with compact support 2 .Under this assumption, a formal derivation of the terms u n and U n is possible for any n ≥ 0.
Let us now show how to obtain error estimates for the asymptotic expansion.For any fixed N ≥ 0, we set By construction, u ε, [N ] solves the following problem Comparing with Problem (1), and using an a priori estimate similar to (3), we immediately get (for any integer N ) This estimate can be improved with the mere observation that leading to the optimal estimates Once this study has been performed for the transmission problem (1), it is possible to do the same for the approximate problem (4).Starting from the ansatz v ε = n≥0 ε n v n , we get the iterative construction of the terms thanks to the following problem.
Again, under the same smoothness assumption on the right-hand side f , this is possible to build the terms v n , and the following optimal estimate is obtained: Now, it becomes possible to compare the solution of the transmission problem u ε − in the limit domain Ω − with the solution of the approximate problem v ε .Indeed, an explicit computation of the first termes shows that Note that u 2 = v 2 in the more general case where the curvature on Γ does not vanish.Then, using estimates (7) and ( 9), we immediately get which corresponds to the numerical evidence.Let us point out that uncoupling the estimates like in (7) between the different part of the domain is the key ingredient to obtain optimal estimates.
Remark 1.1.In the whole paper, we consider as our goal to build approximate boundary conditions which replace the effect of the thin layer.This point can be discussed.Indeed, the construction of the first terms in the asymptotic expansion (see for N = 1 just above) can lead to a numerical method with the same convergence rate.Of course, this requires the resolution of two Dirichlet problems instead of one Robin problem.The difference is not huge, and the advantage would be that, once u 0 and u 1 are computed, the approximation u ε,[1] − = u 0 + εu 1 can be directly computed for any value of ε, without any extra computational effort.

The case of nonsmooth geometries
We investigate now the case where the limit problem is non-smooth.Precisely, we consider the geometry given by Figure 4, where the thin layer lies only on a part of the upper side of the square (Neumann boundary conditions are imposed on Γ N ).At a first glance, it may not be obvious that the situation of Figure 4 is less regular than this of Figure 2. Actually, if we consider the limit case (where ε = 0), then the solution u 0 solves the following problem in Problem ( 11) is a mixed Dirichlet-Neumann problem where singularities at the point 0 (where boundary conditions change) are limiting the regularity of the solution.It is well known -see [25,37] -that u 0 admits the following splitting into singular and regular parts where the regular part u 0 reg belongs to H 2 (Ω − ), the singular coefficient γ is a scalar, and the singular function s 1 2 is given by s(x) = √ r sin θ 2 , where (r, θ) denote the polar coordinates centered in 0 such that Γ corresponds to θ = 0.If c = 0, the function u 0 does only belong to H 3 2 −δ (Ω − ) for any δ > 0. The condition f ∈ C ∞ (Ω − ) with compact support is no more sufficient to ensure u 0 ∈ H ∞ (Ω − ), preventing the construction of the asymptotic expansions (for both the transmission and the approximate problems) as exposed in Section 1.2.

The first terms in the asymptotic expansions
In [22], the asymptotic expansion for the approximate problem (where a Robin condition is imposed on Γ) has been built.At first order, it takes the form where the profile Z solves the following problem in the half-plane Π − (see Figure 5) and v 1 can be computed from u 0 reg and Z, via a problem similar to equations (8).It is enough for the following to point out that • The profile Z has the following behavior as |X| → ∞ (θ denotes the polar angle centered at 0) where with some terms v 1,i which do not depend on any coefficient.
The infinite domains for the definitions.
Meanwhile, it is possible to extend results from [18] to the transmission problem in the domain of Figure 4 (see also [43] for an approach with matched asymptotic expansions of similar problems).The solution u ε admits the following expansion in Ω − where K solves the following transmission problem in the infinite half space with semi-strip of thickness 1, see Fig 5 Again, the following description will be useful • The profile K has the following behavior as |X| → ∞ • The term u 1 admits the splitting

Lack of performance of classical approximate boundary conditions
In the case of a corner domain with a thin layer, is has been proven in [57] that the performance of the approximate boundary conditions is much less than in smooth domains.The situation is the same here: following the same technique, we can show that expansions ( 17) and ( 13) compare as follows Moreover, using ( 15) and ( 19), we immediately get As a direct consequence of ( 21) and ( 22), we obtain the following estimates These estimates are generically optimal, and illustrated by numerical simulations, see Figure 6.

An attempt towards improved approximate boundary conditions
The aim of this section is to present some ideas to improve the approximation of the thin layer problem.The starting point is relation (21), where it is clear that the main terms are as follows The correction using the profile K − Z is the concern of §2.3.1.
• In the L 2 -norm, the terms γε(µ − λ)u 1,2 (x) and γ √ ε(K − Z) x ε are both limiting.But an explicit estimation shows that for some constant c.The (extraordinary) situation where µ = λ would be particularly favorable.This remark will lead to the strategy developed in §2.3.2.

Correction with transmission profiles
A first idea consists in correcting the approximate solution v ε with the difference γ √ ε(K − Z) x ε , the profiles Z and K being computed as a preliminary step.This kind of method has been introduced in [14,24,26] for similar problems, see also a remark in [27,Appendix B.5.2].Here are the steps of the algorithm Compute Z and K by solving problems ( 14) and (18).Since these problems are posed in an infinite domain, a specific numerical strategy is needed.Many possibilities been investigated in the literature, such as the infinite elements method introduced in [13,60], which directly incorporate the infinite part of the domain inside the Galerkin method.The other methods start to bound the domain with an artificial boundary on which some appropriate boundary condition is imposed, see [30,35,40] in the context of absorbing conditions, or [45] for integral representations.We have here chosen to bound the domain with a disk of radius R and impose on the circle the following condition, which is suggested by relation ( 15) D + 2R∂ ν D = 0, with D = Z − s − , This computation is made on a fixed grid thanks to a nodal finite element method in the presented simulations.
• Step 2. For fixed ε, compute the solution v ε to the approximate problem on a mesh of Ω − .
on the mesh where v ε has been computed (this needs a mesh transfer).
which gives a better approximation in Ω − of the solution u ε of the transmission problem.
Of course, the computational cost of Step 1. needs to be relatively cheap!But remember that the computation of the profiles Z and K can be done once and used to approximate different transmission problems with various right-hand sides f , and values of ε.
The above discussion about the asymptotic expansions -see ( 21) -leads to the following error estimates which show that only an improved performance can be expected for the H 1 -norm.This is in accordance with the numerical tests.In Figures 8 and 9 the isovalues of the differences u ε − v ε (without correction) and u ε − ṽε (with profile correction) are plotted.This is clear that the behavior near the point 0 has been correctly taken into account.
Remark 2.1.With the same idea of correction by a profile, we can consider the approximation vε The advantage is to avoid the resolution of the Robin problem, replaced by the limit (Dirichlet) problem.From the computational cost point of view, the difference is not very relevant, but if we wish to perform several simulations for different values of ε, the computation of u 0 is done once and for all.

A multiscale approximate boundary condition
In general, the coefficients µ and λ arising in the behavior of the profiles Z and K at infinity, do not have the same value.Our idea consists in replacing α in the condition Z + α∂ ν Z = 0 on G with a function α • : R + → R + , given by (see also Fig. 10) where ρ 0 and C α are non-negative constants.The profile problem ( 14) is then replaced with and, similarly to (15), we have Remark 2.2.The coefficient α • is not a smooth function.This, however, does not lead to the apparition of strong singularities (only the derivatives of order 3 are really affected).A smoothing would be possible, but is not really necessary.
Our aim is to see if we can adjust ρ 0 in order to have λ • = µ in (22).A direct consequence will be an improvement of the L 2 -norm estimate as follows 2 ), where v ε • solves the following problem As above, we compute the profile Z • on a domain with a disk of radius R and impose, from (25), on the circle the following condition: which is easy to compute once we know Z • and s and .The same expression is available with µ, replacing Z • with K.
A comparison principle on the bounded domain allows us to prove that the mapping α • → Z • is increasing and continuous.For the parameters α = 2 and C α = α 2 , we observe that λ • = µ for the value ρ 0 0.02, see Fig. 11.In Fig. 7, the convergence rate of the quantity u ε − v ε • L 2 (Ω−) is plotted with respect to ρ 0 .We observe that a better convergence rate is obtained for ρ 0 0.02, which is consistent.Let us concede that the rate H 1 -norm has not changed.Nevertheless a closer look to the absolute errors show that for the value ρ 0 0.02, even the H 1 -norms are smaller than for ρ 0 = 0 (standard approximate boundary condition).More details on the presented techniques will be found in the forthcoming article [8].
The presented numerical simulations have been performed thanks to the finite element codes Mélina [47] and FreeFem++ [41].• with respect to ε for the nonsmooth case.

Figure 1 .
Figure 1.Replacing the effect of the thin layer by an approximate/effective boundary condition.

Figure 2 .
Figure 2. The model domain for Problem (1) in the smooth case.

Figure 3 .
Figure 3. Error u ε − v ε inside Ω − with respect to ε for the smooth case, see Figure 2 (f = 1 for the simulation).

Figure 4 .
Figure 4.The model domain for Problem (1) in the nonsmooth case.

Figure 6 .
Figure 6.Error u ε − v inside Ω − with respect to ε for the nonsmooth case, see Figure 4.

Figure 12 .
Figure 12.Convergence rates of the L 2 -norm of u ε − v ε• with respect to ε for the nonsmooth case.