Application of integral equations to simulate local fields in carbon nanotube reinforced composites

by   Mohamed M S Nasser, et al.

We consider the steady heat conduction problem within a thermal isotropic and homogeneous square ring composite reinforced by non-overlapping and randomly distributed carbon nanotubes (CNTs). We treat the CNTs as rigid line inclusions and assume their temperature distribution to be fixed to an undetermined constant value along each line. We suppose also that the temperature distribution is known on the outer boundary and that there is no heat flux through the inner square. The equations for the temperature distribution are governed by the two-dimensional Laplace equation with mixed Dirichlet- Neumann boundary conditions. This boundary value problem is solved using a boundary integral equation method. We demonstrate the performance of our approach through four numerical examples with small and large numbers of CNTs and different side length of the inner square.



There are no comments yet.



Simulating local fields in carbon nanotube reinforced composites for infinite strip with voids

We consider the steady heat conduction problem within a thermal isotropi...

Total heat flux convergence in the calculation of 2d and 3d heat losses through building elements

Heat losses through the building envelope is one of the key factors in t...

Estimación y Análisis de Sensibilidad para el Coeficiente de Difusividad en un Problema de Conducción de Calor

The aim of this article is to discuss the estimation of the diffusivity ...

On the discretization of Laplace's equation with Neumann boundary conditions on polygonal domains

In the present paper we describe a class of algorithms for the solution ...

Non-iterative domain decomposition for the Helmholtz equation using the method of difference potentials

We use the Method of Difference Potentials (MDP) to solve a non-overlapp...

The hot spots conjecture can be false: Some numerical examples

The hot spots conjecture is only known to be true for special geometries...

Solving multiscale steady radiative transfer equation using neural networks with uniform stability

This paper concerns solving the steady radiative transfer equation with ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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


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


For more details, we refer the reader to [15, 13]

Figure 2: The domain exterior to rectilinear slits (right) and the preimage domain exterior to ellipses (left).

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


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


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. 696-697]. 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

Figure 3: The given domain (right) and the preimage domain (left) for .

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


such that is the unique solution of the following boundary value problem on the domain in -plane,


As is harmonic, we can write


where is a single-valued analytic function on the domain . Assuming the boundary values of the function takes the form

on the boundary of , it follows from (7b) that


Similarly, equation (7e) implies that the values of the function are known for with


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.,


where is an undetermined real constant.

Let us define the function on by

Following (9), (10) and (11), the function satisfies the boundary conditions


which are known as the Riemann-Hilbert 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


The function is clearly analytic in . We consider also a complex-valued function defined on by

where is the piecewise constant function given by

Therefore the boundary conditions in (12) implies that solves the Riemann-Hilbert 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


and the function is computed by


Here, is the identity operator and the integral operators and are respectively defined by

For more details, see [12, 14].

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 matrix-vector 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=1e-14, 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




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 Cauchy-Riemann 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


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 right-hand 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 non-overlapping. For these fourth all examples, we use .

Estimating the boundary values of the function , requires two steps:

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

  2. 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 anti-plane 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
Table 1: The total CPU time in seconds required to compute the boundary values of the complex distribution temperature in -plane; and the values of the temperature distribution and the heat flux .

Figure 4: A contour plot of the temperature distribution (left) and a phase portrait of the heat flux (right) for CNTs and .

Figure 5: A contour plot of the temperature distribution (left) and a phase portrait of the heat flux (right) for CNTs and .

Figure 6: A contour plot of the temperature distribution (left) and a phase portrait of the heat flux (right) for CNTs with .

Figure 7: A contour plot of the temperature distribution (left) and a phase portrait of the heat flux (right) for CNTs with .

Figure 8: The (sorted) computed values of the real constants for (top, left), (top, right), (bottom, left), and for (bottom, right).


We would like to thank Vladimir Mityushev for suggesting this research problem and for his helpful comments and discussions.


  • [1] S.M. Lebedev, O.S. Gefle, E.T. Amitov, E.S. Nin, A.E. Bezrodny, M.R. Predtechenskiy, Conductive carbon nanotube-reinforced 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) 1459-1476.
  • [3] S.A.A. Al-Hatemi, 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 two-dimensional 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, 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 Non-linear 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. Springer-Science+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 carbon-nanotube composites using a rigid-line 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 Riemann-Hilbert 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 CNT-based nano-composites. Comput. Meth. Appl. Mech. Eng. 192:5597–5609, 2004.
  • [26] N. Rylko. Fractal local fields in random composites. Computers & Mathematics with Applications 69(3): 247-254 (2015)