A composite material with inextensible-membrane-type interface

We consider a model of a composite material with “inextensible membrane type” interface conditions. An analytic solution of a stationary heat conduction problem in an unbounded doubly periodic two-dimensional composite whose matrix and inclusions consist of isotropic temperature-dependent materials is given. In the case where the conductive properties of the inclusions are proportional to those of the matrix, the problem is transformed into a fully linear boundary value problem for doubly periodic analytic functions. The solution makes it possible to calculate the average properties over the unit cell and discuss the effective conductivity of the composite. We present numerical examples to indicate some peculiarities of the solution.


Introduction
The importance and applications of composite materials are increasing very fast in the last decades. In part, this is due to their very flexible potentialities to distinguish and use a great set of different physical properties which may be described through different meanings (e.g. within mechanical, thermal and electrical senses).
Thus, it is important to propose, understand and analyse different models of composite materials in view of their potential use in different applications. This can be achieved having in mind the unique character of the detailed microstructure of the composites, but their different properties are also significatively dependent on the boundary conditions which describe their interfaces. Namely, in a composite material we typically have several reinforcements such as particles, flakes and fibers which are embedded in a matrix of metals, polymers or ceramics. Although the shape is determined by the matrix when holding all that items, it is the combination and interaction of all them that determine the overall physical properties of the matrix. There, the boundary conditions play a crucial role. Moreover, it is also natural to allow different conditions on different parts of the boundary. Sometimes, it is also relevant to consider the boundary conditions as transmission conditions in some bounded manifold in the interior of the domain. This allows the proposal materials to have thermal and electrical conduction of very different strengths.
From the mathematical point of view, all that need to be formulated in appropriate function spaces within which the boundary conditions must agree (e.g. by using trace theorems, etc.). So, when proposing a model, it should exist a compromise and interaction between the needs from the applications and the convenient theoretical setting; see, e.g., [1,2,6,7,10,11,12,13,16].
In the case of randomly distributed components, effective properties of such composites were successfully studied, for example, in [3,14,15,17,21,23,26,28], while analytical and numerical results for composites with a periodic structure can be found in [4,6,9,18,24,25]. An extensive and complete overview of the employed methods can be found in the fundamental work [22].
Bearing all this in mind, in the present paper we are proposing a model of a composite material which is also determined by some general boundary conditions (and so incorporating different scenarios). We construct an exact solution for the unbounded doubly periodic nonlinear composite under specific assumptions on material properties of the components. Namely, we consider the static thermal conductivity problem of unbounded 2D anisotropic composite materials with circular non-overlapping inclusions in the square unit periodicity cell geometrically forming a doubly periodic structure. We suppose that each component of the composite is imperfectly embedded in the matrix. Namely, the boundary interface conditions are the so-called "inextensible membrane type". It allows very flexible properties on the interfaces upon the choice of different sequences of fixed parameters on the boundary. The conductivities of the matrix and the inclusions depend on the temperature. The key assumption is that ratios of the component conductivities are independent of the temperature. The external flux is assumed to be arbitrarily oriented with respect to the composite symmetry. We determine the temperature and flux distributions and derive the effective conductivity of such composites.
The spaces where we will consider the problem are convenient to prove the existence and uniqueness of a corresponding solution (upon some conditions). This will allow us to understand also the general properties of such solution and to interpret their most useful properties as it concerns the mechanical and physical behaviour. In particular, this will be also done with the help of some software. In our case, we will be particularly concerned with the description of the effective conductivity tensor of a steady-state heat conduction problem in 2D unbounded doubly periodic composite materials with temperature dependent conductivities.

Formulation of the problem
We consider a lattice in the complex plane (identified with C ∼ = R 2 ), with complex variables being denoted by z := x + ıy (for real numbers x and y). As the representative cell, we take the unit square Let E := m 1 ,m 2 {m 1 + ım 2 } be the set of the lattice points, where m 1 , m 2 ∈ Z.
The cells corresponding to the points of the lattice E will be denoted by It is considered the situation when mutually disjoint disks (i.e., inclusions) of (possible) different radii D k := {z ∈ C : |z − a k | < r k } with the boundaries ∂D k := {z ∈ C : |z − a k | = r k } (for k = 1, 2, . . . , N ) are located inside the cell Q (0,0) and periodically repeated in all the cells Q (m 1 ,m 2 ) . Let us denote by We investigate the entire infinite composite where the matrix and inclusions occupy domains with thermal conductivities λ m = λ m (T ) and λ k = λ k (T ), respectively. Here, the temperature T is defined in the whole R 2 . We assume that the conductivities λ m , λ k (k = 1, . . . , N ) are continuous bounded positive functions on R.
We search for the steady-state distribution of the temperature and heat flux within such a composite. The problem is equivalent to determining the function T = T (x, y) satisfying the nonlinear differential equations We assume that the non-ideal contact conditions on the boundaries between the matrix and inclusions are satisfied: (2.4) where γ k is a given function (the so-called resistance coefficient), the vector n = (n 1 , n 2 ) is the outward unit normal vector to ∂D k , the vector s is the outward unit tangent vector to ∂D k , and Let us mention that the usually adopted ideal contact conditions consist in demanding the continuity of the temperature and the thermal flux. Here, we use a relaxation of one of these conditions and allow certain discontinuities. Namely, in (2.3), according to the Fourier's law, we assume that the thermal flux jump across the boundary is proportional to the thermal flux of an inclusion within the tangent direction, caused by heat flow around the inclusions.
It is worth mentioning that our boundary value problem can be used for the characterization of other physical or mechanical processes. We meet analogous boundary contact conditions in classical problems of solid mechanics for elastic media known as "inextensible membrane type" conditions. For more details, we refer to [5] where different types of boundary contact conditions are described.
We assume that the average flux vector of intensity A is directed at an angle θ to axis Ox (see Fig. 1) which does not coincide, in general, with the orientation of the periodic cell. This gives the following conditions Note that, in general, the flux is not periodic. However, since there are no sources and sinks in the composite, the energy conservation law dictates This, in turns, allows us to replace conditions (2.5) and (2.6) with those defined on the opposite sides of the cell.

Reformulation of the problem
To solve the problem, we use the Kirchhoff transformation (cf. [20]) and introduce new continuous functions f m and f k (k = 1, . . . , N ) Then, using the representations (3.1) and changing the dependent variables as The functions f m and f k are monotonic increasing functions of temperature and, therefore, there exist their inverses f −1 m and f −1 k . The contact conditions (2.3) and (2.4) can be rewritten now as follows: where the functions are defined for all ξ ∈ R. Note that, in general, the functions u m and u k may have different values on the interface ∂D k . The derivative of F k can be computed as follows where ξ = f k (T k ). If we suppose that λm(T k ) λ k (T k ) and γ k (T k ) λ k (T k ) are proportional nonlinear coefficients, namely, then all functions F k are linear: for any closed curve Γ in the matrix. Moreover, since there is no source (sink) inside the composite (neither in the matrix nor in any inclusion), the same condition is satisfied for any closed simply connected curve within the inclusion (3.14) Finally, the conditions (2.5) and (2.6) are transformed into the following: Let us introduce new harmonic functions inside the inclusions. Then, the transmission conditions (3.5) and (3.6) become 4 Solution of the problem using the same approach as in [9]. We shall now outline some basic facts which we apply.
Let us introduce complex potentials ϕ(z) and ϕ k (z) which are analytic in D 0 and D k , and continuously differentiable in the closures of D 0 and D k , respectively, by using the following relations where B is an unknown constant belonging to C. Besides, we assume that the real part of ϕ is doubly periodic in D 0 , i.e., It is shown in [18] that ϕ is a single-valued function in D matrix . The harmonic conjugate to u is a function v which has the following form: with the same unknown constant B.
For the determination of the flux ∇u(x, y), we introduce the derivatives of the complex potentials: Applying the Cauchy-Riemann equations ∂um ∂n = ∂vm ∂s , ∂ũ k ∂n = ∂ṽ k ∂s the equality (2.4) can be written as Integrating the last equality on s, we arrive at the relation where C is an arbitrary constant. We put C = 0 since the imaginary part of the function ϕ is determined up to an additive constant which does not impact on the form of u. Using (4.2), we have Im ϕ(t) = −Im Bt + 2 c k + 1 Im ϕ k (t) + 2d k c k + 1 Re ϕ k (t). (4.6) Using now (4.1), we are able to write the equality (2.3) in the following form: Adding the relation (4.7) and (4.6) multiplied by ı, and using Re ϕ k = ϕ k +ϕ k 2 , Im ϕ k = ϕ k −ϕ k 2ı , t − a k = r 2 k t−a k , we rewrite the conditions (2.3) and (2.4) in terms of the complex potentials ϕ(z) and ϕ k (z): Representing the function ϕ k in the form ϕ k (z) = ∞ l=0 α k (z−a k ) l , |z−a k | ≤ r k , and by using the relation t = (4.10) Thus, after differentiating (4.8), we arrive at the following R-linear boundary value problem on each contour |t − a k | = r k : with the unknown constant B = B 1 + ıB 2 . An algorithm for solving such R-linear boundary value problem is developed and described in detail in [18]. We use this approach in our computations.

Effective conductivity
We assume that the effective conductivity tensor Λ e depends on the average temperature T and is defined in the following way: where R e = Λ −1 e is the effective resistance tensor. A similar definition to (5.1) has been used in [27]. Here, the operator · is defined as Note that definition (5.1) needs further justification as the question arises whether the approach is invariant with respect to the averaging cell. We will discuss this issue later during the computations.
We represent all elements involved in (5.1) in terms of solutions u m and u k of the problem (3.3)-(3.6). For the total flux in the x-direction, we have Using the first Green's formula x ∂ũ k ∂s ds, and similarly According to the Cauchy-Riemann condition ∂ũ k ∂s = − ∂ṽ k ∂n , we get x ∂ṽ k ∂n ds Thus, where the components q x and q y can be found via a solution ψ k (cf. also the notations (4.3)) of the R-linear conjugation problem (4.11) as Or, using the mean value theorem for harmonic functions, we get Due to Gauss-Ostrogradsky formula and the boundary condition (2.3), the components of the term ∇T in (5.1) are defined as T m (s) cos(n s , e i ) ds where n s and n k s are the outward unit normal vectors to ∂D 0 + m 1 + ım 2 and ∂D k + m 1 + ım 2 , respectively, and e i is the basis vector. Analogously, y)) cos(n s , e j ) ds.
Finally, the average temperature is given by T (x, y) dxdy Note that it is more convenient to first compute the components of the effective resistance tensor R e from the second formula in (5.1) and then find the effective conductivity tensor Λ e = R −1 e .

Numerical analysis
The algorithm mentioned above is realized in the Maple 14 software. In our computations we consider a composite where four inclusions are situated inside the cell Q (0,0) with the centers (defined in the notations of complex variables): a 1 = −0.18 + 0.2ı, a 2 = 0.33 − 0.34ı, a 3 = 0.33 + 0.35ı, a 4 = −0.18 − 0.2ı. The radii of the inclusions are the same r k = R = 0.145 (cf. Fig. 2). Thus, the volume fraction of the inclusions for such composite is ν = 4πR 2 = 0.2642. For this choice, the inclusion boundaries are situated very close to each other (the minimal distance is 0.02).
Note that the following statement is true. Theorem 6.1 Let λ m , λ k and γ k be periodic functions with the same period τ satisfying (3.9). Then the effective conductivity tensor Λ n defined in (5.1) is also a periodic matrix function with the same period τ .
Since the tensor Λ n is uniquely defined in (5.1), we get for any T n . 2 Note that in the linear case the temperature is defined up to an arbitrary additive constant, and this constant is not involved in the determination of the effective conductivity of a composite material. Generally speaking, this is not the case for nonlinear problems, and the additive constant, appearing during the stage of solving the linear problem (3.3)-(3.16), influences on the computation of the effective conductivity tensor of the equivalent nonlinear composite.
Here, we evaluate the effective resistance tensor R e by the interpolation method and then we get Λ e = R −1 e . Two equivalent procedures are used for obtaining a discrete data set of the effective resistance tensor components R x e , R y e , R xy e , R yx e suggested in [19]. (i) First, one can solve the corresponding linear boundary value problem in a doubly periodic domain preserving its uniqueness by any appropriately chosen condition (for example, here we impose that the function u = u * satisfies the condition u * (0) = 0). Then, to evaluate the properties of the composite material, one can compute the average temperature and the effective resistivity for each particular unit cell presenting the data as the functional relationship R e = R e ( T ).  Tracing cells belonging to the strip, we get a set of points which are not dense enough at the top of the sinusoid considered in Theorem 6.1. Therefore, we used a second method.
(ii) Namely, we consider an arbitrary cell in the original domain and build a set of solutions to the linear boundary value problem in the form u = u * + C, where C is an arbitrary constant. Then, for every constant C, the components of the effective resistance tensor, R e , and the average temperature, T , are functions of the parameter C. Changing the value of C continuously from −∞ to ∞, one receives the sought for the effective conductivity tensor of the composite as a continuous function of the average temperature. It is clear that this procedure does not depend on the chosen cell.
Naturally, for the conductivities of the composite components analyzed in this example and for the fact proved in Theorem 6.1, the nonlinear character of the relationship is observed only within the finite interval of the parameter C.
Note that both methods allow one to determine two components of the effective resistance tensor R e for each one of the two orthogonal flux directions. Thus considering θ = 0, we define R x e = R x e ( T ) and R yx e = R yx e ( T ), and choosing θ = π/2 we find R y e = R y e ( T ) and R xy e = R xy e ( T ). As a result, the entire tensor R e ( T ) is defined.
In view of having a more accurate interpolation procedure, we use both suggested methods and calculate 141 data (R x e , R yx e , R xy e , R y e ) with respect to the average temperature T ∈ [−10, 10]. For the chosen configuration it guarantees a computational error less than 10 −6 . The components Λ x e and Λ y e of the effective conductivity tensor Λ e = R −1 e are presented in Fig. 4 and Fig. 5, respectively, for different d k = 0; 2; 10.