Higher-Order Line Element Analysis of Potential Field with Slender Heterogeneiti

H.-S.Wang,H.Jiang,B.Yang

1 Introduction

Potential field problems described by Laplace(and generally,Poisson)equation are important in many areas such as thermal conduction,potential flow,and electrostatics[Telles and Wrobel(1984);Liu(2009)].Very often slender bodies(i.e.,heterogeneities)are involved,such as nanowires in nano-and micro-electro-mechanical systems,[Chen,Mukherjee,and Aluru(2008);Shen(2009)]and nanotube and glass fibers in composite materials[Nishimura and Liu(2004);Wang and Yao(2013);Chen and Mukherjee(2006);Fiamegkou,Athanasopoulos,and Kostopoulos(2014);Han and Fina(2011)].Many biological materials such as protein,[Simonson(2003);Grochowski and Trylska(2008)]collagens and fibrins,[Dong(2009);Jiang,Yang,and Liu(2013)]and solvated biomolecular system[Lu(2008)]fall in this class as well.The computational analysis of a system of slender bodies poses great challenge due to the high aspect ratio.The literature has seen scarce computational work on the potential field analysis with slender bodies despite its importance,likely due to difficulties arising from the high aspect ratio.

Finite difference(FD)method[Mitchell and Griffiths(1980)]and finite element(FE)method[Reddy and Gartling(2011)]are popular domain-based numerical techniques.By them,the domain is discretized into grids/elements,and a field over the domain is approximated piecewise with basis functions(such as polynomials)over local grids/elements.They are versatile,capable of dealing with virtually any physical problems.However,when they are applied to solve slender-body problems,a fine grid/mesh is required to reasonably resolve the characteristics of a slender body and its surrounding domain.Even if an adaptive mesh is used,the computational task can be prohibitively large.Boundary element(BE)method is another popular numerical technique.It is derived from a boundary integral equation and only involves numerical discretization on the boundary of a domain.The dimensionality reduction from a domain to a surface,typically by applying analytical integral kernels,is advantageous compared to the FD and the FE methods.However,the resulting stiffness matrix is full in contrast to the sparse matrix of the FD and the FE methods.When it is applied to solve slender-body problems,the computational task remains prohibitively large.In order to effectively and accurately solve slender-body problems,further mesh reduction is required.For potential field problems,[Chen and Mukherjee(2006)]introduced a line integral equation model summing the surface charge density around a slender heterogeneity into a line charge density meanwhile neglecting higher-order terms,such as line dipole density due to uneven charge distribution around the perimeter.On the other hand,such higher-order effect may become significant when the material property mismatch between the slender body and the surrounding medium is significant or when slender bodies closely interact with each other or with other objects.A line integral equation formulation has been applied to the hydrodynamics of a slender body with only side net force density[Pozrikidis(2011);Tornberg and Shelley(2004)]and with both net force and moment densities[Jiang,Wu,Zhao,and Yang(2014);Jiang and Yang(2013)].In the present work,we introduce a higher-order line integral equation technique for potential field problems with slender heterogeneities,examine its efficiency and accuracy,and examine the characteristics of the potential field around a slender body.

The rest of the paper is organized as follows.In Sec.2,we derive the line integral equations of potential field from the general boundary integral equations along a slender body by expanding the integral kernels in Taylor series about the slender body central line and retaining the first two integral terms plus a finite term.The finite terms of a cylinder are derived,which carry the local information of a crosssectional shape and represent the core structure behind the singularities.In Sec.3,the numerical treatments to the line integral equation formulation are discussed.In Sec.4,four examples are presented to demonstrate the capability and validity of the present method,in the context of steady-state heat conduction.They are:(a)temperature distribution along two parallel rods,both sources;(b)temperature distribution along two parallel rods,one source and the other sink;(c)temperature distribution along a spiral rod;(d)temperature distribution along a three-dimensional(3D)architecture of rods.In the cases of two parallel rods,the problem is examined with their distance as the parameter to show when the dipolar term is important and when the line model is accurate.In Sec.5,conclusions are drawn.

Figure 1:Schematic of two parallel rods with separation distance.

2 Line integral equations of potential field around a slender heterogeneity

2.1 Line integral equations along a slender cavity

Consider a problem of a scalar potential field around a slender cavityCin a domainΩ,as shown in Fig.1.The fixed Cartesian coordinate system(x1,x2,x3;O)is established for reference.Meanwhile,a Cartesian body frame(X1,X2,X3;xc)is established with the second axis directed along the longitudinal direction and the normal axes residing within a cross section,and origin xcbeing the centroid of the same cross section,ofC.A lower-case subscript index will be used to indicate a component in the global frame,while an upper-case subscript index to indicate a component in the body frame.Whenever a numerical number or a Greek letter is used to indicate a component,it refers to the body frame.(x,y,z)and(X,Y,Z)may also be used to indicate the components,in the global and the body frames,respectively,and they should not be confused with their bold-faced or indexed counterpart of a vector or a component of a vector.

The scalar potential fieldTsatisfies the Laplace equation,∇2T=0,within Ω.By applying the reciprocal theorem,the following boundary integral equation can be derived as

where Γ is the boundary of Ω (including the surface of cavityCand a possible external surface if Ω is finite),J≡−∇T·n n n,and n is the outward unit vector at a point on Γ.Constantntegral kernelsT∗andJ∗are the fundamental solutions satisfying∇2T∗+δ =0 andJ∗≡−∇T∗·n n n,where δ is the Dirac delta function representing a point source,in an infinite space.They are given in terms of source point y and field point x by

In order to derive the corresponding line integral equation,let us approximate the potential field along the perimeter of a cross section of the slender cavity as

where Γscis the slender cavity boundary,and Γexis the external boundary.Examining it in the vicinity of yc,one may realize that the term on the left-hand side of Eq.(4a)is equal toT(yc)plus residual on the order of ρ2/L2,whereLis a characteristic length scale,for instance,the radius of curvature of the slender cavity,the distance between the slender body and an external boundary,and the distance between slender bodies if there are multiple of them–whichever is the smallest.Thus,the above equation may be rewritten as

We then expandT∗in Taylor series(up to the linear term)in the vicinity of xcand use it to approximate Eq.(5)as

where the integrals are taken in the sense of Cauchy or Hadamard principal value depending on the odd or even order of a singularity in the kernels,andTftis a finite term resulting from the limiting process of xc→yc.By realizing that the integral kernels are now independent of local circumferential direction(s)of across section,we rewrite Eq.(6)as

In order to complete the formulation,higher-order integral equations are required.By taking derivative of all terms in Eq.(1)with respect to source point y,and following the same procedure of derivation as above,we can obtain the line integral equation of potential gradient.For the sake of brevity,the derivation is omitted.This equation is given by

The above equations(7)and(8)are termed as the dual line-integral equations for perimeter-average potentialTaand nominal potential gradientat a slender cavity.They are accurate on the order of ρ/L,with higher-order terms neglected.The loading terms of net source and source dipole are taken into account.In the following,we shall fill the slender cavity with a slender body so that the above line integral equations can be further reduced.

2.2 Line integral equation formulation of a slender heterogeneity

Let us fill in the above slender cavityCwith a heterogeneityHof the same shape.The potential field insideHis modeled as one-dimensional,and is characterized by using average potentialTHand nominal gradientsandalong local axesXandYover a cross section,along central linel.The conjugate loading terms are correspondingly heat source and dipoles,QH,and,defined over a cross section.By assuming the continuity conditions of both potential and source fields across the interface between the heterogeneity and surrounding medium,and substituting the heterogeneity fields in Eqs.(7)and(8),we have

where α and β range from 1 to 2,representingXandYin the body frame.

2.3 Finite terms

The line integrals in Eqs.(6)–(10)are singular when the source and field points coincide.They are taken in the sense of Cauchy or Hadamard principal value depending on their odd or even order.It means that these line integrals do not contain any physical information inside the core of the singularity,such as the cross-sectional shape and the interfacial continuity condition.Thus,it must be added back to the line-integral equations to represent the physics correctly.In light of this understanding,these finite terms are defined as

In the present work,we specify the slender cavityC(and hence,the matching heterogeneityH)to be of a circular cross section.Also,it is assumed that the central linelis only slightly curved so that it can be reasonably approximated as straight locally at its any middle location.Under these conditions,the finite terms,Tft,andcan be derived.We have analytically integrated over an in fi nite slender cavity and obtained the following results:

whereRis the radius of the circular cross section.AssigningQandDαto slender heterogeneity,it was implied that the fields are continuous across the interface.

3 Numerical line element method

In this section,we introduce a numerical scheme to solve the problem of potential field with a slender heterogeneity as formulated above.A slender heterogeneity,modeled as a line object,is discretized by using piecewise straight elements.A single node is placed at the middle point of each element.Each node is associated with 3 degrees of freedom,average potentialTH,and nominal potential gradient within a cross-sectional plane,,for β =1,2.The discretized versions of Eqs.(9)and(10)are given by

where superscriptn/mdenotes thenth/mth node,andNis the total number of nodes.

The nodal finite terms are:in whichRis taken to be the slender-body radius at thenth node.Note that only a remote external boundary is considered here whose corresponding terms,if any,are treated as prescribed.

4 Numerical examples

In this section,the proposed LE method is demonstrated by using examples in heat conduction.At first,the case of two parallel rods with identical,higher temperature in an infinite medium is used as a benchmark problem to demonstrate the validity of the method.By the process of retrieving the local information,more physical interpolations of the finite terms are given.It demonstrates that the proposed LE method is valid at the slender heterogeneity limit compared to the BE method.In order to further show the capability of the proposed LE method,three more applications of(b)two anti-symmetric parallel rods with separation distance;(c)a helical rod;and(d)3-dimentional cross architecture of rods in a temperature field,are discussed.Compared with the BE solution,the LE solution is shown to be valid and far more efficient,which can be expected as a powerful tool for the study of potential field problems.

Figure 2:Schematic of two parallel rods with separation distance.

4.1 Two parallel rods,both sources

Two parallel rods as shown in Fig.2 are used in the first case study.In this case,the remote temperature is set to be zero,i.e.T∞=0.The temperature of both rods is specified as one unit,representing a symmetric condition.The rod radiusais used as the characteristic length scale and the length of each rod is 100 times of that.Thus,the high aspect ratio satisfies the slender body assumption.For the purpose of validity demonstration of the proposed LE method,the BE method is used to provide a baseline solution in this case.

Figure 3:Variation of heat flux along a slender rod,predicted with various mesh densities.

For the case of rod separation distanced=3a,the LE results of heat flux along the rod obtained with 10,20,25 elements are shown in Fig.3.The corresponding BE result is also shown,which is due to a mesh of 100 equal divisions in length,24 equal divisions in circumference,and 5 equal divisions in radius,over the rod surface.This mesh was checked to yield converged BE solutions.The results presented in Fig.3 show that the solutions of LE method are stable and agree with the results obtained by using the BE method.

Secondly,we examine the separation distance dependence of the LE solution.Figures 4(a)and(b)show the heat flux and flux gradient,respectively,when the separation distance between the two rods are 3a,5a,and 9a.The results show excellent agreement between the LE and BE solutions,especially when the separation distance is relatively large.As shown in Fig.4(a),the heat flux on the two rods are the same,exhibiting the symmetry.When the separation distance is decreasing,the value of the heat flux decreases.As shown in Fig.4(b),the heat flux nonuniformality around the rod perimeter becomes larger when they get close to each other.

Figure 4:(a)Variation of heat flux(Q)and(b)Variation of heat flux gradient(Q,1)around two slender rods,at separation distance d=3a,5a and 9a.

Upon the above mesh and separation distance dependence analysis,we present some detailed results for further analysis and discussion.Recalling in the process of line integral derivation,the net source and source dipole are the integrations along the perimeter of the cross-section,we are able to retrieve the local information by curve fitting using the least square(LSQ)method.Let us assume that the local heat fl ux can be expressed as

whereJnis the total heat flux at a specific position on the center line of the slender body;aandbare constants,which can be interpolated asandrespectively;∅is the azimuth angle ranging from 0 to 2π.If the position ofyis specified as the middle point on the center line of the rod,the curve fitting results and corresponding residual errors for different separation distances between the two rods are shown in Fig.5(a)and(b).The parameters for curve fitting are listed in Table 1.

Meanwhile,the second order non-linear curve fitting to the full solution is tested,which is expressed by

wherecanddare additional constants,which can be interpolated asandH,respectively.The curve fitting results and corresponding residual errors for various separation distances between the two rods are shown in Figs.6(a)and(b),respectively.The parameters obtained from the curve fitting are listed in Table 2.

Figure 5:The first order nonlinear curve fit result.

Table 1:First order nonlinear curve fit parameters result.

As shown in Figs.5 and 6,both of the first-order and second-order curve fittings can fit the local heat flux qualitatively well.It indicates that the local information dropped in the integration reduction can be retrieved if the geometry of the crosssection is known.When the two rods are far enough away from each other,the high-order term can be neglected,and the first order term can be used attaining enough numerical accuracy.When the separation distance decreases,the fitting results deviate from the full solution,as may be expected.This is because the high order terms become more significant as the two rods get close to each other.By comparing the results in Figs.5(a)and Fig.6(a),it can be seen that when the separation distance is comparable to the diameter of the rod,the high order curve fitting gives more accurate fitting results,suggesting that the higher-order term(＞dipole)should be considered to achieve higher quantitative accuracy.

Figure 6:The second order nonlinear curve fit result.

Table 2:Second order nonlinear curve fit parameters result.

4.2 Two parallel rods,one source and the other sink

In the second example,the thermal distribution of two anti-symmetric parallel rods with various separation distances is analyzed.The geometry of both rods is the same as in the above case.The heat flux condition is prescribed withQ=1 unit on one rod andQ=−1 unit on the other.They represent a source and a sink,respectively.No dipole is prescribed on either rod.The LE results of temperature and temperature gradient along the two rods obtained with 20 elements for selected separation distance are shown in Fig.7.

Figure 7:(a)Variation of temperature(T)and(b)temperature gradient(T,1)distribution around the two anti-symmetry rods when d=3a,5a,7a,9a,11a.

As shown in Fig.7(a),the temperature on the two rods have the same value but different signs,exhibiting the skew symmetry.When the separation distance is small,the absolute value of the temperature approaches to zero.When they get close to each other,the source and sink of equal intensity basically annihilate.As shown in Fig.7(b),the temperature nonuniformality around the rod perimeter is trivial when the rods are far away from each other but becomes significant when they get close to each other.It might be worth noting that if the rods are thermoelastic,they would bend due to the varying temperature around the surface and inside the rods.

4.3 A helical rod

In the third case study,a helical rod as shown in Fig.8 is considered.The helical rod is of radius R,and of a center-line shape defined by

where ϕ is the azimuth angle,(x,y,z)are the coordinates in the global coordinate system,andc1,c2andc3are the constants.Herec1=c2=R,c3=50Rand ϕ =[0,6π]are used.A uniform heat fluxQequal to one unit is applied through the helix.The spiral slender body is discretized into 90 segments and each segment is treated as a straight rod.The line represents the helical rod,and the black solid dots show the nodes.The LE solution is superimposed on the rod curve,as shown in Fig.8.Each arrow starts from the location of a node,and the end of the arrow is at(x+T,1x,y+T,1y,z+T,1z).It represents the calculated temperature “gradient”at the specific position and illustrates the temperature nonuniformity.It can be seen that the temperature gradient vectors(i.e.,arrows)along the solid surface imply that the heat is released more into the center of the spiral.

Figure 8:A slightly curved slender spiral.The dots show the nodes used in the simulation and the arrows represent the predicted temperature gradient vectors along the spiral when heat flux is considered as a constant.

Figure 9:A 3d architecture of rods.The dots show the nodes used in the simulation and the arrows represent the predicted temperature gradient vectors along the spiral when heat flux is considered as a constant.

4.4 A 3D architecture of rods

A 3D cross architecture of rods in a temperature field is discussed in the fourth case.The rods in a layer are overlapped with right angles to the rods in the neighboring layers,which is defined by

where ϕx, ϕyand ϕzare the azimuth angles alongxaxis,yaxis andzaxis,separately.(x,y,z)are the coordinates in the global coordinate system,andc1,c2andc3are the constants.Herec1=c2=c3.ϕy=ϕx+π/2,ϕz=π/2.The geometry of the 3dcross architecture of rods is presented in Fig.14.The 3dcross architecture is stacked by 7 layers of lattices separated by the distance of 3.Each layer consists of 7 parallel rods and the neighboring rods are separated by the distance of 3.Each cross rod is discretized into 20 segments.The line with solid dots represents the 3D cross architecture of rods,and the black solid dots show the nodes as shown in Fig.9.

When heat flux is assumed as a constant(Q=1),the LE result of temperature gradient(T,1)distributed on the 3D cross architecture of rods in global frame is shown in Fig.9.As shown in this figure,the temperature gradient vectors(i.e.,arrows)of the cross architecture rods imply heat releasing into the center of 3D cross architecture.The amplitude of the temperature gradient is decreasing from the outer layer to the center layer but for each rod,the amplitude of the temperature gradient is increasing from the ends to its middle point.

5 Conclusion

We have presented a numerical analysis of the potential field problem by applying a second-order LE method.The LE method is generally a mesh reduction method further from the BE method,specially devised for slender heterogeneities.The derivation is rigorous by reducing from the boundary integral equation.It can be applied to obtain the potential difference and flux difference(dipole)on opposing sides of a slender heterogeneity as well as the average potential and total flux around it.Such higher-order terms can be important when to study close interaction of slender heterogeneities with walls and with themselves.In particular,we have examined in detail the benchmark case of two parallel rods by comparing the LE and the BE solutions.The parameterized LE solution is reconstructed by using sinusoidal functions.The reconstructed fields agree with the BE solution well when the rods distance is large(compared to the rod radius).They become different when the rods get very close to each as expected.This comparative study demonstrates the validity of the present LE method.Furthermore,three more examples are presented to demonstrate the capability of the present LE method,including a long spiral heterogeneity and a 3D cross architecture of rods emitting heat.The present LE method may find itself for various applications in MEMS,fiber composites,etc.involving slender heterogeneities.

Acknowledgement:This work is sponsored by the State Scholarship Fund from the China Scholarship Council(CSC),whose support(Grant No.201203070476)is gratefully acknowledged.

Chen,H.;Mukherjee,S.(2006):Charge distribution on thin conducting nanotubes–reduced 3-D model.International Journal for Numerical Methods in Engineering,vol.68,no.5,pp.503–524.

Chen,H.;Mukherjee,S.;Aluru,N.(2008):Charge distribution on thin semiconducting silicon nanowires.Computer Methods in Applied Mechanics and Engineering,vol.197,no.41,pp.3366–3377.

Dong,B.(2009):Electrospinning of collagen nanofiber scaffolds from benign solvents.Macromolecular Rapid Communications,vol.30,no.7,pp.539–542.

Fiamegkou,E.;Athanasopoulos,N.;Kostopoulos,V.(2014):Prediction of the effective thermal conductivity of carbon nanotube-reinforced polymer systems.Polymer Composites,vol.35,no.10,pp.1997–2009.

Grochowski,P.;Trylska,J.(2008):Continuum molecular electrostatics,salt effects,and counterion binding–a review of the Poisson-Boltzmann theory and its modifications.Biopolymers,vol.89,no.2,pp.93–113.

Han,Z.;Fina,A.(2011):Thermal conductivity of carbon nanotubes and their polymer nanocomposites:a review.Progress in Polymer Science,vol.36,no.7,pp.914–944.

Jiang,H.;Wu,Y.T.;Zhao,Y.-P.;Yang,B.(2014):Force-moment line element method for Stokes flow around a slender body.Engineering Analysis with Boundary Elements,vol.44,pp.120–129.

Jiang,H.;Yang,B.(2013):Force-moment line element method for flexible slender bodies in Stokes flow.Physical Review E,vol.88,no.3.

Jiang,H.;Yang,B.;Liu,S.(2013):Effects of transverse shear on strain stiffening of biological fiber networks.CMC:Computers,Materials&Continua,vol.38,no.2,pp.61–77.

Liu,Y.(2009):Fast multipole boundary element method:theory and applications in engineering,Cambridge university press.

Lu,B.(2008):Recent progress in numerical methods for the Poisson-Boltzmann equation in biophysical applications.Commun Comput Phys,vol.3,no.5,pp.973–1009.

Mitchell,A.R.;Griffiths,D.F.(1980):The finite difference method in partial differential equations:John Wiley.

Nishimura,N.;Liu,Y.(2004):Thermal analysis of carbon-nanotube composites using a rigid-line inclusion model by the boundary integral equation method.Computational Mechanics,vol.35,no.1,pp.1–10.

Pozrikidis,C.(2011):Shear flow past slender elastic rods attached to a plane.International Journal of Solids and Structures,vol.48,no.1,pp.137–143.

Reddy,J.N.;Gartling,D.K.(2011):The finite element method in heat transfer and fluid dynamics:CRC press.

Shen,G.(2009):Devices and chemical sensing applications of metal oxide nanowires.Journal of Materials Chemistry,vol.19,no.7,pp.828–839.

Simonson,T.(2003):Electrostatics and dynamics of proteins.Reports on Progress in Physics,vol.66,no.5,pp.737.

Telles,J.;Wrobel,L.(1984):Boundary element techniques.Theory and applications in engineering,Springer-Verlag,Berlin.

Tornberg,A.K.;Shelley,M.J.(2004):Simulating the dynamics and interactions of flexible fibers in Stokes flows.Journal of Computational Physics,vol.196,no.1,pp.8–40.

Wang,H.;Yao,Z.(2013):Large-scale thermal analysis of fiber composites using a line-inclusion model by the fast boundary element method.Engineering Analysis with Boundary Elements,vol.37,no.2,pp.319–326.