- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Author Guidelines ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Table of Contents

Advances in Mechanical Engineering

Volume 2013 (2013), Article ID 942073, 10 pages

http://dx.doi.org/10.1155/2013/942073

## Mixed Volume and Boundary Integral Equation Method with Elliptical Inclusions and a Circular/Elliptical Void

Department of Mechanical and Design Engineering, Hongik University, 2639 Sejong-Ro, Jochiwon-Eup, Sejong City 339-701, Republic of Korea

Received 31 May 2013; Accepted 7 September 2013

Academic Editor: Indra Vir Singh

Copyright © 2013 Jung-Ki Lee et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A mixed volume and boundary integral equation method (mixed VIEM-BIEM) is used to calculate the plane elastostatic field in an unbounded isotropic elastic medium containing multiple isotropic/orthotropic elliptical inclusions of arbitrary orientation and a circular/elliptical void subjected to remote loading. In order to investigate the influence of a circular/elliptical void on the interfacial stress field, a detailed analysis of the stress field at the interface between the matrix and the central isotropic/orthotropic inclusion is carried out for the square packing of eight inclusions and one void, taking into account different values for the orientation angles and concentration of the inclusions. The mixed method is shown to be very accurate and effective for investigating the local stresses in composites containing isotropic/anisotropic fibers and a circular/elliptical void.

#### 1. Introduction

A micrograph of the cross section of a phosphate glass fiber/polymer composite is shown in Figure 1 [1]. Figure 1 indicates that the fibers are close to ellipses, and the major axis of the elliptical fibers is not aligned in any one direction. If voids are caused by manufacturing and/or service induced defects, they may have a significant effect on the local stresses in the composite.

A number of analytical techniques are available for solving stress analysis of inclusion problems when the geometry of the inclusions is simple (i.e., cylindrical, spherical, or ellipsoidal) and when they are well separated [2–4]. However, these approaches cannot be applied to more general problems where the inclusions are of arbitrary shape and their concentration is high. Thus, the stress analysis of heterogeneous solids often requires the use of numerical techniques based on the finite element method (FEM) or boundary element method (BIEM). Unfortunately, both methods encounter limitations in dealing with problems involving infinite media or multiple inclusions. However, it has been demonstrated that a recently developed numerical method based on a volume integral formulation can overcome such difficulties in solving a large class of inclusion problems [5–10]. One advantage of the volume integral equation method (VIEM) over the boundary integral equation method (BIEM) is that it does not require the use of Green’s functions for both the matrix and the inclusions [11, 12]. In addition, the VIEM is not sensitive to the geometry or concentration of the inclusions. Moreover, in contrast to the finite element method (FEM), where the full domain needs to be discretized, the VIEM requires discretization of the inclusions only [13, 14].

However, the VIEM cannot be directly applied to problems involving voids since the field quantities are undefined in the domain of the integral equation. By contrast, problems involving voids can be most effectively solved using the BIEM. However, if the material contains a combination of elliptical inclusions and a circular/elliptical void, it is most advantageous to apply a combination of the two methods [7–9].

Therefore, in this paper, the mixed volume and boundary integral equation method (mixed VIEM-BIEM) developed by Lee et al. [7–9] is used to calculate the elastic field in composites consisting of an isotropic matrix containing multiple isotropic/anisotropic elliptical inclusions and a circular/elliptical void subjected to remote uniaxial tension. The elliptical inclusions are assumed to have different orientations relative to the loading direction. In order to investigate the influence of a circular/elliptical void on the interfacial stress field, a detailed analysis of the stress field at the interface between the matrix and the central isotropic/orthotropic inclusion is carried out for square packing of eight inclusions and one void, taking into account different values for the orientation angles and concentration of the inclusions. It should be noted that, in the previous mentioned papers [7–9], the inclusions were circles composed of isotropic or orthotropic materials. It is demonstrated that the mixed volume and boundary integral equation method is very accurate and effective in calculating the local stresses for these types of composites.

#### 2. Volume Integral Equation Method

The geometry of the general elastostatic problem considered here is shown in Figure 2 where an unbounded isotropic elastic solid containing a number of isotropic or anisotropic inclusions of arbitrary shape is subjected to prescribed loading at infinity. The symbol denotes the elastic stiffness tensor of the inclusion. denotes the elastic stiffness tensor of the matrix material. The matrix is assumed to be homogeneous and isotropic so that is a constant isotropic tensor, while can be arbitrary; that is, the inclusions may, in general, be inhomogeneous and anisotropic. The interfaces between the inclusions and the matrix are assumed to be perfectly bonded ensuring continuity conditions at the interface.

Let denote the th component of the displacement vector due to the applied remote stress at in the absence of the inclusions. Also, let denote the same in the presence of the inclusions. It has been shown by Mal and Knopoff [15] and Lee and Mal [5] that the elastostatic displacement in the composite satisfies the volume integral equation where the integral is over the domain occupied by the inclusions, , and is the static Green’s function (or Kelvin’s solution) for the unbounded matrix material; that is, represents the th component of the displacement at due to unit concentrated force at in the th direction. In (1), the summation convention and comma notation have been used and the differentiations are with respect to . It should also be pointed out that the integrand is nonzero only within the inclusions, since , outside the inclusions.

If , then (1) is an integrodifferential equation for the unknown displacement vector . It can, therefore in principle, be determined through the solution of the equation. An algorithm for the solution of (1) was developed by Lee and Mal [5] by discretizing the inclusions using conventional finite elements. Once within the inclusions is determined, the displacement field in the matrix can be calculated from (1) by evaluating the integral for . The stress field within and outside the inclusions can also be determined in a similar manner. The details of the numerical treatment of (1) for plane elastostatic problems can be found in [5]. Further explanation of the volume integral equation method for isotropic inclusions in an isotropic matrix can also be found in Section 4.3 “Volume Integral Equation Method” by Buryachenko [16].

#### 3. Mixed Volume and Boundary Integral Equation Method

In this section a mixed volume and boundary integral equation method is introduced as an effective numerical scheme for the solution of elastostatic problems in an unbounded isotropic matrix containing isotropic or anisotropic elliptical inclusions of arbitrary orientation and voids subjected to remote loading. The VIEM is effective for problems with isotropic or anisotropic inclusions due to the fact that it is not necessary to apply continuity conditions at each interface, and only the Green’s function for the unbounded isotropic matrix is needed in the formulation. However, the VIEM cannot be directly applied to problems involving voids, since the field quantities are undefined in the domain of the integral equation. Problems involving voids can be most effectively solved using the BIEM. However, if the material contains a combination of voids and inclusions, it is most advantageous to apply a combination of the two methods. The mixed volume and boundary integral equation for the problem shown in Figure 3 can be expressed in the following form [7, 9] where is the volume of the inclusion, is the surface of the void, and is the outward unit normal to . The superscript implies that the quantities involved are for the matrix. Details of the numerical treatment of (2) for plane elastostatic problems can be found in [7–9].

#### 4. Elastostatic Problems for Multiple Elliptical Inclusions of Arbitrary Orientation and a Void

To investigate the influence of a circular/elliptical void on the interfacial stress field, we consider plane strain problems for multiple isotropic or orthotropic elliptical inclusions of arbitrary orientation and a circular/elliptical void in an unbounded isotropic matrix under uniform remote tensile loading, , as shown in Figure 3. For elliptical cylindrical inclusions or a void, when the major axis coincides with the -axis at an inclined angle , the aspect ratio for each elliptical inclusion or void is assumed to be 0.5 (see Figure 3). is defined as the orientation angle of each elliptical inclusion or void measured counterclockwise from the -axis.

A detailed analysis of the stress field at the interface between the matrix and the central inclusion is carried out for the square packing of eight isotropic/orthotropic inclusions and one void in the unbounded matrix. The fiber volume fraction, defined as that of the circular inclusions (indicated by the dotted lines) circumscribing the elliptical inclusions, is 0.3 and 0.5, respectively (see Figure 3). The inclusions are also rotated from to in increments counterclockwise from the -axis. It should be noted that, since the state of stress in the range of to is a mirror image of that in the range of to , only the orientation angle in the first quadrant is considered here.

The elastic constants for the matrix and the isotropic/orthotropic inclusions are listed in Table 1. Table 2 shows fiber separation distances (i.e., distances between the centers of each elliptical inclusion or circular/elliptical void) for different fiber volume fractions. As an example, square packing leads to a fiber separation distance for , where “” represents the radius of each circular fiber circumscribing each elliptical fiber. It should be noted that while the spacing between the elliptical fibers is independent of their number, aspect ratio, and orientation angle, it does however decrease as the fiber volume fraction increases.

For plane strain problems, the mixed volume and boundary integral equation (2) for multiple orthotropic inclusions and a void can be expressed in the following form: where and are the in-plane displacement components, , , , and , are the stiffness constants of the inclusions, while and are the Lamé constants of the matrix. Also, are given as where is the outward unit normal to .

In (3), represents the Green’s function for the unbounded isotropic matrix material and is given by Banerjee [17] as where and . If the material of the inclusions is isotropic with Lamé constants and , then ), , and .

It should be noted that for uniformly applied remote stresses, the displacement vector is of the form where the constants are related to the tensile and shear components of the applied remote stress.

In (3), is the Green’s function for the unbounded isotropic matrix material and is given in (5). Thus, the mixed volume and boundary integral equation method does not require the use of Green’s function for the anisotropic material of the inclusions. This is in contrast to the BIEM, where the infinite medium of Green’s functions for both the matrix and the inclusions are involved in the formulation of the equations. In general, the Green’s function for an anisotropic material is much more complex than that for isotropic materials [18] and makes the solution of the integrodifferential equation extremely difficult, if not impossible, to carry out. Furthermore, the VIEM is not sensitive to the geometry of the inclusions and in contrast to FEM, where the full domain needs to be discretized, the mixed volume and boundary integral equation method requires discretization of the inclusions (VIEM) and the void (BIEM) only. Specifically, when the fiber volume fraction changes from 0.3 to 0.5, the full domain including the inclusions, the void, and the matrix needs to be discretized from the beginning for each separate case in the FEM model. By contrast, only the position of the elliptical inclusions of arbitrary orientation and the void needs to be changed in the mixed volume and boundary integral equation method.

##### 4.1. Validation of the Mixed VIEM-BIEM Solution for Multiple Elliptical Inclusions and a Void Problem

In order to check the accuracy of the mixed volume and boundary integral equation method, we first compare the numerical solution using the mixed VIEM-BIEM for the square packing of 8 isotropic elliptical inclusions and one elliptical void with an aspect ratio of 0.5, an inclination angle , and a volume concentration of with the solution obtained from the commercial finite element code NASTRAN (2010). The elastic constants for the matrix and the isotropic inclusion are listed in Table 1. Table 2 shows fiber separation distances for different fiber volume fractions.

Figure 4 shows a typical finite element model using NASTRAN, while Figure 5 shows a typical discretized model used in the mixed volume and boundary integral equation method (Patran, 2010). The number of standard eight-node quadrilateral and six-node triangular elements in each isotropic elliptical inclusion was 576, while the number of standard quadratic elements on the surface of the void was 180. Figure 4(a) shows the entire discretized domain, while Figures 4(b), 4(c), and 4(d) show expanded views of the mesh surrounding the elliptical inclusions and the elliptical void (represented as (A)) and the central inclusion (represented as (B)). Standard four-node quadrilateral and three-node triangular elements were used in the discretization. The total number of elements used in the model was 203,200. To capture the accurate stress distribution at the interface between the matrix and the central elliptical inclusion, very refined finite elements were used near the interfaces. The infinite dimensions of the matrix were approximated to be , where “” represents the side length of the square in (A).

It can be easily seen that the discretized model in the mixed volume and boundary integral equation method in Figure 5 is much more effective than the finite element method in Figure 4. This is due to the fact that, in contrast to FEM where the full domain needs to be discretized, the mixed VIEM-BIEM requires discretization of the multiple inclusions and the void only. Thus, for variations in volume fraction and orientation angle, only the position of the inclusions and the void needs to be changed in the mixed VIEM-BIEM model. Figure 6 shows a comparison between the numerical results using mixed VIEM-BIEM and NASTRAN for the normalized von Mises stress component at the interface between the matrix and the central elliptical inclusion ( = 0°~360°).

The von Mises stress can be expressed in the form

For plane strain problems, we have

Thus, the von Mises stress for plane strain problems can be written as or

Figure 6 shows that there is excellent agreement between the mixed VIEM-BIEM and the FEM results. It can be concluded that the mixed volume and boundary integral equation method provides an accurate solution for this class of problems.

##### 4.2. Square Packing of Eight Elliptical Inclusions and One Circular/Elliptical Void at an Inclined Angle

Next, we examine the square packing of 8 isotropic/orthotropic elliptical inclusions and one circular/elliptical void with an aspect ratio of 0.5, at an inclination angle . The fiber volume fraction is 0.3 and 0.5, respectively. It should be noted that the fiber volume fraction of the elliptical inclusions is defined as that of the circular inclusions (indicated by the dotted lines) circumscribing the elliptical inclusions (see Figure 3). The elastic constants for the matrix and the isotropic/orthotropic inclusions are listed in Table 1. Table 2 shows fiber separation distances (i.e., distances between the centers of each elliptical inclusion or circular/elliptical void) for different fiber volume fractions.

Figure 7 shows a typical discretized model of 8 orthotropic elliptical inclusions and one circular void at an inclined angle used in the mixed VIEM-BIEM [19]. The number of standard eight-node quadrilateral and six-node triangular elements in each orthotropic elliptical inclusion was 576, while the number of standard quadratic elements on the surface of the circular void was 180. Figure 8 shows the normalized tensile stress component () at the interface between the matrix and the central elliptical inclusion for models (a) containing a square array of 8 isotropic elliptical inclusions and one circular/elliptical void and (b) containing a square array of 8 orthotropic elliptical inclusions and one circular/elliptical void. For these cases, the fiber volume fraction is 0.3 and 0.5, where the orientation angle () is ( = 0°~360°). Here, indicates the angle of the central elliptical inclusion at the matrix-inclusion interface as measured counterclockwise from the -axis. Figure 8 indicates that the interfacial stress of the central elliptical inclusion with an elliptical void is greater than that of the central elliptical inclusion with a circular void, and the interfacial stress of the central isotropic elliptical inclusion is greater than that of the central orthotropic elliptical inclusion, since in the isotropic inclusion is greater than that in the orthotropic inclusion. It should also be noted that the stress at the interface between the matrix and the central isotropic inclusion appears to increase as the fiber volume fraction increases, while the stress at the interface between the matrix and the central orthotropic inclusion appears to decrease as the fiber volume fraction increases.

##### 4.3. Square Packing of Eight Elliptical Inclusions and One Circular/Elliptical Void at an Inclined Angle

Next, we examine a square packing of 8 isotropic/orthotropic elliptical inclusions and one circular/elliptical void with an aspect ratio of 0.5, at an inclination angle . The fiber volume fraction is 0.3 and 0.5, respectively. The elastic constants for the matrix and the isotropic/orthotropic inclusions are listed in Table 1. Table 2 shows fiber separation distances (i.e., distances between the centers of each elliptical inclusion or circular/elliptical void) for different fiber volume fractions.

Figure 9 shows the normalized tensile stress component () at the interface between the matrix and the central elliptical inclusion for models (a) containing a square array of 8 isotropic elliptical inclusions and one circular/elliptical void and (b) containing a square array of 8 orthotropic elliptical inclusions and one circular/elliptical void. For these cases, the fiber volume fraction is 0.3 and 0.5, where the orientation angle () is 30° ( = 0°~360°). It should be noted that the interfacial stress of the central elliptical inclusion with an elliptical void is greater than that of the central elliptical inclusion with a circular void, the stress at the interface between the matrix and the central inclusion appears to decrease as the fiber volume fraction increases, and the interfacial stress of the central isotropic elliptical inclusion is greater than that of the central orthotropic elliptical inclusion, since in the isotropic inclusion is greater than that in the orthotropic inclusion.

##### 4.4. Square Packing of Eight Elliptical Inclusions and One Circular/Elliptical Void at an Inclined Angle

Figure 10 shows the normalized tensile stress component at the interface between the matrix and the central elliptical inclusion for models (a) containing a square array of 8 isotropic elliptical inclusions and one circular/elliptical void and (b) containing a square array of 8 orthotropic elliptical inclusions and one circular/elliptical void. For these cases, the fiber volume fraction is 0.3 and 0.5, where the orientation angle () is 60° ( = 0°~360°). It should be noted that the general features of the interfacial stress are similar to those observed in the cases where the inclined angle , as discussed in the preceding subsection.

##### 4.5. Square Packing of Eight Elliptical Inclusions and One Circular/Elliptical Void at an Inclined Angle

Figure 11 shows the normalized tensile stress component () at the interface between the matrix and the central elliptical inclusion for models (a) containing a square array of 8 isotropic elliptical inclusions and one circular/elliptical void and (b) containing a square array of 8 orthotropic elliptical inclusions and one circular/elliptical void. For these cases, the fiber volume fraction is 0.3 and 0.5, where the orientation angle () is 90° ( = 0°~360°). It should be noted that the general features of the interfacial stress are similar to those observed in the cases where the inclined angle , as discussed in the preceding subsection.

Figures 8, 9, 10, and 11 show that the stress at the interface between the matrix and the central inclusion appears to decrease as the orientation angle increases for all cases considered. Given these results, it is therefore necessary to model arbitrarily oriented isotropic or anisotropic fibers with arbitrary geometry and various void shapes in order to accurately predict the failure and damage mechanisms in real composites.

#### 5. Concluding Remarks

Cross sections of certain composites have been shown to contain fibers of near elliptical shape with varying orientation angles. If voids are caused by manufacturing and/or service induced defects, they may have a significant effect on the local stresses in a composite. For this reason, the stress field at the interface between the matrix and the central isotropic/orthotropic inclusion was investigated through the square packing of eight isotropic/orthotropic elliptical inclusions of arbitrary orientation and one circular/elliptical void with an aspect ratio of 0.5 under uniform remote tensile loading, taking into account different values for the orientation angles and concentration of the inclusions. It was observed that the stress at the interface between the matrix and the central inclusion decreased as the orientation angle increased for all cases considered, the interfacial stress of the central elliptical inclusion with an elliptical void was greater than that of the central elliptical inclusion with a circular void, and the interfacial stress of the central isotropic elliptical inclusion was greater than that of the central orthotropic elliptical inclusion, since in the isotropic inclusion was greater than that in the orthotropic inclusion. In addition, it was also discovered that as the fiber volume fraction increased, the interaction effect became weaker except for the square packing of eight isotropic elliptical inclusions and one circular/elliptical void at an inclined angle .

Through the analysis of plane elastostatic problems in an unbounded isotropic matrix with multiple isotropic or anisotropic elliptical inclusions of arbitrary orientation and a circular/elliptical void under uniform remote tensile loading, it was further noted that in order to accurately predict the failure and damage mechanisms in real composites, it is necessary to investigate multiple arbitrarily oriented isotropic or anisotropic fibers with arbitrary geometry and various void shapes under arbitrary loading conditions [20, 21].

#### Acknowledgments

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no. 2012-0007636). The authors would also like to acknowledge the support from Korea Institute of Science and Technology Information (KISTI) supercomputing center through the strategic support program for the supercomputing application research (no. KSC-2012-C1-06) and appreciate their assistance in parallelizing and optimizing the VIEM code cited in this paper.

#### References

- D. S. Brauer, C. Rüssel, S. Vogt, J. Weisser, and M. Schnabelrauch, “Degradable phosphate glass fiber reinforced polymer matrices: mechanical properties and cell response,”
*Journal of Materials Science: Materials in Medicine*, vol. 19, no. 1, pp. 121–127, 2008. View at Publisher · View at Google Scholar · View at Scopus - R. M. Christensen and K. H. Lo, “Solutions for effective shear properties in three phase sphere and cylinder models,”
*Journal of the Mechanics and Physics of Solids*, vol. 27, no. 4, pp. 315–330, 1979. View at Scopus - B. P. Maslov, “Stress-strain state in the matrix of stochastically reinforced composites,”
*Mechanics of Composite Materials*, vol. 25, no. 3, pp. 282–288, 1989. View at Publisher · View at Google Scholar · View at Scopus - R.-B. Yang and A. K. Mal, “The effective transverse moduli of a composite with degraded fiber-matrix interfaces,”
*International Journal of Engineering Science*, vol. 33, no. 11, pp. 1623–1632, 1995. View at Scopus - J. Lee and A. Mal, “A volume integral equation technique for multiple inclusion and crack interaction problems,”
*Journal of Applied Mechanics*, vol. 64, no. 1, pp. 23–30, 1997. View at Scopus - V. A. Buryachenko and V. T. Bechel, “A series solution of the volume integral equation for multiple-inclusion interaction problems,”
*Composites Science and Technology*, vol. 60, no. 12-13, pp. 2465–2469, 2000. View at Publisher · View at Google Scholar · View at Scopus - J. Lee, S. Choi, and A. Mal, “Stress analysis of an unbounded elastic solid with orthoropic inclusions and voids using a new integral equation technique,”
*International Journal of Solids and Structures*, vol. 38, no. 16, pp. 2789–2802, 2001. View at Publisher · View at Google Scholar · View at Scopus - J. Lee, H. Lee, and A. Mal, “A mixed volume and boundary integral equation technique for elastic wave field calculations in heterogeneous materials,”
*Wave Motion*, vol. 39, no. 1, pp. 1–19, 2004. View at Publisher · View at Google Scholar · View at Scopus - J. Lee, H. Han, and A. Mal, “Effects of anisotropic fiber packing on stresses in composites,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 195, no. 33–36, pp. 4544–4556, 2006. View at Publisher · View at Google Scholar · View at Scopus - J. Lee and H.-R. Kim, “Volume integral equation method for multiple circular and elliptical inclusion problems in antiplane elastostatics,”
*Composites: Part B*, vol. 43, no. 3, pp. 1224–1243, 2012. View at Publisher · View at Google Scholar · View at Scopus - S. G. Mogilevskaya and S. L. Crouch, “A Galerkin boundary integral method for multiple circular elastic inclusions with homogeneously imperfect interfaces,”
*International Journal of Solids and Structures*, vol. 39, no. 18, pp. 4723–4746, 2002. View at Publisher · View at Google Scholar · View at Scopus - F.-Z. Kong, X.-P. Zheng, and Z.-H. Yao, “Numerical simulation of 2D fiber-reinforced composites using boundary element method,”
*Applied Mathematics and Mechanics*, vol. 26, no. 11, pp. 1515–1522, 2005. View at Scopus - M.-C. Chen and X.-C. Ping, “A novel hybrid finite element analysis of inplane singular elastic field around inclusion corners in elastic media,”
*International Journal of Solids and Structures*, vol. 46, no. 13, pp. 2527–2538, 2009. View at Publisher · View at Google Scholar · View at Scopus - H. H. Zhang and L. X. Li, “Modeling inclusion problems in viscoelastic materials with the extended finite element method,”
*Finite Elements in Analysis and Design*, vol. 45, no. 10, pp. 721–729, 2009. View at Publisher · View at Google Scholar · View at Scopus - A. K. Mal and L. Knopoff, “Elastic wave velocities in two-component systems,”
*Journal of the Institute of Mathematics and Its Applications*, vol. 3, no. 4, pp. 376–387, 1967. View at Publisher · View at Google Scholar · View at Scopus - V. A. Buryachenko,
*Micromechanics of Heterogeneous Materials*, Springer, New York, NY, USA, 2007. - P. K. Banerjee,
*The Boundary Element Methods in Engineering*, McGraw-Hill, New York, NY, USA, 1993. - K. J. Lee and A. K. Mal, “A boundary element method for plane anisotropic elastic media,”
*Journal of Applied Mechanics*, vol. 57, no. 3, pp. 600–606, 1990. View at Scopus - PATRAN User's Manual, Version 7. 0, MSC/PATRAN, 1998.
- J. A. Koltay and D. L. Feke, “Preparation of continuous fiber ceramic composites using a combination of steric-stabilization and depletion-flocculation phenomena,”
*Composites A*, vol. 30, no. 3, pp. 231–237, 1999. View at Scopus - Y. T. Lee and J. T. Chen, “Null-field approach for the antiplane problem with elliptical holes and/or inclusions,”
*Composites: Part B*, vol. 44, no. 1, pp. 283–294, 2013. View at Publisher · View at Google Scholar