1 Introduction
Nanofibers embedded in polymer matrices have attracted attention as one of the reinforcements for composite materials. Carbon nanotubes (CNTs) reinforced polymer nanocomposites are considered as conventional micro and macrocomposites [15]. Their thermal, mechanical, and electric properties are determined by experimental and theoretical investigations [3, 6, 14]. CNTs are considered as perfectly conducting inclusions, which suggests imposing Dirichlet boundary conditions on the boundary of CNTs. On the other hand, the classical problems for materials with holes in porous media and materials with voids and insulting inclusions are modeled by the Neumann boundary condition [1, 4].
The present paper is devoted to the heat conduction within a 2D (twodimensional) thermal isotropic and homogeneous nanocomposite, which takes the form of an infinite strip, when it is reinforced by nonoverlapping and randomly distributed CNTs and contains defects and voids. In particular, we are interested in studying the effect of CNTs as well as of the presence of voids on the macroscopic conductive and mechanical properties of this composite. Owing to the superconductivity of CNTs and the extremely low conductivity of voids, we can assume that the conductivity of CNTs, of the polymer host and of voids to be governed by the inequalities . Such an assumption leads to a mixed boundary value problem where the latter inequality becomes . The host conductivity can be normalized to unity, i.e., .
Theoretical investigation of mixed boundary value problems by integral equations can be found in [13, 12]. In the same time, implementation of numerical methods for large number of inclusions and holes is still a challenging problem of applied and computational mathematics. We propose in this work a fast and effective algorithm for the numerical solution of the formulated mixed boundary value problem. The method is based on the boundary integral equation with the generalized Neumann kernel. The integral equation has been used in [12] to solve a similar mixed boundary value problem related to the capacity of generalized condensers. The proposed method can be even employed when the number of perfectly conducting inclusions and holes is very large.
As a result of simulations, we first study the 2D local fields for three types of media. In the first type, we consider the case of pure void cracks with . The second type consists of pure CNT inclusions with and . Finally, we treat the case of a large number of combined inclusions and holes by considering either of one of the two or of each. Afterward, we take up the systematic investigation of the effective conductivity of the considered composites. It is important in applications to predict the macroscopic properties of composites which depend on the concentration of perfectly conducting CNTs as well as on the concentration of holes and voids. It is worth noting that the notation of concentration are different for slit shapes of CNTs and circular shapes of holes. The performed simulations of local fields and computation of their averaged conductivities for various concentrations allows to establish the dependence of the macroscopic conductivity on the main geometrical parameters.
2 Problem formulation
Let us consider a channel medium embedding inhomogeneities in the form of nanofillers and holes (voids). As many nanofillers (e.g, carbon nanontubes [7]) have cross sections of elliptical shapes, we model them as ellipses . Furthermore, we represent the nonconducting holes by inner circles . The top and bottom infinite walls of the channel are denoted respectively by and , which yields a multiply connected domain of connectivity with a boundary set where . An example of this domain for the case of and is illustrated in Figure 1.
The medium matrix without inhomogeneities is supposed to be homogeneous and isotropic with a constant thermal conductivity . We also assume that conduction is the only dominating mechanism of heat transfer in the medium. Except being nonoverlapping, no other restriction is imposed on the inhomogeneities as they can be placed at random orientation and position.
The nanofillers are treated as heat superconductors with an almost uniform temperature distribution within each one. Therefore the temperature is assumed to be fixed to an indeterminate constant value along each ellipse for . This assumption is consistent with the numerical results reported in [17] for CNT reinforced polymer composites. Furthermore, by the law of energy conservation in steadystate heat conduction, there should be no net thermal flow through each nanofiller. This constraint is written by means of the net heat flux boundary condition (1e).
On the other hand, the curves are assumed to be perfect insulators and therefore they act as barriers to heat flow. Henceforth, the Neumann boundary condition (1f) is imposed along the holes contours. Finally, isothermal conditions are imposed on the external boundaries by assuming that the lower infinite wall is a heater of temperature , and the upper wall acts as a heat sink, which can be held at a fixed temperature . Thee two values and of the temperature on the external boundaries are normalized to and , respectively.
Under steadystate conditions, Fourier’s law of heat conduction and the above specified heat boundaries conditions yield the temperature distribution governed by the following mixed DirichletNeumann boundary value problem:
(1a)  
(1b)  
(1c)  
(1d)  
(1e)  
(1f) 
where denotes the normal derivative of , and are undetermined real constants that need to be found alongside the distribution temperature .
3 The integral equation method
The boundary integral equation with the generalized Neumann kernel is not directly applicable to the above boundary value problem (1) because of the external boundary component. However, the boundary value problem (1) is invariant under conformal mapping. The mapping function
conformally maps the unit disk onto the infinite strip . Thus, the inverse mapping
conformally maps the infinite strip onto the unit disk , the real axis onto the lower half of the unit circle, the line onto the upper half of the unit circle, and satisfies . Consequently, the function maps the multiply connected domain in the plane (the physical domain) onto a multiply connected domain in the plane interior of the unit circle and exterior of smooth Jordan curves (the computational domain). In Figure 2, we display the result of the conformal mapping of the example shown in Figure 1.
It follows that the harmonic function can be written as
in which the function is the solution of the following boundary value problem in the plane:
(2a)  
(2b)  
(2c)  
(2d)  
(2e)  
(2f) 
where , , and for . Note that the restriction of the function on the external boundary is discontinuous at . However, the function can be cast into the form
where is a harmonic function in , and
The function is harmonic in with on the upper half of the unit circle and on the lower part. The function is the solution of the boundary value problem
(3a)  
(3b)  
(3c)  
(3d)  
(3e) 
where is the unit circle.
For the orientation of the boundary components of , we assume that is oriented counterclockwise and the other curves are oriented clockwise. We assume that each boundary component , , is parametrized by a periodic function , such that . Let be the disjoint union of the intervals , the whole boundary is parametrized by the complex function defined on by [16, 11]
Note that the unit circle is parametrized by , .
Let
be the unit outward normal vector at
and let be the angle between the normal vector and the positive real axis. Then, for ,(4) 
Thus
(5) 
The harmonic function is the real part of a singlevalued analytic function , i.e., , where
(6) 
and the branch of the logarithm function is chosen such that . Then by the CauchyRiemann equations, we have
which, in view of (4) and (5), implies that
(7) 
Since
it follows that for and ,
(8) 
The harmonic function can be assumed to be a real part of an analytic function , . The boundary conditions (3b) and (3c) give the real parts of the function on for . Specifically, we have
(9) 
and
(10) 
For the remaining boundary components for , we use the condition (3e) to determine the boundary condition on . By the CauchyRiemann equations, we can show using the same arguments as in (7) that
(11) 
Thus, for for , it follows from (3e), (8), and (11) that
Integrating with respect to the parameter for , , we obtain
(12) 
where the integration constants are undetermined. The constants , in (10) and (12) are determined so that is a singlevalued analytic function.
Since we are interested in the function , the real part of , we may assume that is real for some given point in . Define an analytic function in the domain through
(13) 
Define also
(14) 
where is the piecewise constant function given by
(15) 
Thus
which implies that
On the basis of the conditions (9), (10), and (12), the function satisfies the RiemannHilbert problem
(16) 
where
(17) 
It is clear that the function is known and the piecewise constant function is unknown and should be determined. Let , i.e., the boundary values of an analytic function are given by
(18) 
Thus, in order to find the boundary values of the analytic function , we need to determine the two unknown functions and . These two functions can be computed using the boundary integral equation with the generalized Neumann kernel [16, 11, 10].
Let be the space of all real Hölder continuous functions on , let be the identity operator, and let the integral operators and are defined on by
The kernel of the operator is known as the generalized Neumann kernel. For more details, see [16, 11, 10]. On account of [10], we have is the unique solution of the integral equation
(19) 
Additionally, the piecewise constant function is given by
(20) 
We compute approximations to the functions in (19) and in (20) by the MATLAB function fbie
from [11]. This function employs a discretization of the integral equation (19) by the Nyström method using the trapezoidal rule [2] to obtain an algebraic linear system of size where is the number of discretization points in each boundary component. The resulting system is solved
by applying the generalized minimal residual method through
the MATLAB function . The matrixvector multiplication in is computed using the MATLAB function from the MATLAB toolbox [5]. The values of the other parameters in the function fbie
are chosen as in [9]. For more details, we refer the reader to [11].
4 Computing the temperature distribution and the heat flux
By computing and , we obtain the boundary values of the function through (18).
The values of the function for can be computed by the Cauchy integral formula. For the numerical computation of for , we use the MATLAB function fcau
from [11].
Then, the values of can be computed by (13) and hence the values of the solution of the boundary value problem (2) is given for by
We deduce the values of the temperature distribution for any by
Moreover, by computing the piecewise constant function , we can compute as well the values of the undetermined real constants from (16).
The function is the real part of the function
According to the CauchyRiemann equations, it follows that the derivative of the complex potential on is given by
One the other hand,
(21) 
where the denominator does not vanish in the domain since is a conformal mapping. Therefore the heat flux can be expressed for in terms of by the formula
(22) 
Hence
(23) 
The derivatives and in (21) can be computed analytically. So, the values of the heat flux
can be estimated on the domain
by first approximating the derivatives of the boundary values of the analytic function on each boundary components. This can be done by approximating the functionusing trigonometric interpolating polynomials then differentiating. The values of
, in the righthand side of (21), can be then computed for using the Cauchy integral formula.5 Computing the effective thermal conductivity
The medium matrix without inhomogeneities is assumed to be homogeneous and isotropic. We will assume that the CNTs and the circular voids are in the part of the domain between and . Thus, the effective conductivity of a layer in the direction is calculated by the formula (3.2.33) from the book [4, p. 53], which in our case on account of (23) becomes
(24) 
Since
where for , is on the line and for , is on the real line. Thus, for , we have
Since and where
(25) 
Consequently, (24) can be written as
(26) 
In combining (21) with the fact that , we can see that
Hence,
(27) 
which implies that
(28) 
The second term in the righthand side of (27) does not depend on the CNTs or the voids. In view of (6) and (25), we have
and hence
(29) 
Since and are on the unit circle , the external boundary of , and taking into account (13), (14), (15), and (18), Equation (29) can be written as
(30) 
By solving the integral equation (19), we obtain approximate values of at the discretization points. These values are employed to interpolate the approximate solution on by a trigonometric interpolation polynomial, which is then used to approximate the values of and .
6 Numerical results
The above proposed method with is applied to compute the temperature field and the heat flux for several examples. We will choose the CNTs and the circular voids within the part of the domain between and . To compute the values of the temperature distribution and the heat flux , we discretize part of the domain , namely for and . Afterwards, we compute the values of the temperature distribution and the heat flux at these points as described in Section 4.
6.1 The domain with only circular voids
In this subsection, we consider the domain with nonoverlapping circular holes and without any CNT (i.e., ). We also assume that all circular holes have the same radius with the parametrization
where are the centers of the circular holes. As these circular holes are chosen in the part of the domain between and , we define the concentration of these voids to be the area of these circular holes over the area of the rectangle , i.e.,
(31) 
The ClausiusMossotti approximation (CMA) also known as Maxwel’s formula can be applied for dilute composites when the concentration (31) is sufficiently small. Below, we write this formula for a macroscopically isotropic media with insulators of identical circular holes within the precisely established precision in [8]
(32) 
Example 1
We consider circular holes with the radius for . For Case I, we assume the centers of the holes to be set to , , , , and . The contour plot of and for are shown in Figure 3 (first row). The approximate value of the effective thermal conductivity for is
When is close to , the circular holes become adjacent to each other. To show the effects of the radius on the effective thermal conductivity , we compute the values of for several values of , . The obtained results are presented in Figure 4 where, by (31), the concentration of these holes is for and . The values of the estimated effective conductivity is given also in Figure 4. As one can expect, there is a good agreement between and for small values of . In the same time, the divergence of and is observed for the concentrations greater than .
For Case II, the centers of the holes become , , , , and , which means the centers are not anymore horizontally aligned as the second and fourth centers are now shifted by 0.2 up and down, respectively. This is displayed in Figure 3 (second row). The curve showing the obtained values of as a function of the concentration is depicted in Figure 4.
Figure 4 illustrates that the values of depend on the position of the circular holes centers while the values of are the same for both cases since it depends only on the concentration of the circular holes and not on their positions. We notice a better agreement between and in Cases II when comparing to Case I.
Example 2
We consider circular holes with centers , , and , where for , and with radius for . The contour plot of and for are shown in Figure 5. The approximate value of the effective thermal conductivity for is
Example 3
We take up here the case of circular holes with centers , , , , and , where for , and with radius for . On the basis of (31), the concentration of these holes is . For , we have . When is close to , the circular holes become adjacent to each other, and the concentration is almost equal to . The obtained results showing the behavior of as a function of the radius , for , are presented in Figure 6 (right).
6.2 The domain with only CNTs
In this subsection, we consider the domain with nonoverlapping elliptic CNTs without any circular holes (i.e., ). We assume that all CNTs have equal sizes and are of elliptic shape where the ellipses have the parametrization
(33) 
where is the center of the ellipse, and are the length of the ellipses axes in the and directions, respectively. If , the major axis of the ellipses is horizontal, if , the major axis of the ellipses is vertical, and if , the ellipses reduced to circles. Here, we choose and such that their ratio satisfies . These elliptic shape CNTs are chosen in the part of the domain between and . So, we define the concentration of these CNTs to be
(34) 
If , instead of (34) the plane slits density is considered in the theory of composites and porous media
(35) 
For a macroscopically isotropic media with only perfectly conducting identical circular inclusions (CNTs), an approximation of the effective conductivity is given by the inverse to (32) value (see [8])
(36) 
Example 4
We consider elliptic CNTs with centers , , , , and , and where and . Figure 7 (first row) presents the contour plot of and for and (the ellipses are horizontal). For these values of and , the approximate value of the effective thermal conductivity is
For and , the contour plot of and are shown in Figure 7 (second row). The approximate value of the effective thermal conductivity for these values of and is
The CNTs in Figure 7 have the same concentration. However, the value of is larger for the vertical ellipses case.
When approaches , the ellipses get adjacent to each other. On the other hand, they come close to the upper and lower walls when approaches . We compute the values of for several values of , , with . The obtained results are presented in Figure 8 (left) where, by (34), the concentration of these ellipses is . Note that, for and , we have .
Example 5
We consider elliptic CNTs with centers for and where and , and with and .
We compute the values of for several values of , , and (i.e., the ellipses are horizontal) where the ellipses become close to each other when approaches . The obtained results are presented in Figure 9 (left) where the concentration of these ellipses is . For and , we have . Then, we compute the values of for several values of , , and . The obtained results for versus the the plane slits density are presented in Figure 9 (right) where for and .
6.3 The domain with CNTs and/or circular voids
We are concerned in this section with the study of a large number of perfect conductors and/or insulators. We consider two example where in the first both perfect conductors and insulators have the same circular shape, while in the second, conductors have an elliptic shape and insulators have a circular shape. The present investigation is useful when studying the impact of geometric shapes on the macroscopic properties of threephases high contrast media.
Example 6
We take circular holes of equal size with radius . In this example, the concentration is constant and the locations of these holes are chosen randomly. In this case, the following extension of CMA may be used
(37) 
where denotes the conductor concentration, the insulator concentration, and . Three cases are considered:
 Case I:

We assume that half of the holes are CNTs and the other half are voids (see Figure 11). For this case, and are given by .
 Case II:

All holes are voids, and hence while .
 Case III:

All holes are CNTs, and hence while .
For each case, we run the code for 20 times, so that to get 20 different locations for these circular holes. In each of these 20 experiments, we compute the value of the effective thermal conductivity by the presented method and the values of the estimated effective conductivity by (32). As we can see from Figure 12, is a constant and the values of depend on the locations of the holes.
Example 7
We consider elliptic and circular holes with elliptic perfect conductors and circular insulators of equal area (see Figure 13). The radius is chosen to be the same as in the previous example, i.e., . The locations of both elliptic and circular holes are chosen randomly. For the ellipses, we assume that the ratio between the length of the major axis and the minor axis is , and the angles between the major axis and the axis are chosen randomly. As in the previous example, we run the code for 20 times. In each of these 20 experiments, we compute the value of the effective thermal conductivity by the presented method. The computed values are shown in Figure 13 (left).
Since we have the same number of elliptic perfect conductors and circular insulators of equal area , the conductor concentration and the insulator concentration are equal and given by
Comments
There are no comments yet.