Numerical method for solving joint thermo-diffusive problems in an infinite combined domain with thin resistant interphase

This work deals with a class of Boundary Value Problems describing joint thermo-diffussive fields in an infinite combined domain, which consists of two subdomains, matched by a thin intermediate layer. The main problem is reduced to an equivalent one given in the bounded subdomain, with non-local boundary condition on the transmission surface. Such a condition incorporates all the information about the infinite subdomain and the intermediate layer. The equivalent problem is solved by means of Finite Element Method in frames of Matlab package. As it is not possible to introduce the non-local boundary conditions along a part of the boundary directly into FEM code, a dedicated iterative subroutine is constructed. Effectiveness of the method has been checked on selected benchmarks. Accuracy and convergence of the procedure have been addressed in the analysis.


INTRODUCTION
Nowadays, in many problems of modern engineering and technologies a multitude of advanced, composite materials are dealt with.Increasingly more sophisticated devices and constructions are designed.All these circumstances give rise to a need of effective numerical analysis of different multistructures.Such problems constitute a serious challenge, as they usually involve cases where different parts of computational domains exhibit various geometrical and material properties.Obviously, the shapes of respective subdomains themselves may cause solution singularities in the event of edges, wedges or conical points of the domain boundary (extensive literature can be found for example in [1]- [2]).This creates necessity for effective numerical evaluation of the singularity and related coefficients (stress intensity factors in mechanics e.g.[3], [4]).In case of inhomogeneous structure, the Domain Decomposition Methods proved their effectiveness [5], [6].In this approach the main problem is split into a few component problems in respective subdomains, solved simultaneously with appropriate compatibility conditions at the division surfaces.Such conditions are considered in the subdomains as the boundary conditions [7].Subdomain partition is carried out depending on the types of governing equations and/or the material properties of different parts of the domain.Additionally, when the domain boundary shape is very complicated, the geometrical decomposition can improve the accuracy and efficiency of computations.Thanks to the aforementioned domain division, for each component problem, the most suitable numerical method may be employed.
Unfortunately, in case of infinite domain DDM approach with standard FEM elements fails.To tackle such problems a specific group of DDM methods has been developed [8], [9].However, the most effective way here is coupling FEM-BEM [10], [12].Each of these component methods has its own range of application, where the maximal effectiveness is obtained.FEM is especially useful for strongly non-linear problems in the domains of very complicated shapes, but finite dimensions.On the other hand, BEM proves its effectives for the linear problems in infinite domains.
Another problem which appears while considering bonded structures is related to thin adhesive layers situated between the components.As a result, different size elements have to be employed, introducing an additional computational challenge.This becomes even more difficult if one deals with an unbounded combined domain, where the interphase thickness is appreciably smaller than the dimensions of the bounded part of the domain.In such cases usually a simplified numerical model is introduced, in which the thin layer is replaced by some transmission conditions suitable for DDM.These transmission conditions usually result from appropriate phenomenological assumptions related to the physical peculiarities of the problem under consideration.
In case of joint fields all the aforementioned features of the problem become even more pronounced.Assumptions valid for one component filed may lead to incorrect result with respect to another component, especially if different regimes are realized for each of them.
In order to avoid such a drawback, the problem with thin interphase between bounded and unbounded subdomains considered in this paper is reduced by means of asymptotic analysis to two limit problems matched with each other [7], [13], [14].One of them (external problem) is defined in the domain composed of the bounded and unbounded part only, where the interphase is replaced by specific imperfect transmission conditions.Solution to the second one is determined within the rescaled interphase.
The structure of the paper is as follows.In the section 2 a stationary problem of joint thermo-diffusive fields in combined unbounded domain is formulated.The domain is composed of three subdomains.One of them is a half-space connected to the bounded region through an interphase layer, which is very thin in comparison with the latter (Fig. 1).Different materials are attributed to respective subdomains.Material properties are highly nonlinear and the fields are joint through the diffusion and heat transfer parameters as well as by possible heat/mass sources.Moreover, material properties of the intermediate layer are essentially different from those of both adherents, namely it is heat and mass resistant.By means of asymptotic analysis, in chapter 2.2 we evaluate imperfect transmission conditions, which incorporate the information about the intermediate layer (by virtue of limited space this description is rather cursory, and obviously it does not contain extensive explanation of all its aspects).These conditions are used in the numerical model to replace the original thin interphase.Naturally, such transmission conditions are also nonlinear and exhibit joint field feature.On the contrary to the linear case, the operators involved in the conditions essentially depend on unknown distribution of the component fields within the interphase.The later in turn, are defined from the asymptotic procedure based on the values of the thermodiffusive field at the interphase boundaries.
For the semi-infinite subdomain analytical methods are applied to obtain a relationship between the solution and its normal derivative at the external boundary.Next, in section 2.3 this equation is implemented in the imperfect transmission conditions along the intephase boundaries to finally obtain a non-classical, non-local boundary condition for the problem in the bounded domain.Such a BVP is solved by means of Matlab FEM block, with use of a dedicated iterative procedure.The feedback between respective fields is provided by the material properties dependence on the temperature and mass concentration.This feature also makes the problem essentially non-linear.The algorithm of solution is described in chapter 3. Proposed method is, in fact, further development of the methods described in [12], [15], [16] where only single physical fields for analogous domain geometry without intermediate layer were considered.
The method used in this paper refers to DDM techniques as the problem has been divided into three component problems in respective subdomains and for each of them different analytical and numerical tools are used.In chapter 4 the numerical examples are presented.All the computations are carried out on the assumption of axial symmetry of the phenomena, what, in fact, constitutes a 2D problem.
The accuracy and convergence rate of constructed procedure are examined on the basis of an analytical benchmark in section 4.1.Next (4.2) a numerical example is presented, and compared with a similar problem, where the feedback between respective fields is neglected.Thus, one may asses the influence of thermo-diffusive coupling on the final result.The final conclusions and remarks are presented in section 5.

GENERAL DESCRIPTION
We consider a stationary multiphysics problem with joint thermo-diffusive fields described by the system of equations given in the semi-infinite domain Ω∪Ω h ∪Π: where the components of the vector-function u(x) = [u 1 (x), u 2 (x)], u 1 and u 2 , denote temperature and mass concentration, respectively, while f is a known bounded vectorfunction f = [f 1 (x, u), f 2 (x, u)], which defines the heat, f 1 , and mass, f 2 , sources.Matrixfunction µ = µ(x, u) is a diagonal one: with the components µ j (x, u) assumed to be positive, bounded and separated from zero for arbitrary x and any vector u.They describe the heat conductivity, µ 1 (x, u), and the mass transfer coefficient, µ 2 (u,x).Heat and mass fluxes through any surface S are defined in the following manner: where n is a unit outward normal vector.The domain geometry considered in this work is shown in Fig. 1.We denote by Ω the bounded part of the domain, whereas the half-space Π is connected to the domain Ω through a thin intermediate layer In what follows we assume that the thickness of the thin layer Ω h is negligible in comparison with the dimensions of the bounded subdomain Ω.The upper and lower surfaces of the interphase are defined in the following manner: z = h 1 and z = −h 2 .Everywhere below we assume that: where 0< ε<<1 is a small parameter, and D Ω is a linear characteristic dimension of the subdomain Ω.Although the asymptotic method allows to consider a general case (z = h 1 (x, y), z = −h 2 (x, y)) under additional constraint of small curvature, we assume for simplicity that h j are constants.Inside the interphase the same equation ( 1) is satisfied.In order to make a distinction between main subdomains and the thin interphase situated between them, we denote respective heat conductivities and mass transfer coefficients in the following manner: µ int = µ int (x, u int ) for x ∈ Ω h , µ(u) = µ Ω (x, u) for x ∈ Ω and µ(u) = µ Π (x, u) for x ∈ Π. Moreover everywhere in the paper we assume that for any possible solution the interphase coefficients are essentially smaller that those of both adherents.
where u int = [u 1int , u 2int ] defines temperature and mass concentration inside the interphase.As a result, the thermo-diffusive phenomena within the intermediate layer has also essentially non-linear character.Additionally, the intermediate layer may incorporate internal heat and mass sources, described by a given vector function: whose contribution to the problem may also be significant.This can be implemented to the asymptotic procedure by the following assumption: At the external boundary of the combined domain classical types boundary conditions are imposed: where g is a known vector-function, whereas A Γ is one of the classical boundary operators constituting Dirichlet, Neumann or Robin type boundary conditions (see for example [17] ): Here n is a unit outward normal to the boundary Γ, while α and β are some positive constants.One can also consider the so-called mixed boundary conditions when the different operators from ( 9) are defined on different parts of the boundary .Along the transmission surfaces (that are z = h 1 and z = −h 2 ) the so-called ideal or perfect transmission conditions are assumed to be true: Here u ± and q ± are considered limit values of the respective component fields while approaching the interphase from Ω and Π side, respectively.
In the following, the formulated multiphysics problem in the combined, semi-infinite domain containing thin interphase is reduced to an equivalent one defined in the bounded subdomain Ω only (see Fig. 2), with a non-local boundary condition at the transmission surface.Such condition, in turn, incorporates all the basic information about the semi-infinite subdomain Π and the intermediate layer Ω h .
To this end, the transmission conditions between Π and Ω are evaluated in terms of component solutions' jumps across the interphase under assumptions (4) and ( 5) and (7).These conditions together with some operator relation defined previously for subdomain Π [18], are used in 2.3 to formulate the aforementioned non-local, non-classical boundary condition.

INTERMEDIATE LAYER -IMPERFECT TRANSMISSION CONDITIONS
In view of the presence of a small parameter (interphase thickness), direct numerical FEM computation poses a challenging task.On the other hand, this geometrical peculiarity may be successfully utilized to simplify the problem by means of asymptotic analysis.Namely, instead of the main problem, one can consider a sequence of auxiliary problems of simpler structures.Each of them might be solved by either numerical or analytical methods.In such an approach the interphase contribution to the external BVPs is reduced to the imperfect transmission conditions between subdomains Π and Ω, defined along artificial interface.Distribution of the field within the interphase can be derived from the solution to the internal BVPs.
To formulate the aforementioned imperfect transmission conditions one needs to utilize the ideal contact conditions on either side of the intermediate layer (10) together with the assumptions (4), ( 5) and (7).After appropriate rescaling z = εξ, the governing equation within the intermediate layer takes form: in the rescaled domain: Note that if the sources are weak enough we will have X * = 0, what essentially simplifies the analysis (see [14] ).
Everywhere below we neglect the discrepancy on the right-hand side of ( 11), taking into account that it is of two order smaller than the rest of the equation.
In order to determine distribution of the physical fields within the interphase and to evaluate the transmission conditions, let us assume that in each successive iteration temperature and mass concentration fields are known.Thanks to this temporary assumption, one may attribute pertinent material properties to the appropriate points of the domain, what means that: After formal integration of the equation ( 11) we obtain: , .
and a = [a 1 (x, y), a 2 (x, y)] are unknown parameters, which may be defined by substituting ( 13) into (10) 3 and (10) 4 : On the other hand, subtracting (10) 3 from (10) 4 , one can derive the jump of the flux across the interphase: .( 16) Integrating ( 13) within the interval −h * 2 to h * 1 we evaluate the solution jump: where: Note that relations ( 16) and ( 17) constitute together imperfect transmission conditions between subdomains Π and Ω.Unfortunately, they cannot be implemented directly, as, due to ( 18) and ( 12), they require information about temperature and mass concentration distribution within the intermediate layer.Unlike linear case, the non-linear behaviour makes it impossible to derive analytical formulae defining these quantities in a closed form (compare [19]).On the other hand, one may compute the temperature and mass concentration within the interphase using the same asymptotic results.In fact, let us analyse equation ( 13) again, integrating it from −h 2 to an arbitrary z, to receive: Parameter a can be equivalently computed from the ideal transmission condition (10) 1 .Thus, (compare with ( 15)): where matrix E and vector F are described in (18).Finally, substituting (20) into (19) we obtain integral equation, which defines temperature and mass concentration within the intermediate layer:

EQUIVALENT BVP WITH NON-LOCAL BOUNDARY CONDITIONS
Let us assume that along the external part of the half-space S A (see Fig. 2) the Neumann type boundary conditions are imposed.On the assumption that µ ∏ is constant, solution of the respective BVP in the half-space can be found analytically by various techniques (integral transform, Green function e.g.).That allows to obtain respective solution in the form: and, as a result, the following relationship between the solution and its normal derivative along the boundary S T : here G 0 q − (x, y) = Gq − (x, y, 0) is an integral operator acting on the flux q − , whereas function x y h F F s 0 (x, y) = s(x, y, 0) incorporates all the information about given sources and fluxes at the domain boundary.If µ ∏ depends on the possible solution, one may built an iterative process which enables to deal with such cases as well.Substituting ( 22) and ( 16) into ( 17) we obtain a non-local boundary condition, for subdomain Ω, which incorporates all the basic information about Π subdomain as well as about the intermediate layer: .
Note here that the operators E and F depend on the solution within the interphase (see ( 16)) as well as function h: whereas function u int should be found from (21).
In this way, the original BVP defined by equation ( 1) under boundary conditions (8) in the semi-infinite, combined domain Π∪Ω h ∪Ω may be reduced to its asymptotically equivalent form in Ω subdomain only.The equivalent BVP is made up by equation ( 1) under boundary conditions (8) at the external boundaries Γ Ω of the domain Ω, and the non-local, nonclassical boundary condition (24) at the transmission surface S T (see Fig. 2) together with the integral relationship (21).

ALGORITHM OF SOLUTION
As mentioned previously, the equivalent BVP in Ω subdomian is solved by means of Matlab FEM block.Obviously when Γ Ω = ∂Ω, one can solve BVP given by equation (1) and boundary conditions (8) with any FEM commercial package, as all of them allow users to deal with the classical boundary conditions.Unfortunately, while tackling BVP (1), ( 8), (24), the problem does not lend itself to simple solutions in frames of commercial software.For this reason, we construct below a numerical iterative algorithm capable of solving this kind of problems, which enables to implement indirectly the non-local boundary condition (24) into the computational scheme.As it will be shown later, the algorithm has proved to be effective in a wide range of the parameters under consideration.Moreover, such an approach allows us to tackle possible solution singularities near the corner point of the domain, without introducing any special elements into the FEM procedure itself [15], thus without interference into the FEM stiffness matrix.
Note also that we determine distribution of the temperature and mass concentration field within the intermediate layer directly from (21) at each successive iteration.

ITERATIVE PROCEDURE
The proposed iterative algorithm is straightforward.At the first step, we take an arbitrary initial vector (superscript denotes the number of iteration).This function is considered as a first type boundary condition applied at S T instead of condition (24).Thus, in the first step the equivalent BVP incorporates only classical types boundary conditions.Additionally, we assume equation (1) to be linear (constant heat conductivity and mass transfer coefficient).From the solution obtained under these assumptions, heat and mass fluxes at the boundary S T are approximated by means of a

) x y h x y h
x y h dedicated procedure.In the next iterations, the non-local boundary condition (24) is used.Namely, at each successive (i + 1) -iteration we define the non-local boundary condition (24) in a form of the classical one (9) as: where: , , .
The governing equation ( 1) is linearised in the following manner: Distribution of the field components within the interphase, which is necessary for formula (21), is introduced in our computations also by means of iterative procedure, where in each successive iteration the following assumption is accepted (superscripts denote number of iteration): . (29) As a result, the nonlinear problem under consideration is made up by equation (28) under boundary conditions ( 8) and ( 27) with taking into account (29) at each next step.The whole process is reiterated until the normalized difference between two consecutive numerical solutions δ i is smaller than some prescribed value γ.

NUMERICAL RESULTS AND DISCUSSION
In the following numerical examples, a simplified model of domain geometry will be used.Namely, the axial symmetry of the phenomena will be assumed, what allows us to reduce considered problem from 3D to 2D.The domain geometry considered here is shown in Fig. 3.In Fig. 4 the equivalent BVP is presented.
Respective subdomains are defined as: , and .As for the boundary conditions at the external boundaries, we assume them to be of Dirichlet type at S U and S R , while the natural symmetry condition at S S (axial symmetry) gives us a homogeneous Neumann condition.Note that Γ Ω ≡ S S ∪ S U ∪ S R .
Integral operator G 0 from (21) was derived for presented geometry in [18]: where K is a complete elliptic integral of the first kind (see.g. [20]).Moreover, to simplify the problem even more we assume that heat and mass sources within the intermediate layer are described by a function of only r variable (X = X(r)).Then ( 14) is reduced to the form: and consequently (18) (2) may be written as: Knowing these values one can apply the procedure described in 3.1.The algorithm is tested on the basis of analytical benchmark constructed in section 4.1.Section 4.2 is dedicated to a practical problem describing joint thermo-diffusive field.

VERIFICATION OF THE APPLIED PROCEDURE
In order to investigate the convergence rate and accuracy of the procedure, the following analytical benchmark is used.Namely, we assume that material coefficients within the bounded subdomain are constants, and there is no source term in equation ( 1): As a result, the system of partial differential equations (1) splits into two independent Laplace equations in a cylindrical coordinate system.Moreover, we normalize the boundary conditions ( 8) and (24) on the value of µ, so this effectively means that µ 1 = 1 value can be assumed.
As the solution to the problem ( 1)-( 8), (24) we choose the following functions: (34) with the same constants c k , d k for both fields.Note that u 1(2) relates here only to the solution in subdomain Ω. Respective boundary conditions on the boundaries S R , S U are to be taken in accordance with solutions (34).In physical interpretation formulae (34) describe the joint thermo-diffusive field induced by a set of point heat and mass sources, respectively.We assume that d k <0, thus all of them are located outside the domain Ω.
Under such assumptions, the feedback between thermal and diffusive fields is realized only through the intermediate layer, or in other words via the boundary condition (24), where the field components remain joint by the nonlinear boundary conditions.
It is natural that the original boundary conditions (24) do not satisfy solution (34).For this reason we have corrected them slightly, adding a known component to the right-hand side: where R denotes a right hand side of (24), q w is the heat/mass flux, computed in the same way as in (3), while functions E w , h w are determined in the same way as in ( 18), ( 24) on the assumption that the solution within Ω is expressed by (34).The formulae for material constants describing average heat conductivity and mass transfer coefficient within intermediate layer are: , Unfortunately, experimental data describing the heat conductivity and mass transfer coefficient as functions of temperature and mass concentration simultaneously are hardly available.We decided to use relations (36)-(37) as each successive multiplier incorporated in them characterizes material properties in a standard way (for example as for crystalline material [21]).Pertinent constants were taken in such a way to cause the maximal coefficients deviation from their average values less than 20 percent (within predicted range of solution values).All the values are considered normalized, thus dimensionless.Subdomain Ω is defined by: a = 1, b = 0.9.For (34) the following parameters values are assumed: Computations were carried out for the mesh of triangular linear finite elements, composed of 2145 nodes.The convergence rates for component thermal and diffusive problems are shown in Fig. 5a.Here δ 1 and δ 2 denote the L 2 -norms of difference between two consecutive solutions.As the PDEs and the boundary conditions at external part of the boundary are the same for the heat and mass components, it is obvious that the difference in convergence is caused by different character of relations (36) and (37).At first glance one may conclude from Fig. 5a that the accuracy of the computation can be increased to any requested order.However, this is not in the case, as the values of δ 1 and δ 2 only indicate the convergence rate, but not the accuracy itself.To clarify this fact, we have drawn the relative error of the thermal component solution (Fig. 5b) defined as for the 12th iteration.(Error distribution for the diffusive component problem has similar character and order).As .u   anticipated, the maximum error appears near the transmission boundary what is in accordance with the maximum principle for elliptic PDEs.If we consider variation of the maximal relative error values in the bounded domain Ω in each successive iteration, we can see that the accuracy level stabilizes near e 1 = 7.90⋅10 −4 , e 2 = 7.56⋅10 −4 , and cannot be improved with further iterations as it follows from Fig. 5c.It is notable that the maximal relative error of FEM Matlab block itself for the same mesh and boundary conditions is of the same order, and the character of its distribution is also similar to that shown in Fig. 5b.In fact the accuracy limit is reached very quickly so e 1 and e 2 stabilize on their minimal level after sixth iteration, what can be seen in Fig. 5c.Further reiteration of the procedure does not improve the accuracy of solution.Note that other nonlocal boundary conditions taken in form , , provide the same accuracy and character of convergence as (35).

NUMERICAL EXAMPLE
In this section a numerical example of axisymmetrical thermo-diffusive problem (with unknown solution) is considered.We take advantage of the domain geometry defined in the previous chapter.The character of relations describing material coefficients remains the same as in ( 36) and (37).For Ω subdomain we assume: , . (40) Note, that these assumptions introduce a non-linear effect into Ω.Material properties of the intermediate layer are as follows: , Finally, constant values of heat conductivity and mass transfer coefficient are assumed within the subdomain Π: µ 1Π = 1.9, µ 2Π = 5.10 −9 .The values of the heat conductivities were taken in the range specific for building construction materials, separated by a thermoresistant layer [22].In case of diffusion coefficients, the range of their values corresponds to the porous media infiltrated by liquid [23].
As for the boundary conditions in Ω, the following vector functions are imposed (subscript corresponds to the part of the boundary, at which condition is imposed -see Fig. 4 w w e w u w As explained in section 2.3, considered BVP must be supplemented with Neumann type boundary conditions on the part S A of boundary surface (see Fig. 2) of the infinite domain Π For the numerical example we assume that:    (44) Computations were carried out for the same as previously mesh of finite elements.The intermediate layer thickness was assumed 0.01.
As it follows from the previous benchmark example, computation can be stopped when both values δ 1 and δ 2 are less than a maximal possible FEM procedure accuracy, achieved for the same mesh and standard boundary conditions.Thus, results presented in Fig. 6b-6e correspond to the fourth iteration.
However, further reiteration can be treated here as an additional test of the convergence and accuracy.The maximal relative difference between solutions obtained in the fourth and tenth iterations is of the order of 10 −5 , what is less then the accuracy of FEM for the chosen mesh (which is 10 −4 ).
The component solutions are shown in Fig. 6b and Fig. 6c.Respective solutions gradients are presented in Fig. 6d and Fig. 6e.The convergence rate for considered numerical example is shown in Fig. 6a.
In order to investigate the influence of thermo-diffusive feedback on the final result, we consider other problem, analogous to the previous one.Namely, we still assume the same geometry and boundary conditions, but pertinent material coefficients are taken constant.Their values are calculated as average values (arithmetical mean of the nodal values) of heat conductivities and mass transfer coefficients obtained in the previous example: µ 1Ω = 1.18, µ 2Ω = 8.98⋅10 −10 , µ 1int = 1.85⋅10 −2 , µ 2int = 4.60⋅10 −12 .Material coefficients in subdomain Π remain the same.The BVP formulated in this way does not exhibit any relationship between different components of the physical field.Thus, it may be treated as two independent BVPs (thermal and diffusive).Obviously, each of the formulated problems becomes a linear one.The relative differences between respective linear and non-linear problems, calculated as where u L is a solution of linear problem, are shown in Fig. 7a-Fig.7b.Their maximal values amount to: 0.16 for thermal component problem and 0.47 for diffusive component, respectively (note that maximal discrepancies are situated at the transmission surface).From the presented example it follows that the thermo-diffusive feedback, which in our case is the reason for the non-linear character of BVP, has essential influence on the final solution.It is especially clear while considering the diffusive component problem.Obviously, such situation takes place due the strong non-linearity caused by the exponential multiplier in relations (40), (42).Thus, one may conclude that the interaction between respective fields may be neglected only under condition that pertinent material coefficients undergo very little changes within possible solution range.

CONCLUSIONS
It has been shown that proposed method enables to solve problems of joint physical fields in combined, semi-infinite domains.Moreover, such an approach allows to take into account the presence of a thin resistant layer situated between subdomains.Although the initial BVP is reduced to the equivalent one defined only in the bounded subdomain, the physical field distribution within the intermediate layer and the semi-infinite subdomain may also be reconstructed.Namely, one can find it: for interphase directly from (19) and for subdomain Π using representation (20).The convergence rate and accuracy of the built iterative procedure were investigated on the basis of analytical benchmark.It turned out that the accuracy reached within a few first iterations is comparable with that guaranteed by used mesh of finite elements under standard boundary conditions.
In the class of physical problems considered in this paper, we do not encounter possible field singularities.This is a direct consequence of the form of evaluated imperfect transmission conditions.Such singular fields appear near the corner points in case of ideal contact between the different domains.However, the presented approach may be successfully used for such kind of problems as well.It was proved for a single-field phenomena in [12].It is worth mentioning here, that in case of ideal interface when material parameters are constant, the flux singularity exponent can be easily computed analytically.Unfortunately this assertion cannot be extended to non-linear problems with ideal contact.However, in the latter case the exponent can be successfully found numerically.This is not of any surprise that the proposed numerical procedure should have its own limitations.It has been discovered by direct simulation that when the interphase thickness essentially decreases, the procedure does not converge to the solution.The reason for this is quite simple.Let us remind here that the asymptotic procedure used to evaluate the transmission conditions ( 16)- (17) was carried out on the assumption that interphase thickness and its material properties are small parameters of similar order.Thus, the model with the imperfect interface was derived and valid for comparable values of these two small parameters (see (4), (5)).If the geometrical parameter is essentially smaller than respective parameter related to the material properties of the interphase and adherends, then the limit external problem exhibits ideal transmission conditions between the subdomains.This, in turn, immediately results in flux singularity at the corner point.Obviously, such a singularity has not been taken into account in the presented numerical procedure.However, if one corrects the algorithm in a way described in [15], transforming the direct iteration to a projective-iterative method, the computations become stable again.

Figure 1
Figure 1 Original BVP in the infinite combined domain.

Multiphysics Volume 3 Figure 2
Figure 2 Domain decomposition and equivalent BVP.

Figure 5b
Figure 5bRelative error of the thermal component of the solution for the benchmark.

Figure 5c
Figure 5cMaximal relative errors of component solutions for the benchmark.

Figure 6a
Figure 6a Procedure convergence for numerical example.

Figure 6b
Figure 6b Thermal component solution.

Figure 6e
Figure 6d Absolute vale of u 1 gradient.

Figure 7a Figure 7b
Figure 7aRelative difference between thermal component solutions.