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

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.

## Authors

• 4 publications
• 2 publications
12/28/2021

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

We consider the steady heat conduction problem within a thermal isotropi...
10/11/2020

### 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...
11/03/2021

01/15/2020

### 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 ...
03/22/2021

### 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...
01/04/2021

### The hot spots conjecture can be false: Some numerical examples

The hot spots conjecture is only known to be true for special geometries...
10/13/2021

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

##### 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

 ηk(t)=zk+0.5eiβkak(cost−irsint),t∈Jk,k=1,2,…,m, (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

 Φ(∞)=∞,limw→∞(Φ(w)−w)=0. (3)

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

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

 z=Φ(w)=w+12πi∫^ΓΦ(τ)−ττ−wdτ.

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]

 Φ−1(∞)=∞,limz→∞(Φ−1(z)−z)=0.

Therefore the values of for can be computed by

 w=Φ−1(z)=z+12πi∫^LΦ−1(ξ)−ξξ−zdξ. (4)

For the parameterization of the boundary , we adopt the same notations used in [8, 12, 13, 14, 15, 22]. Let

 ^J=m⨆k=1Jk=m⋃k=1{(t,k):t∈Jk},

where and means for . In this way, we define a parameterization of the boundary on by

 ^η(t,k)=ηk(t),t∈Jk,k=1,…,m.

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

 ^ζ(t)=Φ(^η(t)),t∈^J.

Thus, by making use of (4), we can determine the values of the inverse map for through

 w=Φ−1(z) =z+12πi∫^JΦ−1(^ζ(t))−^ζ(t)^ζ(t)−z^ζ′(t)dt =z+12πi∫^J^η(t)−^ζ(t)^ζ(t)−z^ζ′(t)dt (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. 696-697]. Hence for , the curve is parameterized by

 ηk(t)=Φ−1(ζk(t)),t∈Jk,

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

 U(z)=u(Φ−1(z)),z∈D, (6)

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

 Δu =0,in G, (7a) u =δk,on Γk,k=1,2,…,m, (7b) ∫Γk∂u∂nds =0,k=1,2,…,m, (7c) ∂u∂n =0,on Γm+1, (7d) u(τ) =Re[τ]forτ∈Γm+2. (7e)

As is harmonic, we can write

 u(w)=Re[f(w)] (8)

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

 f(η(t))=ϕ(t)+iψ(t),t∈J

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

 ϕ(t)=δk,t∈Jk,k=1,2,…,m. (9)

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

 ϕ(t)=Re[ηm+2(t)],t∈Jm+2. (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.,

 ψ(t)=δm+1,t∈Jm+1. (11)

where is an undetermined real constant.

Let us define the function on by

 ^A(t)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1,t∈J1,⋮1,t∈Jm,−i,t∈Jm+1,1,t∈Jm+2.

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

 (12)

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

 g(w)=f(w)−cw−α,w∈G∪Γ. (13)

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

 A(t)=e−iθ(t)(η(t)−α),t∈J,

where is the piecewise constant function given by

 θ(t)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩0,t∈J1,⋮0,t∈Jm,π/2,t∈Jm+1,0,t∈Jm+2.

Therefore the boundary conditions in (12) implies that solves the Riemann-Hilbert problem

 Re[A(t)g(η(t))]=γ(t)+h(t),t∈J,

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

 A(t)g(η(t))=γ(t)+h(t)+iμ(t),t∈J,

where uniquely solves the integral equation

 (I−N)μ=−Mγ (14)

and the function is computed by

 h=[Mμ−(I−N)γ]/2. (15)

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

 Nμ(s)=∫J1πIm(A(s)A(t)η′(t)η(t)−η(s))μ(t)dt,s∈J, Mμ(s)=∫J1πRe(A(s)A(t)η′(t)η(t)−η(s))μ(t)dt,s∈J.

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

 f(η(t))=(η(t)−α)g(η(t))+c.

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

 U(z)=Re[F(z))] (16)

where

 F(z)=f(Φ−1(z)),z∈D, (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

 q=−λ∇U=−λ(∂U∂x,∂U∂y)

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

 F′(z)=∂U∂x−i∂U∂y.

One the other hand (17) yields

 F′(z)=f′(Φ−1(z))Φ′(Φ−1(z)),z∈D,

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

 q(z)=−(∂U∂x,∂U∂y)∣∣∣z=−¯¯¯¯¯¯¯¯¯¯¯¯¯F′(z)=−¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯(f′(Φ−1(z))Φ′(Φ−1(z))),z∈D. (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 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.

## 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 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, 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 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)