1 Conformally mapping the physical domain
Consider the unbounded multiply connected domain being exterior to the rectilinear slits (see Figure 2 (right)). We assume each slit make an angle with the positive axis for . The iterative method presented in [15] computes an unbounded preimage domain in the exterior of ellipses (see Figure 2 (left)). The ellipses are parameterized by
(2) 
where is the ratio between minor and major axes’ lengths of the ellipses. The parameters and , , are computed by the iterative method. The method provides also the boundary values of the unique conformal mapping from the domain in plane onto the domain in plane with the normalization
(3) 
Putting , the domain is on the left of because the curves are clockwise oriented in view of (2). According to (3) the mapping function is analytic in with . Thus, the Cauchy integral formula implies that the values of the mapping function can be obtained for by
We can also use the Cauchy integral formula to determine the values of the inverse mapping function for (see e.g., [24]). The inverse mapping function satisfies also the normalization [23, p. 114, p. 127]
Therefore the values of for can be computed by
(4) 
For the parameterization of the boundary , we adopt the same notations used in [8, 12, 13, 14, 15, 22]. Let
where and means for . In this way, we define a parameterization of the boundary on by
We will drop the index in the notation of as the value of such that will be clear from the context. Therefore we simply write
The iterative method in [15] provides us with the boundary values , which are used to get a parameterization of the boundary of
Thus, by making use of (4), we can determine the values of the inverse map for through
(5) 
To compute the values of , we first approximate the real and imaginary parts of on each interval ,
, by trigonometric interpolating polynomials, and then differentiate.
Now we are ready to compute the preimage of the given domain (see Figure 3 (right)). The inverse mapping function maps the domain outside the slits onto the domain outside of the ellipses . Furthermore, maps the squares and onto piecewise smooth Jordan curves and , respectively, such that the ellipses are in the ring domain between the curves and (see Figure 3 (left)). Thus the whole boundary of the domain is .
For each , we parameterize the square by , , where is chosen using the same approach used in [8, pp. 696697]. Hence for , the curve is parameterized by
such that the values of are computed through (1). The values of the derivative , for , are computed numerically again by using trigonometric interpolating polynomials as explained above. Recall that the parameterization of the ellipses are given by (2). Henceforth if we take as the disjoint union of the intervals , then we parameterize the whole boundary by
2 Computing the temperature distribution
Note first that the boundary value problem (1) is invariant under conformal mapping. Indeed, the unique solution of the problem (1) is given by
(6) 
such that is the unique solution of the following boundary value problem on the domain in plane,
(7a)  
(7b)  
(7c)  
(7d)  
(7e) 
As is harmonic, we can write
(8) 
where is a singlevalued analytic function on the domain . Assuming the boundary values of the function takes the form
on the boundary of , it follows from (7b) that
(9) 
Similarly, equation (7e) implies that the values of the function are known for with
(10) 
No information can be found from (7) about the values of the function for . Nevertheless, the condition (7d) yields for [9, 16]. Thus, is a constant function on the interval , i.e.,
(11) 
where is an undetermined real constant.
Let us define the function on by
Following (9), (10) and (11), the function satisfies the boundary conditions
(12) 
which are known as the RiemannHilbert problem. To solve the problem (12), we tale an auxiliary given point in the domain . Since we are only interested in determining the real function , we can assume that is real. We define an auxiliary function on the domain by
(13) 
The function is clearly analytic in . We consider also a complexvalued function defined on by
where is the piecewise constant function given by
Therefore the boundary conditions in (12) implies that solves the RiemannHilbert problem
given that
where for , , and . The piecewise constant function is an unknown and need to be determined together with the function . Both of these two functions can be found using a boundary integral equation method based on the generalized Neumann kernel. More precisely, following [12] the boundary values of the analytic function are determined through
where uniquely solves the integral equation
(14) 
and the function is computed by
(15) 
Here, is the identity operator and the integral operators and are respectively defined by
We compute approximations to the functions in (14) and in (15) by the MATLAB function fbie
from [14] as follows
[mu,h] = fbie(et,etp,A,gamk,n,iprec,restart,tol,maxit).
This functions employs a discretization of the integral equation (14) by the Nyström method using the trapezoidal rule to obtain an algebraic linear system [4, 20]. Each boundary component is discretized by nodes yielding a linear system of size . This system is solved by the MATLAB function
in which the matrixvector multiplication is computed using the MATLAB function
from toolbox [7]. In our numerical calculations, we choose ; which means the tolerance in the FMM is . The GMRES performs a maximum number of iterations maxit=100 with an accuracy tolerance tol=1e14, and the method is used without restart by choosing restart=[ ].By computing the real functions and , we obtain the boundary values of the analytic function through (2). Then, in view of (13), we obtain the boundary values of the analytic function by
The values of the function can be computed for by the Cauchy integral formula. Finally, following (6) and (8), the values of the temperature distribution is determined for by
(16) 
where
(17) 
is an analytic function on the medium called complex distribution temperature. The derivative if this function is usefully related to the heat flux as we will see in the next section.
3 Computing the heat flux
According to Fourier’s law of conduction the heat flux vector related to the gradient temperature is given by the equation
where is the equivalent thermal conductivity which we assume to be normalized to unity.
From (16) and using the CauchyRiemann equations, it follows that the derivative of the complex temperature distribution on is given by
One the other hand (17) yields
where the denominator does not vanish in the domain since is a conformal mapping. Therefore the heat flux can be expressed in terms of by the formula
(18) 
The values of the heat flux can be estimated on the domain by first computing the derivatives of the boundary values of the analytic functions and on each boundary components using trigonometric interpolating polynomials as explained above. Then, the values of and , in the righthand side of (18), are computed for using the Cauchy integral formula.
4 Examples
We apply the method presented above to compute the temperature field and the heat flux for four different number of CNTs. In the first example, we consider the case of CNTs of different lengths with an inner square of side length (see Figure 4). In the second example, we have CNTs of different lengths and (see Figure 5). The third example involves CNTs of equal length with (see Figure 6). In the fourth example, we consider CNTs of random lengths between and with (see Figure 6). The centers and angles of the rigid line inclusions are chosen randomly so that all lines are nonoverlapping. For these fourth all examples, we use .
Estimating the boundary values of the function , requires two steps:

Computing the domain . For this step, the total number of discretization points is .

Compute the boundary values of the function . The total number of discretization points is .
The total CPU time required by the two steps to compute is shown in Table 1.
To compute the values of and , we discretize the domain using a matrix W
with approximately million points in . Thus is a matrix of points that discretize the domain . Afterwards, we compute the values of the temperature distribution and the heat flux at the points of the matrix using the formulas (16) and (18), respectively. For the temperature distribution , we use the MATLAB function contourf
to plot the contours of the function on the domain . To visualize the heat flux , we present a phase portrait with modulus contour lines of [21]. The colors represent the phase of where red means is in the direction of the positive real axis, green for the positive imaginary axis, cyan for the negative real axis, and violet for the negative imaginary [21]. We show also the contour lines of the amplitude of .
The total CPU time required to compute the values of and are shown in Table 1.
For the four examples, we plot in Figure 8 the computed values of the real constants (sorted from the smallest to the largest). We see from the figures that the sorted values of these constants almost lie on a straight line for large .
These examples demonstrate the computational effectiveness of the developed method of integral equations. We select a particular type of 2D structures, namely, a square with a large square hole. Such a hole models a cavity or a defect. The slits model perfectly conducting nanotubes randomly distributed in the bulk medium. The number of slits increasing from to clearly display the local fields in the considered composites. Different inclinations of slits (nanotubes) in Figure 4 show that the slit perpendicular to the external flux does not essentially disturb it contrary to the slit parallel to the external flux. The flux intensity can increase in 2 times due to the nanotube. Such an irregularity of the flux was noted in [19, 26] for analogous problems. With increasing the number of slits the perturbations of the local fields become smaller and a heterogeneous structure can be homogenized. The same examples demonstrate the elastic stress perturbations in the antiplane problem. For elasticity problem, the modulus of the flux has to be replaced by the stress intensity where
stands for the components of the stress tensor.
Number of rigid  Side length of  and  

line inclusions  the inner square  
4  0.5  6.69  46.76 
10  0.4  26.41  63.10 
123  0.3  209.03  63.10 
1005  0.1  1989.21  90.7 
Acknowledgments
We would like to thank Vladimir Mityushev for suggesting this research problem and for his helpful comments and discussions.
References
 [1] S.M. Lebedev, O.S. Gefle, E.T. Amitov, E.S. Nin, A.E. Bezrodny, M.R. Predtechenskiy, Conductive carbon nanotubereinforced polymer composites and their characterization, IEEE Transactions on Dielectrics and Electrical Insulation Volume: 23 , Issue: 3 , June 2016
 [2] Robert J. Young, Ian A. Kinloch, Lei Gong, Kostya S. Novoselov, The mechanics of graphene nanocomposites: A review. Composites Science and Technology 72 (2012) 14591476.
 [3] S.A.A. AlHatemi, A.H.M. Murid and M.M.S. Nasser. A boundary integral equation with the generalized Neumann kernel for a mixed boundary value problem in unbounded multiply connected regions. Bound. Value Probl., Article number: 54, 2013.
 [4] K.E. Atkinson. The Numerical Solution of Integral Equations of the Second Kind. Cambridge University Press, Cambridge, 1997.
 [5] P. Drygaś. New approach to mathematical model of elastic in twodimensional composites. In P. Drygaś and S. Rogosin, editors, Modern Problems in Applied Analysis, pages 87–100. Birkhäuser, 2018.
 [6] F.D. Gakhov. Boundary Value Problems. Pergamon Press, Oxford, 1966.
 [7] L. Greengard and Z. Gimbutas. FMMLIB2D: A MATLAB toolbox for fast multipole method in two dimensions. Version 1.2, http://www.cims.nyu.edu/cmcl/fmm2dlib/fmm2dlib.html. Accessed 1 Jan 2018.
 [8] J. Liesen, O. Séte and M.M.S. Nasser. A fast and accurate computation of the logarithmic capacity of compact sets. Comput. Methods Funct. Theory, 17:689–713, 2017.
 [9] V. Mityushev. Mixed problem for Laplace’s equation in an arbitrary circular multiply connected domain. In P. Drygaś and S. Rogosin, editors, Modern Problems in Applied Analysis, pages 135–152. Birkhäuser, 2018.
 [10] V.V. Mityushev and S.V. Rogosin. Constructive Methods to Linear and Nonlinear Boundary Value Problems of the Analytic Function. Chapman & Hall/CRC, Boca Raton, 2000.
 [11] N.I. Muskhelishvili. Some Basic Problems of the Mathematical Theory of Elasticity. SpringerScience+Business Media, Dordrecht, 1977.
 [12] M.M.S. Nasser. Numerical conformal mapping of multiply connected regions onto the second, third and fourth categories of Koebe’s canonical slit domains. J. Math. Anal. Appl., 382:47–56, 2011.
 [13] M.M.S. Nasser. Numerical computing of preimage domains for bounded multiply connected slit domains. J. Sci. Comput., 78:582–606, 2019.
 [14] M.M.S. Nasser. Fast solution of boundary integral equations with the generalized Neumann kernel. Electron. Trans. Numer. Anal., 44:189–229, 2015.
 [15] M.M.S. Nasser and C.C. Green. A fast numerical method for ideal fluid flow in domains with multiple stirrers. Nonlinearity, 31:815–837, 2018.
 [16] M.M.S. Nasser, A.H.M. Murid, M. Ismail and E.M.A. Alejaily. Boundary integral equations with the generalized Neumann kernel for Laplace’s equation in multiply connected regions. Appl. Math. Comput., 217:4710–4727, 2011.
 [17] N. Nishimura and Y.J. Liu. Thermal analysis of carbonnanotube composites using a rigidline inclusion model by the boundary integral equation method. Comput. Mech., 35:1–10, 2004.
 [18] E. Pesetskaya, R. Czapla and V. Mityushev. An analytical formula for the effective conductivity of 2D domains with cracks of high density. Appl. Math. Model., 53:214–222, 2018.
 [19] N. Rylko. Edge effects for heat flux in fibrous composites. Comput. Math. Appl., 70:2283–2291, 2015.
 [20] L.N. Trefethen and J.A.C. Weideman. The exponentially convergent trapezoidal rule. SIAM Review, 56:385–458, 2014.
 [21] E. Wegert. Visual complex functions. An introduction with phase portraits. Birkhäuser/Springer Basel AG, Basel, 2012.
 [22] R. Wegmann and M.M.S. Nasser. The RiemannHilbert problem and the generalized Neumann kernel on multiply connected regions. J. Comput. Appl. Math., 214:36–57, 2008.
 [23] G.C. Wen. Conformal Mapping and Boundary Value Problems. Amer. Math. Soc., Providence, 1992.
 [24] A.A.M. Yunus, A.H.M. Murid and M.M.S. Nasser. Numerical conformal mapping and its inverse of unbounded multiply connected regions onto logarithmic spiral slit regions and straight slit regions. Proc. R. Soc. A 470:20130514, 2013.
 [25] J. Zhang, M. Tanaka and T. Matsumoto. A simplified approach for heat conduction analysis of CNTbased nanocomposites. Comput. Meth. Appl. Mech. Eng. 192:5597–5609, 2004.
 [26] N. Rylko. Fractal local fields in random composites. Computers & Mathematics with Applications 69(3): 247254 (2015)
Comments
There are no comments yet.