Nonlinear transmission conditions for thin curvilinear low-conductive interphases

In this paper, we consider the heat transfer problem in a cylindrical composite material with an adhesive interphase of closed curvilinear shape. The interphase exhibits nonlinear behaviour and at the same time has physical characteristics (size, thermal conductivity) that significantly vary from the properties of the surrounding components. As the latter complicates the direct use of FEM, we use asymptotic method to replace the interphase with an imperfect interface of zero thickness, preserving the essential features of the thermal behaviour of the interphase through the evaluated transmission conditions. Carefully designed numerical simulations verify their validity. We place special emphasis on the impact of geometric properties of the interphase, in particular, the curvature of its boundaries, on the accuracy of the conditions.

Thin arbitrary curvilinear interphase idea, the thin soft elastic layer between two solids was represented as a zero thickness object to study into surface/interfacial waves in [8]. Similar methodology was also developed for thin coatings, for example, see [4].
Interestingly, such type of problem in application to various physical phenomena has also inspired the development of mathematical theory [5].
The transmission conditions derived in this paper can be used for modelling heat transfer or mass transfer in a thin layer [15]. In the setting of mass transfer, transmission conditions can be used by analogy once one replaces Fourier's law with first Fick's law and heat equation with second Fick's law. However, we shall from now on consider the problem in terms of heat transfer only. The generation or absorption of heat in the thin layer can happen due to internal chemical or electrochemical reactions, for example, in SOFC cells [10,24] or as a result of mechanical work [1]. In the latter case, the problem of heat transfer is coupled with that of plastic deformation.
Physical problem similar to the one described in this paper has been considered in [20] and verified in [21] for materials with thin layers with flat, rectangular surface for the case of monotonic temperature distribution in the interphase. Eliminating this assumption of temperature monotonicity, more general transmission conditions were obtained in [18,22]. Real-life applications and use of composite materials, however, make it necessary to consider media with interphases of closed curvilinear rather than flat shape. For the first time, the transmission conditions for such a geometry were obtained and verified in [2], where the interphase was assumed to be a thin ring, or at least have a form close to circular, which significantly simplified the formulae. At the same time, the increase in curvature of the interphase led to an observable decrease in the reliability of the conditions. The present paper achieves the necessary extension of [2] to a more general geometry concentrating on the case of thin interphases being represented as domains bounded by two smooth closed curves of small curvature. To be precise, we consider an interphase of curvilinear shape close to that of a four-ray star, as its geometry allows us to investigate an effect similar to the edge effects for planar domains, where singularities in the solution may occur due to the failure of initial assumptions or mismatch of the physical properties of the materials at the free edge [16,19,26]. The inclusion here consists of a material of low conductivity as compared to the two surrounding (not necessarily the same) materials. The following section contains the problem formulation and gives some preliminary steps to the process of evaluating transmission conditions, which is presented in Sect. 3. The results of the numerical computations that support our conclusions are presented in Sect. 4. Finally, Sect. 5 provides the summary of the results.

Problem formulation
Consider a composite consisting of two materials joined by a thin curvilinear interphase (see Fig. 1). Within the latter, the heat transfer equation is satisfied Here, T (x, y) is the unknown temperature, k(T, x, y), c, ρ are, respectively, the thermal conductivity, the heat capacity and the density of the material, and Q(T, x, y) the thermal source or sink spread within the interphase. Note that it is assumed that the sign of Q does not change within the interphase, as there is no point in placing multiple heat sources and sinks within such a thin region.
Along the boundaries of the interphase (smooth closed curves denoted as Γ ± ), perfect transmission conditions have the form: Here, n is the normal vector to the surface and q the heat flux, defined according to Fourier's law: q = −k(T )∇T.
Within each of the three parts of the domain, the temperature distribution is different and is described with the functions T (x, y, t) within the interphase and T ± (x, y, t) within the surrounding media adjacent to the boundaries Γ ± of the interphase, respectively.
The thermal conductivity of each material is denoted as k + , k, k − , where k k ± and, in general, k + = k − .

Preliminary steps
Transition to polar coordinates (r, φ) is a natural and convenient step considering the geometry of the domain [28]. We also introduce r = r 0 (φ)-the centre line Γ between the boundaries of the interphase and h(φ) = r + (φ) − r − (φ)-the width of the interphase.
To obtain the transmission conditions for the problem considered, the convenient method [2,18,20] is to rescale the thin interphase, as it is much thinner than the adherent media. This is done by rescaling the width of the interphase and simultaneously introducing the new coordinate ξ : while the second coordinate φ is unaltered (as in [2]). Taking into account the proportions between the physical parameters, the thermal conductivity of the interphase and the heat source are also rescaled where ε 1. After the rescaling, the temperature distribution within the interphase is also defined in terms of the new coordinates: while the normal vectors to the boundaries Γ ± are described by the formula Following all these transformations (for details, see [2]), the heat equation and the transmission conditions take the form Here and further, the notation h = h(φ), r 0 = r 0 (φ) is adopted for brevity.

Evaluation of the transmission conditions
Using the asymptotic method, we are looking for the solution within the interphase in the form At the same time, the normal vectors to the interphase boundaries are also expanded asymptotically as where Substituting the expansion (11) into (8) gives a chain of boundary value problems for the corresponding terms T j (ξ, φ, t): Here, we provide just the first two of the problems, demonstrating that only the first of them is nonlinear while other ones are linear. As we shall from this point limit ourselves to the leading terms of the temperature and the heat flux, only the equation (14) needs to be solved. Substituting the expansions (11), (12) into (9-10), and, again, restricting to the leading terms only, results in the boundary value problem Observe that this formulation of the problem takes into account the geometry of the interphase, which may be very different from a circle, unlike the analogous system of equations in [2].
We have obtained the transmission conditions in an explicit form for several special cases presented in the following subsections.

Transmission conditions for coordinate-dependent thermal conductivity and heat source
In the first instance, neither the thermal conductivity nor the heat source depends on the temperature, i.e.
Integrating (16) with respect to ξ gives Let us introduce auxiliary functions Then, substituting ξ = 1 2 into (19) and making use of (18) gives or which is the first transmission condition.
For the other condition, we rewrite (19) as Then, we introduce two more auxiliary functions and, by integrating (23), obtain the second transmission condition Supposing that both the thermal conductivity and the heat source are constant, the transmission conditions (22) and (25) are simplified to the form Remark In the forthcoming cases, assuming that Q, k are constant leads, after some transformations, to the same set of conditions (26).

Transmission conditions for temperature-dependent thermal conductivity and coordinate-dependent heat source
In this subcase, the assumption is that This case also includes the instance of k = k(T, ξ, φ, t) as a particular case.
Here, we shall evaluate two different sets of transmission conditions, depending on the assumption whether the temperature distribution within the interphase is monotonic or not. Thus, this instance is split into two subcases.

Transmission conditions with the additional assumption of monotonic temperature in the interphase
In this case, the first transmission condition, as well as its evaluation process, is exactly the same as in the previous one. To obtain the second condition, we write by analogy with (23) Then, we introduce auxiliary functions . (29) Integrating (28), with the change in integration variable in the left side of the equation, we get the second transmission condition

Transmission conditions with no additional assumptions on the temperature in the interphase
We suppose in this subcase that there exists a unique point of extremum ξ * for T (ξ ) in the interval − 1 2 , 1 2 (the uniqueness of the extremum has been demonstrated in [2]). Then, using the introduced auxiliary functions in the appropriate parts of the domain and integrating the initial differential equation (16) in the respective subdomains will lead to r 2 Here, we have taken into account that ∂ T ∂ξ ξ * = 0.
Simplifying (31) leads to Each of these formulae can also be used for evaluating the extremum point ξ * . Then, equating the two expressions for ξ * , the first transmission condition is obtained in the form Introducing the notation T * = T (ξ * ), we get by means of integration Then, the second transmission condition is where ξ * is substituted from (33). Meanwhile, both sides of (34) provide a formula for the maximum/minimum temperature value T * .

Special case (2): Q
The substitution of such expressions for k, Q into the transmission conditions provides the same set of equations regardless of whether we make the additional assumption on monotonicity or not. While the first transmission condition remains the same as in 26, the second one is given by the same equation as was obtained in the analogous case in [2]:

Transmission conditions for temperature-dependent thermal conductivity and heat source
Here, both heat conductivity and heat source are nonlinear, i.e.
This case, like the previous one, is split into two subcases, depending on whether we assume the monotonicity of the temperature or not.

Transmission conditions with the additional assumption of monotonic temperature in the interphase
In the first place, we transform the equation to solve into where Lowering the order of the equation, we obtain where we've introduced the auxiliary functions as Taking the value at the boundary, we get or Then, the first transmission condition takes the form To get the second transmission condition, (38) is transformed to the view − r 0 h Introducing the auxiliary functions and integrating (41), we have which is, in fact, the second transmission condition. By analogy, an equivalent condition can be obtained in the form It should be noted that the definition of the functions Υ ± (T ) implies that The latter condition is equivalent to having within the entire domain, i.e. at every point with every value of T.

Transmission conditions with no additional assumptions on the temperature in the interphase
As we don't assume anymore that the temperature within the interphase is monotonic, we reintroduce the notations T * , ξ * that correspond to the maximum/minimum value of the temperature and the point at which it is achieved. Following a similar process to the one described in the previous subcase with assumed monotonicity (37-38), we get an expression analogous to (39) where we've taken into account that k 2 (T * )(T * ) 2 = 0. Equation (46) is equivalent to q 2 ± = 2Φ ± (T * ). Expressing T * from these formulae, we obtain the first transmission condition To derive the second transmission condition, we again use the previously introduced functions Υ ± (T ) and the same logical approach. Then, in the corresponding parts of the domain, we have From this, we can get the formula to evaluate the extremum point ξ * as and the second transmission condition

Special case (3): Q(T
In this instance, we are omitting cases of negative β to avoid the situation when Q 0 = 0. Besides, the second transmission condition is obtained in different forms for the cases of heat source ( Q 0 > 0) and heat sink ( Q 0 < 0). In the former case, an appropriate substitution leads to the following equations under the assumptions of monotonicity. Without this additional assumption, the first one transforms into the same equation as (51 1 ), while the second condition has the form If we consider the case of heat sink, however, the second condition is with or without assuming temperature monotonicity, respectively. Appendix provides the details on how these conditions were evaluated.

Numerical examples
To validate the accuracy of the transmission conditions we consider a numerical example in which the domain is a cylinder with a curvilinear interphase of shape resembling a four-vertices star (see Fig. 2). The boundaries of the interphase are formed by segments of circles of radius r − = 1.01 (for the inner boundary) and r + = 1 (for the outer boundary) with centres at (1; 1), (−1; 1), (−1; −1), (1; −1), where the outset of the coordinates is the centre of the cylinder. The sharp ends of the interphase were approximated with smooth curves in order to avoid singularities. The boundaries of the cylinder are circles of radius R = 1.1 for the outer boundary, r = 0.1 for the inner one. The temperatures set at these boundaries are, respectively, 300 and 320. Like in [2], the thermal conductivities are taken as 237 for the outer media and 0.2 for the material of the interphase, while the heat source is Q(T ) = 180(50 + T ). We take ε = r − − r + = 0.01 for all further computations.
It should be mentioned that the choice of the size of the interphase enabled us to have accurate COMSOL simulations as well as approximations through transmission conditions, while for narrower region FEM failed. At the same time, as has been observed for a circular interphase in [2], the model with transmission conditions has even better accuracy for smaller ε.

Verification of the transmission conditions
For the described case, we cannot carry out the conditions verification explicitly, as we cannot obtain the exact analytical solution to the stated problem. Therefore, we check the accuracy of the transmission conditions indirectly, using the numerical solution from the FEM software COMSOL Multiphysics 4.3b.
One of the ways of such assessment, similar to the procedure used in [17], is to use the transmission conditions to evaluate the physical parameters along the interphase boundaries and compare them to the numerical results from COMSOL. Knowing that the latter are more accurate and reliable in terms of the temperature rather than the heat flux, we take the COMSOL-evaluated values of T ± and substitute them into the conditions. Solving thus obtained systems of equations, we get the values q (m) ± , q (nm) ± (corresponding to the conditions derived with and without assumptions of monotonicity), and compare them to q ± from COMSOL.
In Fig. 3a,b, where the heat fluxes evaluated from the two sets of transmission conditions are plotted alongside the numerical solutions for heat flux, we can observe that the difference between such solutions is barely distinguishable by sight, apart from at the points approaching the ends of the interphase. The reason for these seemingly unreasonable values is the geometry of the domain in COMSOL. As the sharp notches of the interphase were "cut off" by the approximating curve, there is a mismatch between the points of the outer and inner boundaries, and therefore, the transmission conditions cannot be used in those regions. Eliminating such probe points, we see that they are followed by a small part of the domain where the error is different from the error in the more central part of the interphase. To demonstrate this edge effect, the absolute and relative errors for heat flux are shown in Fig. 3c,d. Both the vertex region with huge error and the zone with edge effect can be clearly seen.
The type of the mesh, as turned out while trying several options, from extremely fine to extremely coarse, doesn't have much impact on the eventual result. The relative error always remains at roughly 10 −2 −10 −3 (i.e. the order of ε) in the central parts of the boundary yet visibly changing in the edge zones. This proves that our initial assumption on the curvature of the boundaries was indeed a crucial one.
This phenomenon can also be seen if we assess the transmission conditions by taking the values of both T ± and q ± as evaluated by COMSOL and substituting them into the conditions. Then, we calculate the absolute and relative errors between the left and right sides of the expressions, again observing the significant growth of errors closer to the vertices of the interphase. Table 1 gives the maximum absolute ( ) and relative (δ) errors evaluated for the entire set of probe points (with the edge effect zones) and for the cropped set (without the edge zones). Here, as before, the "+" and "−" signs are used for the outer and inner boundaries, respectively, while "m" and "nm" stand for the set of conditions used to evaluate the heat flux (the ones derived with or without assumptions of monotonicity). The problematic points at the very ends are not taken into consideration at this stage. To determine the edge effect zone, we applied the 1 % criterion like in [16] and removed the points until the mean relative error δ (m) 1 became  smaller than 0.01. In our case, the size of this zone is roughly 5% of the length of the interphase boundaries.
In the central part of the interphase, outside of the edge effect zones, the relative error is no greater than of the expected order, which is equivalent to the small asymptotic parameter. In the edge effect zones, the relative error is different from the error in the central part. Curiously, it does not necessarily increase, like it does at the inner boundary, but may also go down, as it happens at the outer boundary. This can be explained by the fact that in the regions closer to the vertices of the interphase our assumption on the curvature fails.

Conclusions and discussion
We have obtained the transmission conditions for the general case of a domain with a thin curvilinear interphase under minimal assumptions. The transmission conditions were assessed indirectly, showing good precision (at least the order of the small parameter ε as predicted by the evaluation) for the physical parameters along the interphase. In addition, we have demonstrated that our initial assumption on the curvature of the interphase, despite being fairly general, is both sufficient and necessary, and that the accuracy of the conditions is impaired in the regions where this assumption fails. As a result, special transmission elements can be implemented into FE code to deal with these rather complex interphases. In the edge regions, special FEM elements should be however introduced taking into account the singular edge shapes.
so they eliminate each other in (59), and (50) reduces to the form that is easily turned into the view (52) by a trivial transformation and applying function sin to both sides of the equation. For the case of Q 0 < 0, we follow a similar line of reasoning. However, evaluating the integral like in (57), we get Substituting this into (47) and taking the exponent of both sides of the resulting equation leads to a condition which is equivalent to (53).