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

We consider the steady heat conduction problem within a thermal isotropic and homogeneous infinite strip composite reinforced by uniformly and randomly distributed non-overlapping carbon nanotubes (CNTs) and containing voids. We treat the CNTs as thin perfectly conducting elliptic inclusions and assume the voids to be of circular shape and act as barriers to heat flow. We also impose isothermal conditions on the external boundaries by assuming the lower infinite wall to be a heater under a given temperature, and the upper wall to be a cooler that can be held at a lower fixed temperature. The equations for the temperature distribution are governed by the two-dimensional Laplace equation with mixed Dirichlet-Neumann boundary conditions. The resulting boundary value problem is solved using the boundary integral equation with the generalized Neumann kernel. We illustrate the performance of the proposed method through several numerical examples including the case of the presence a large number of CNTs and voids.

There are no comments yet.

## Authors

• 2 publications
• 2 publications
• 1 publication
• 1 publication
10/13/2019

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

We consider the steady heat conduction problem within a thermal isotropi...
11/03/2021

11/17/2020

### Laplace Green's Functions for Infinite Ground Planes with Local Roughness

The Green's functions for the Laplace equation respectively satisfying t...
06/04/2021

### Unconditionally energy stable and first-order accurate numerical schemes for the heat equation with uncertain temperature-dependent conductivity

In this paper, we present first-order accurate numerical methods for sol...
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...
04/11/2022

### Linear Moment Models to Approximate Knudsen Layers

We propose a well-posed Maxwell-type boundary condition for the linear m...
03/31/2020

### Controlling Rayleigh-Bénard convection via Reinforcement Learning

Thermal convection is ubiquitous in nature as well as in many industrial...
##### This week in AI

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

## 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 macro-composites [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 (two-dimensional) thermal isotropic and homogeneous nanocomposite, which takes the form of an infinite strip, when it is reinforced by non-overlapping 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 non-conducting 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 non-overlapping, 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 steady-state 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 steady-state conditions, Fourier’s law of heat conduction and the above specified heat boundaries conditions yield the temperature distribution governed by the following mixed Dirichlet-Neumann boundary value problem:

 ΔT =0in Ω, (1a) T =0on C′0, (1b) T =1on C′′0, (1c) T =δkon Ck,k=1,2,…,ℓ, (1d) ∫Ck∂T∂nds =0k=1,2,…,ℓ, (1e) ∂T∂n =0on Ck,k=ℓ+1,ℓ+2,…,m, (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

 z=Φ(ζ)=1πlog1+ζ1−ζ+i2

conformally maps the unit disk onto the infinite strip . Thus, the inverse mapping

 ζ=Φ−1(z)=tanh(πz2−πi4)

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

 T(z)=U(Φ−1(z))

in which the function is the solution of the following boundary value problem in the -plane:

 ΔU =0in G, (2a) U =0on Γ′0, (2b) U =1on Γ′′0, (2c) U =δkon Γk,k=1,2,…,ℓ, (2d) ∫Γk∂U∂nds =0k=1,2,…,ℓ, (2e) ∂U∂n =0on Γk,k=ℓ+1,ℓ+2,…,m, (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

 U(ζ)=u0(ζ)+u(ζ)

where is a harmonic function in , and

 u0(ζ)=1πImlog1−ζ1+ζ+12.

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

 Δu(ζ) =0if ζ∈G, (3a) u(ζ) =0if ζ∈Γ0, (3b) u(ζ) =δk−1πImlog1−ζ1+ζ−12if ζ∈Γk,k=1,2,…,ℓ, (3c) ∫Γk∂u∂nds =0k=1,2,…,ℓ, (3d) ∂u∂n∣∣∣ζ =−∂u0∂n∣∣∣ζif ζ∈Γk,k=ℓ+1,ℓ+2,…,m, (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 ,

 n(ζ)=eiν(ζ)=−iη′(t)|η′(t)|. (4)

Thus

 ∂u0∂n=∇u0⋅n=cosν∂u0∂x+sinν∂u0∂y=Re[eiν(∂u0∂x−i∂u0∂y)]. (5)

The harmonic function is the real part of a single-valued analytic function , i.e., , where

 f0(ζ)=1πilog1−ζ1+ζ+12, (6)

and the branch of the logarithm function is chosen such that . Then by the Cauchy-Riemann equations, we have

 f′0(ζ)=∂u0(ζ)∂x−i∂u0(ζ)∂y,

which, in view of (4) and (5), implies that

 |η′(t)|∂u0∂n∣∣∣η(t)=Re[−iη′(t)f′0(η(t))],η(t)∈Γk,k=ℓ+1,…,m. (7)

Since

 f′0(ζ)=iπ(11−ζ+11+ζ),

it follows that for and ,

 |η′(t)|∂u0∂n∣∣∣ζ=η(t)=1πRe[η′(t)1−η(t)+η′(t)1+η(t)]. (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

 Re[f(η(t))]=0,η(t)∈Γ0 (9)

and

 Re[f(η(t))]=δk−1πImlog1−η(t)1+η(t)−12if η(t)∈Γk,k=1,…,ℓ. (10)

For the remaining boundary components for , we use the condition (3e) to determine the boundary condition on . By the Cauchy-Riemann equations, we can show using the same arguments as in (7) that

 |η′(t)|∂u∂n∣∣∣η(t)=Re[−iη′(t)f′(η(t))] (11)

Thus, for for , it follows from (3e), (8), and (11) that

 Re[−iη′(t)f′(η(t))]=−1πRe[η′(t)1−η(t)+η′(t)1+η(t)].

Integrating with respect to the parameter for , , we obtain

 Re[−if(η(t))]=1πlog∣∣∣1−η(t)1+η(t)∣∣∣+δk, (12)

where the integration constants are undetermined. The constants , in (10) and (12) are determined so that is a single-valued 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

 f(ζ)=(ζ−α)g(ζ)+c. (13)

Define also

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

where is the piecewise constant function given by

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

Thus

 e−iθ(t)f(η(t))=A(t)g(η(t))+e−iθ(t)c,

which implies that

 Re[A(t)g(η(t))]=Re[e−iθ(t)f(η(t))]−ccosθ(t).

On the basis of the conditions (9), (10), and (12), the function satisfies the Riemann-Hilbert problem

 Re[A(t)g(η(t))]=γ(t)+h(t), (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

 g(η(t))=γ(t)+h(t)+iμ(t)A(t),t∈J. (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

 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.

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

 (I−N)μ=−Mγ. (19)

Additionally, the piecewise constant function is given by

 h=[Mμ−(I−N)γ]/2. (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 matrix-vector 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

 U(ζ)=Re[f(ζ)+f0(ζ)].

We deduce the values of the temperature distribution for any by

 T(z)=Re[f(Φ−1(z))+f0(Φ−1(z))].

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

 F(z)=f(Φ−1(z))+f0(Φ−1(z)),z∈Ω.

According to the Cauchy-Riemann equations, it follows that the derivative of the complex potential on is given by

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

One the other hand,

 F′(z)=f′(Φ−1(z))Φ′(Φ−1(z))+f′0(Φ−1(z))Φ′(Φ−1(z)),z∈Ω, (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

 q(z)=−(∂T∂x,∂T∂y)∣∣∣z≡−¯¯¯¯¯¯¯¯¯¯¯¯¯F′(z). (22)

Hence

 ∂T∂y=−ImF′(z). (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 function

using trigonometric interpolating polynomials then differentiating. The values of

, in the right-hand 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

 λy=−12∫1−1∂T∂y(x,0)dx=12Im[∫1−1F′(x)dx]. (24)

Since

 ξ0(t)=Φ(η0(t))=Φ(eit)=1πlog1+eit1−eit+i2,0≤t≤2π,

where for , is on the line and for , is on the real line. Thus, for , we have

 ξ0(t)=1πlog∣∣∣cott2∣∣∣.

Since and where

 π

Consequently, (24) can be written as

 λy=12Im[∫t2t1F′(ξ0(t))ξ′0(t)dt]. (26)

In combining (21) with the fact that , we can see that

 F′(ξ0(t))=f′(Φ−1(ξ0(t)))Φ′(Φ−1(ξ0(t)))+f′0(Φ−1(ξ0(t)))Φ′(Φ−1(ξ0(t)))=f′(eit)Φ′(eit)+f′0(eit)Φ′(eit).

Hence,

 λy=12Im[∫t2t1[ieit(f′(eit)+f′0(eit))]dt], (27)

which implies that

 (28)

The second term in the right-hand side of (27) does not depend on the CNTs or the voids. In view of (6) and (25), we have

 12Im[f0(eit2)−f0(eit1)]=1,

and hence

 λy=1+12Im[f(eit2)−f(eit1)]. (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

 λy=1+12[μ(t2)−μ(t1)]. (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 non-overlapping circular holes and without any CNT (i.e., ). We also assume that all circular holes have the same radius with the parametrization

 ηj(t)=zj+re−it,0≤t≤2π,j=1,2,…,m,

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

 c(m,r)=mr2π2. (31)

The Clausius-Mossotti 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]

 λe=1−c1+c+O(c3). (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

 λy=0.8533491.

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

 λy=0.1519156.

When is close to , the circular holes become adjacent to each other. We compute the values of for several values of , . The obtained results are depicted in Figure 6 (left) where, by (31), the concentration of these holes is . Note that for .

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

 ηj(t)=zj+acost−ibsint,0≤t≤2π,j=1,2,…,m, (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

 c(m,a,b)=mabπ2. (34)

If , instead of (34) the plane slits density is considered in the theory of composites and porous media

 ϕ=mb2|Ω|=mb22, (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])

 λe=1+c1−c+O(c3). (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

 λy=1.0272480.

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

 λy=1.2804116.

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 .

The values of are also computed for several values of for with . Since is small, the obtained values of are plotted versus the values of , given by (35), where for . The obtained results are presented in Figure 8 (right).

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

When , the ellipses reduce to circles. We compute the values of for several values of , . The obtained results are presented in Figure 10 where the concentration of these ellipses is . For and , we have . Figure 10 presents also the values of the estimated effective conductivity given by (36).

### 6.3 The domain Ω with 2000 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 three-phases 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

 λe=1+c1−c21−c1+c2+O(c3), (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