Composite materials with thin coating layers (interphases) are encountered in a variety of engineering applications. This explains a large number of publications on the topic of interphase modeling in which various problems were studied in the context of heat conduction, elasticity, thermoelasticity, etc., see comprehensive reviews in, e.g., [1][2][3][4][5][6][7]. As characteristic scales in such problems may vary across many orders of magnitude, direct simulation of coating layers could be, in general, challenging. As an alternative, the so-called imperfect interface approach is available. It reduces a typical three-phase configuration (matrix, reinforcement, and coating) problem to a two-phase configuration (matrix and reinforcement) problem, in which presence of the layer is mimicked using jump conditions in relevant fields across some boundaries.
Here we focus on the so-called asymptotic-based imperfect interface modeling approach in which jump conditions across either a zero thickness interface (typically located at mid-surface of a layer) or the layer are obtained as the result of formal asymptotic analysis of a three-phase configuration problem. The publications relevant to such an approach could be classified into two major groups.
The publications of the first group use the matched asymptotic expansion method. The reviews of earlier literature on the method could be found in [8][9][10]. More recently proposed models of that type are discussed in, e.g., [7,11,12] and in the references therein. In the method, thickness of the layer is chosen as a small (perturbation) parameter and the order of an interface model is defined by highest power of that parameter in the asymptotic series. The method results in a recursive procedure in which resulting jump conditions across the interface are mostly presented in implicit forms and the jump conditions for the zero order models coincide with those across the layer. Most publications on the matched asymptotic expansion method deal with problems involving planar or straight interfaces. We are only aware of a few relevant papers that deal with curved interfaces, see [13,14]. However, it is clear from the plots presented on Figs. 2, 3 of [13] that the reported model only works well for very thin layers and for a rather limited range of material parameters. Even though some authors (e.g., [7,14]) refer to their models as "higher order models," they, in fact, only report jump conditions of zero and first orders.
The second group of publications are based on the use of Taylor series expansions of the fields involved. In that approach, the order of the interface model is defined by the truncation order of Taylor expansions. Significant contributions to the development of the approach are due to Hashin [15,16], Bövik [17,18], and Benveniste [19]. However, in Appendix B of [19], clear distinction was made between the approach of [15,16] and that of [17,18]. In the former papers, the fields related to the layer material are expanded in the normal direction about the points located at its inner boundary and the expansions are used to evaluate the fields at the points located at the outer boundary of the layer. At both boundaries, the perfect bond conditions are enforced to exclude the fields related to the layer. The obtained jump conditions across the layer are treated as the jump conditions across the interface, which coincides with the trace of the inner boundary of the layer.
The approach of Bövik and Benveniste (referred in the literature as the Bövik-Benveniste methodology) was first developed in [17,18] and further generalized in [19]. This approach is based on the following two-step procedure. At the first step, the fields at the points located at the mid-surface of the layer are expanded in the normal direction about the points located at both boundaries and, after enforcing the perfect bond conditions, the jump conditions across the layer are obtained. At the second step, additional Taylor series expansions are used for the jumps of the first step to eventually obtain the expressions for the jump conditions across the interface, which coincides with the trace of the mid-surface. Later on, this methodology was used by several authors, see, e.g., [20][21][22][23][24][25]. Most of those papers reported first order models. For two-dimensional elasticity, the second order model was derived in [26] for interfaces of arbitrary shapes; the third order model was derived for the layers of circular shapes in [27]. For potential problems, the interface model of arbitrary order was constructed for interfaces of arbitrary smooth curvatures in [6] (for two-dimensional problem) and in [28] (for three-dimensional problem). The models were only tested for two-dimensional problems with circular interfaces. The attractive feature of the higher order models obtained with Bövik-Benveniste methodology is formulation of jump conditions in explicit forms, which is convenient for their incorporation in numerical models, see [29].
In [26,[30][31][32] yet another model was suggested and studied. This model consists in the formulation of the jump conditions across the layer using the first step of the Bövik and Benveniste methodology and using them in two-phase configuration problems as jump conditions across the traces of the layer boundaries. The latter model was labeled as Model I and the model obtained at the second step of the methodology was referred to as Model II. While the boundary value problems associated with Model II were somewhat studied, see e.g. [33][34][35][36], we are not aware of similar studies for the boundary value problems associated with Model I. However, for the problems with circular and spherical interfaces for which analytical solutions can be readily constructed, those solutions are unique. Additionally, it was demonstrated in [30,32] that the both models with the jump conditions of the first order satisfy the reciprocity theorem.
The analysis of relevant literature also reveals that, unlike some first order models, see, e.g., [23,24], higher order models have only been tested for problems involving a single reinforcement. Even for such problems, rigorous benchmark solutions are mostly available for two-dimensional cases with circular interfaces.
In the present paper, two types of higher order imperfect interface models for three-dimensional steady-state conduction problems involving coated spherical particles are derived. The paper has the following goals: i) to present arbitrary order jump conditions for imperfect interface models with spherical surfaces, (ii) to incorporate these conditions in the unit cell model of particulate composite for accurate solutions for the local fields and overall properties, (iii) to create reliable benchmark examples for high order imperfect interface modeling of steady-state conduction problems involving spherical interphases.
Consider a steady heat conduction problem involving a heterogeneous solid composed of the following three homogeneous and isotropic materials: (i) infinite matrix of conductivity k 2, (ii) spherical inhomogeneities of radius R 1 and conductivity k 1 , and (iii) interphase layers
The temperature φ and the heat flux q in each phase of the composite system are related by the Fourier law
Unit cell model of the coated spherical particle composite in which the conductivity corresponding to that phase is used. Assume that all conductivities are constant and thus, the temperature φ in each phase satisfies the Laplace equation
In what follows, the classical Rayleigh unit cell model (cube with a single inhomogeneity, Fig. 1) is used as geometry model of a coated spherical particle composite. All the fields related to the inhomogeneity, the interphase layer, and the matrix are identified by the superscripts '1', '0', and '2', respectively. The volume fraction of the compound (core plus interphase layer) inhomogeneities is given by c = 4 3 π (R 2 /a) 3 , where a is the linear size of the unit cell. A single inhomogeneity in an unbounded solid is a specific case of this model corresponding to c = 0.
In the case of macroscopically uniform heat flux, the temperature field inside the matrix of a periodic composite is a superposition of the linear far field φ f ar = G • x and the periodic perturbation field φ per :
where x = x j i j and G = G j i j = grad φ is the macroscopic temperature gradient. This is fulfilled by imposing the constraint equation
The spherical coordinate system (r, θ, ϕ)
is introduced in a way that its origin is located at the center of the inhomogeneity. The inner and outer surfaces of the interphase layer coincide with the coordinate surfaces r = R 1 and r = R 2 , respectively. The following boundary conditions of perfect contact are imposed:
where φ (α) and q (α) = q (α) • n = -k α ∂φ (α) /∂r are respectively the temperature and normal heat flux in the medium 'α'. In what follows, this boundary-value problem is referred to as the Interphase Model, Fig. 2a. Next, consider the model of Fig. 2b, referred to as Model I in [26,[30][31][32]. In that model, the domain occupied by the interphase layer (medium '0') is not directly involved and only accounted for via the jump conditions for the relevant fields across its boundaries. Finally, consider a two-phase configuration problem in which the matrix and core phases are expanded up to the interface R 0 = (R 1 + R 2 ) /2, see Fig. 2c. In this, so-called imperfect interface model of [17] and [19], referred to in [26,[30][31][32] as Model II, the interphase layer is approximated by the jump conditions in the relevant fields across r = R 0 .
In the next Section, the jump conditions for Model I and Model II are derived.
The derivations for the models follow the corresponding steps of the Bövik-Benveniste methodology [17,28]. As it was discussed in the Introduction, the methodology consists of two steps.
In the first step, the jump conditions in the relevant fields across the layer are obtained in such a way that the fields related to medium '0' are eliminated due to the use of perfect bonding conditions, Eq. ( 4). The resulting jump conditions contain the fields related to media '1' and '2' only and they constitute the boundary conditions for the boundary value problems of Model I.
In the second step, the two-phase configuration is considered in which the phases are separated by the interface R 0 = (R 1 + R 2 ) /2. This effectively means that the medium '1' and medium '2' are expanded up to that interface, thus, eliminating the layer (medium '0'). Then, the jumps in the relevant fields obtained in the first step are assumed to be identical to the corresponding jumps across the traces of layer's boundaries in two-phase configuration. After that, the fields involved are let to approach the interface from both sides, resulting in the final jump conditions across the interface r = R 0 . The latter jumps constitute the boundary conditions for the boundary value problems of Model II.
Following [28], we expand φ (0) and q (0) at r = R 0 in a vicinity of the interphase boundary r = R α (α = 1, 2) as
where m α = (-1) α-1 . Importantly, the normal derivatives of φ and q can be expressed via the following surface differential operators
n , and
The explicit expressions for
n , and R (α) n are provided in Appendix A. Here, we note only that they are linear operators with variable coefficients containing only the tangential derivatives and depending on the conductivity k α .
It follows from Eq. ( 5) that
Now, the use of Eqs. ( 5) and the continuity conditions of Eq. ( 4) yields the interrelations between the temperature fields at the surfaces r = R 1 and r = R 2 :
where
The counterpart expression for the normal heat flux is
where
Equations ( 8) and ( 10) are analogous to those derived in [26] in the elasticity context (for N = 2). Those equations allow one to avoid solving for the fields inside the interphase and to carry out the analysis in media '1' and '2' only. It is noteworthy that, by contrast with the classical boundary conditions, they are prescribed at the spherical surface in the bulk rather than at the interface.
Using the expressions for the corresponding surface operators given in Appendix A, one can obtain explicit expressions for the temperature and normal heat flux jump across the interphase
and
Here, ∆ S is the surface Laplacian defined in Appendix C. The expressions for N = 1, 2 are retrieved from ( 12), ( 13) by omitting the terms with the corresponding powers of h, i.e. with h 2 and h 3 for N = 1 and with h 3 for N = 2.
To formulate the jump conditions across the interface r = R 0 , we expand the conditions of Eqs ( 8) and ( 10) into Taylor series in a vicinity of R 0 . This re-expansion procedure is discussed in detail in [28] and re-written here to provide compact form of resulting expressions.
Taylor series expansion of φ (α) | Rα (9) in a vicinity of r = R 0 up to h N has the following form:
The inner sum is appropriately truncated in order to omit the terms with n + m > N .
Similarly, Taylor series expansion of q (α) | Rα of Eq. ( 11) is written as
The normal derivatives in Eqs ( 14) and ( 15) are expressed via surface differential operators using the following differentiation rules:
in which the expressions for
mn , and R
mn (α = 1, 2) will be provided in the next Section.
Using ( 14)-( 17), one case express the imperfect interface conditions of Model II as follows:
× P (α) mn φ (α) + R (α) mn q (α) | R0 = 0. For N = 3, the explicit expressions for the jumps in temperature and normal heat flux across the interface r = R 0 are
The explicit expressions for the jumps related to the lower order models, N = 1, 2, can be recovered from ( 20), ( 21) by dropping the terms with the corresponding powers of h.
For n = 0 and n = 1, the explicit expressions for
n , and R (α) n of Eqs. ( 6) and ( 7) are given by Eqs. ( 46) and (47) of Appendix A, respectively. For n ≥ 2 , they are derived in a recursive manner. It follows from Eq. ( 6) that
Comparison with Eq, ( 6) yields
Another pair of the recursive relations that follow from Eq. ( 7) is
The normal derivatives in Eqs ( 23) and ( 24) are equal to zero for n = 1, see Eq. (46). For n = 2, they obey the following relations:
1 .
For n = 3, the relations
follow directly from Eqs. ( 23) and (24) with n = 2, and at cetera. The explicit expressions for the mixed surface operators of Eqs ( 16) and ( 17) are reduced to those considered above in the specific cases of m = 0
and n = 0:
All remaining mixed operators of Eqs. ( 16) and ( 17) are derived in a recursive manner. For example, term wise differentiation in the right hand side of the identity
yields, in terms of operators,
It also follows from Eq. ( 17) that
Again, the normal derivatives are the only terms in Eqs ( 29) and ( 30) to be determined recursively.
An alternate way to evaluate the mixed surface operators of Eqs. ( 14), (15) consists in using the differentiation rule
where m k are the binomial coefficients, and the analogous formulas for the derivatives
The multipole expansion method [37] provides a complete solution to the boundary value problems associated with the models developed in the previous Section.
The series expansion of the regular far-field temperature φ f ar is
where y s t (x) are the inner spherical solid harmonics (Eq. (57) of Appendix C). The series expansion coefficients c ts are found as
where χ s t are the spherical surface harmonics, Eq. ( 58). In the case of linear andc (q) ts = 0 otherwise. The temperature field φ (1) is regular inside the inhomogeneity. Its series expansion is analogous to Eq. ( 31) and contains the inner harmonics only:
where d ts are the unknown complex numbers. As the temperature is a real quantity, the following relation between the series expansion coefficients must be enforced: d t,-s = (-1) s d ts .
The temperature field φ (0) inside the spherical interphase layer involves a full set of solid harmonics (inner and outer):
where e ts and f ts are the expansion coefficients.
The local expansion of the temperature φ (2) in the matrix part of the unit cell is given by the formulas (e.g., [37])
where a ts are the unknown coefficients. Also,
where η ktls are the triple (lattice) sums defined in such a way as to satisfy Eq. ( 2). For the explicit form of η ktls and more details, see [37]. In the single inhomogeneity case, b ts ≡ 0 and hence φ per (x) → 0 when x → ∞.
Due to the orthogonality property of the spherical surface harmonics χ s t , fulfilling the perfect bond conditions of Eq. ( 4) is a straightforward task. Substitution of φ (0) Eq. ( 34), φ (1) of Eq. ( 33) and φ (2) of Eq. ( 35) into Eq. (4) yields the following set of linear algebraic equations for t ≥ 0 and |s| ≤ t:
Equations ( 37) and ( 36) constitute the infinite linear system from which all the series expansion coefficients can be found with any desirable accuracy. In the single inhomogeneity case, system of Eq. ( 37) is solved independently for each (t, s) pair. In what follows, this solution is used as the benchmark for the approximate models of interphase layer.
The solution protocol for this model is analogous to that discussed in Section .1. The temperature fields inside the inhomogeneity and matrix are given by Eq. ( 33) and (31), respectively. Their substitution into Eqs. ( 8), (10) yields the linear system for the series expansion coefficients, analogous to Eq. ( 37). However, this task is not straightforward because Eqs. ( 8) and (10) involve the surface differential operators. It was mentioned already that these operators are linear ones. For example,
which means that one needs to evaluate L
n and all other operators for χ s t only. Moreover, it appears that
In other words, the spherical surface harmonics are eigenfunctions of the surface differential operators. Application of these operators to χ s t is equivalent to their multiplication by some Prepared using sagej.cls real number (eigenvalue). Equations (38) are remarkable in that they reduce manipulations with surface differential operators to standard algebra. Evaluation of the eigenvalues
nt , and R (α) nt of Eq. ( 38) is based on the use of Eqs. ( 50) and (51) of Appendix A. With aid of these formulas, one obtains
t -
where (α = 1, 2)
The use of Eq. ( 39) and analogous expansion for φ (2) | R2 reduces Eq. ( 8), due to the orthogonality properties of surface harmonics, to the linear system
t -
= L
t -
The second condition, Eq. ( 10), transforms into another set of linear algebraic equations, this time
where
For the explicit form of coefficients for N ≤ 3, see Eq. (55) of Appendix B.
The resolving linear system for the Model II is obtained from Eqs. ( 18) and ( 19) by applying the following formulas for surface spherical harmonics:
The eigenvalues
mnt , and R (α) mnt are found from Eqs. ( 29), (30) using the recursive procedure. Application of Eq. (42) to Eqs. ( 18) and ( 19) again results in Eqs. ( 40) and (41), in which now
The formulas for P (α) t and R
(α) t are obtained by replacing in Eq. ( 43) symbol L with P and Q with R. The explicit expressions for these coefficients for N ≤ 3 are provided in Appendix B, Eq. (56).
The solutions obtained in the preceding Sections provide evaluation of the temperature and heat flux fields in every point inside the unit cell as well as the effective conductivity k ef f of heterogeneous solid. Periodic composite with the unit cell shown in Fig. 1 is macroscopically isotropic and its effective conductivity is given by the Rayleigh's formula [38]
where a 10 is the series expansion coefficient in Eq. ( 35) found from the unit cell problem. Alternatively, one may prefer using the Maxwell homogenization scheme for effective conductivity (e.g., [39]) which reads
where
and a 10 found from the single inhomogeneity problem. The estimates of Eqs (44) and ( 45) are valid for all three models considered in the present paper.
The theoretical solutions to the boundary value problems associated with the developed interface models involve the infinite series expansions and the infinite linear systems for the series expansion coefficients. Its practical implementation consists in numerical inversion of the truncated linear system of finite size and yields an approximate results whose accuracy depends on the number of harmonics retained for numerical study. Applicability of the truncation method is justified for the infinite linear systems with the normal type determinant [40]. In view of the complexity of expressions entering the resolving system, analytical proof of the determinant normality is an overwhelmingly difficult problem. An alternate, commonly accepted approach consists in studying the practical convergence of numerical solution, either in terms of the series expansion coefficients or the elastic fields. In this work, all the numerical data were checked for convergence, and only the practically convergent results are reported in all the plots but Figs 14 and 15, illustrating the convergence rate and the applicability limits of the developed models.
The numerical results presented in this Section illustrate the validity range and accuracy of the developed interface models. The original problem of Fig. 1 involves a number of parameters and its full parametric study is beyond the scope of this paper. Unless otherwise specified explicitly, we consider a relatively thick layer h = 0.1R 2 and keep the following parameters fixed: R 2 = 1, k 2 = 1, and k 2 /k 1 = 10. In all figures, the solid curves represent the practically convergent solution to the three phase configuration problem (Interphase Model), given by Eq. (37). The dashed curves and open symbols represent the approximations of that solution by the first order (squares), second order (circles) and third order (triangles) Models I and II.
In this problem c = 0 and the far field temperature is defined as φ f ar = y s t (x)/y s t (x 0 ), where x 0 = R 0 i 3 . The temperature φ is evaluated at the points x 1 = R 1 i 3 and x 2 = R 2 i 3 of the inner and outer boundaries of the interphase layer, respectively.
Model I By analogy with [28], φ(R 1 ) and φ(R 2 ) is plotted as a function of k 0 /k 1 . These data are shown in Fig. 3(a,b) for t = 1, s = 0 and t = 5, s = 0, respectively. There, the φ(R 1 ) and φ(R 2 ) values for the Interphase Model obtained with the use of (37) are shown by the solid curves 1 and 2, respectively. The open circles represent the values predicted by the first order Model I. As seen, the results for the models considered are visually indistinguishable. Second and third order approximations of Model I practically coincide with the solution to the Interphase Model and therefore are not shown. To better illustrate the performance of Model I, the approximation error δφ (difference between the Model I and Interphase Model) is plotted in Figs. 4,5 as a function of k 0 /k 1 for t = 1, s = 0 and t = 5, s = 0, respectively. As seen from these plots, the relative error is below 1% for all models and tends to zero when k 0 /k 1 → 0 and k 0 /k 1 → ∞. Also, the presented data clearly demonstrate convergence tendency with an increase of the approximation order. In all tests, the first order model provided satisfactory approximation. The second order model is even more accurate and the third order model gives practically convergent solution.
Model II The data shown in Figs 6 and7 are analogous to those in Fig. 3, but obtained using Model II with t = 1, s = 0 and t = 5, s = 0, respectively. As it would be expected, this model is less accurate because of the second Taylor series expansion involved in its derivation. In fact, the first order model (most widely discussed in literature) works more or less satisfactory only for a narrow range of k 0 /k 1 and fails when k 1 /k 2 is small and k 0 /k 1 tends to zero. On the contrary, -5,0 -2,5 0,0
-5,0 -2,5 0,0 the second and, especially, the third order Model II provides reasonably good approximation of the actual potential field. In Fig. 8, the approximation error δφ(R 1 ) of Model II is plotted as a function of interphase thickness h for (a) t = 1, s = 0 and (b) t = 5, s = 0, respectively. These data correspond to k 0 /k 1 = 0.001 where, according to Figs. 6, 7, the approximation error reaches a maximum. It is seen that an accuracy of the first order Model II is satisfactory only for h/R 0 ≤ 0.001 whereas the second and third order models are accurate up to h/R 0 = 0.05 and h/R 0 = 0.10, respectively.
Maxwell scheme for effective conductivity In Figs. 9 and 10, the effective conductivity k ef f of spherical particle composite, predicted by Eq. ( 45) for c = 0.5, is shown for Model I and Model II, respectively. Recall that c is a volume fraction of the coated (core plus interphase) inhomogeneities. Here too, Model I approaches the Interphase Model-based predictions better than Model II. However, the Maxwell scheme is itself the O (c) approximation to the effective conductivity of composite. To get more reliable/accurate data, one has to use the unit cell model.
In this set of tests, an accuracy of the Models I and II in the multi-particle environment is studied.
To get the convergent solution to the unit cell problem, the harmonics up degree t max = 25 were taken into account.
Model I In Fig. 11, φ(R 1 ) and φ(R 2 ) are plotted for c = 0.5 as a function of k 0 /k 1 . Noteworthy, this geometry is close to dense packing of the cubic array of spheres, c max = π/6 ≈ 0.524. The solid curves represent an accurate solution to the Interphase Model, obtained with the use of To investigate them in more detail, we check convergence of the solution for k 0 = 0.0032; h = 0.05R 2 and h = 0.1R 2 ; c = 0.4 and c = 0.5. The results are summarized in Fig. 15. For h = 0.05R 2 , the solution for the first order Model II converges -but more slowly than for the third order model. For c = 0.4 and h = 0.1R 2 , the solution for the first order model is still convergent but very inaccurate. For c = 0.5 and h = 0.1R 2 , the solution for the first order model is divergent and convergence of the solution for the third order model is not reached for t max = 25. So, this is the applicability limit for the third order model. It is believed, however, that this limitation can be released by further increasing the order of Model II.
Rayleigh scheme for effective conductivity In Fig. 16, the effective conductivity k ef f is shown as a function of k 0 /k 1 for c = 0.5. These results are obtained using the Rayleigh homogenization scheme of Eq. ( 44 well whereas the first order Model II greatly underestimate effective conductivity of composite with highly conducting interphase. Comparison of these data with Figs 9 and 10 reveals that Maxwell scheme underestimates k ef f to a large extent.
In this paper, the two-step Bövik-Benveniste methodology is analyzed in the context of steadystate conduction problems involving composite systems with isotropic matrices and isotropic coated spherical particles. Two higher order interface models are derived that allow for the solutions of these problems without a need to explicitly resolve the fields inside the coating layers. In one model, the effect of the layer is accounted for via jumps in the field variables across the traces of its boundaries, while in the other -via corresponding jumps across the trace of its mid-surface. The paper presents several novel contributions to the topic of interface modeling. First, for the two models, compact expressions for the arbitrary order jump conditions are derived and the explicit expressions are presented for the jumps associated with the models up to the third order. It is proven that the spherical surface harmonics are eigenfunctions of all considered surface differential operators which greatly facilitates both analytical and numerical effort by reducing manipulations with surface differential operators to standard algebra. Second, the jump conditions up to the third order are incorporated into the unit cell model of particulate composite to study their performance in a multi-particle environment. Finally, a number of reliable benchmark examples is provided for the problems under study that could be useful for the future investigator.
Based on the results obtained in the paper, few major conclusions can be made. First, comparative analysis of Model I and Model II suggests that Model I leads to more accurate results (both for the local fields and overall properties) than Model II. This could be expected since Model II is derived from Model I using additional truncated series expansions of the fields involved. While the boundary value problems governed by Model I have unique solutions in case of spherical interfaces, it is not clear if that will remain true for interfaces of non-constant curvatures. Second, it is found that the first order Model performs satisfactory only in case of very thin layers and rather narrow intervals of the remaining problems parameters. By contrast, the third order Model II gives the accurate results for a wide range of interphase parameters. Also, the applicability limits of both models are estimated for high-filled spherical particle composite with relatively interphase layers.
The present work can be extended in several directions. One of them is the derivation of the higher order models for elastic problems involving particulate composites with spherical interphases. For both conductivity and elasticity problems, it is also worthwhile to study the influence that the location of the interface has on the performance of Model II. It is clear that using mid-surface is only justified in the case when the core and the matrix materials have the same properties.
For all subsequent n values, (α)
n-1,t .
The explicit formulas for the low order mixed operators entering Eqs. ( 14) and ( 15) are Appendix B. Matrix coefficients of the linear system of Eqs ( 40) and (41)
For the Model I, the matrix coefficients of the linear system of Eqs. ( 40) and ( 41) are (α = 1, 2) : Recall that m 1 = 1 and m 2 = -1.
For the Model II, the matrix coefficients of the linear system of Eqs. ( 40) and (41) are
Prepared using sagej.cls
Journal Title XX(X)