A microstructure modeling scheme for unidirectional composites using signed distance function based boundary smoothing and element trimming

Advances in Engineering Software 109 (2017) 1–14

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

Research paper

A microstructure modeling scheme for unidirectional composites using signed distance function based boundary smoothing and element trimming Jae Hyuk Lim a, Hobeom Kim b, Sun-Won Kim c, Dongwoo Sohn d,∗ a

Department of Mechanical Engineering, College of Engineering, Chonbuk National University, 567 Baekje-daero, Deokjin-gu, Jeonju, Jeonbuk 54896, Republic of Korea Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea c Satellite Mechanical Department, Korea Aerospace Research Institute, 169-84 Gwahak-ro, Yuseong-gu, Daejeon 34133, Republic of Korea d Division of Mechanical Engineering, College of Engineering, Korea Maritime and Ocean University, 727 Taejong-ro, Yeongdo-gu, Busan 49112, Republic of Korea b

a r t i c l e

i n f o

Article history: Received 11 August 2016 Revised 19 January 2017 Accepted 12 February 2017 Available online 22 March 2017 Keywords: Trimmed meshes Image processing Signed distance function Unidirectional composites

a b s t r a c t A simple and accurate scheme for modeling microstructures is proposed with the help of element trimming combined with signed distance function based boundary smoothing. To accommodate randomly distributed fibers in unidirectional composites, digital image processing is used. The interfaces of multimaterials are identified by introducing a signed distance function, and then, square background elements crossing the interfaces are simply trimmed and divided to represent a single material behavior by a single element. After element trimming, the elements that are polygon-shaped in the two-dimensional domain are split into conventional three-node triangle elements (six-node prism elements in the threedimensional domain) available in many commercial software packages. The present modeling scheme was verified through benchmark examples in terms of the accuracy and efficiency and then applied to the modeling of unidirectional composites based on real microscopic images to evaluate the equivalent elastic properties. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction A variety of composite materials (reinforced plastics, metal composites, ceramic composites and so forth) have been widely adopted in various applications in the automotive, offshore plant, construction, sporting, and aerospace industries due to their high stiffness and strength-to-weight ratio. Unidirectional (UD) fiberreinforced polymer composites containing microstructures have been routinely modeled by two constituents (the fiber and matrix) or by more than two constituents (fiber, matrix, void and interphase/interface). The microstructure conditions at the small-scale level (constituent-level), such as dimensions, shapes, spatial distributions, material properties of the constituents, strongly affect the interactions between the constituents and thus, have a critical role in the performance at the large-scale level (lamina-, laminate, and structure-levels) [1,2]. Therefore, there have been numerous

∗

Corresponding author. E-mail address: [email protected] (D. Sohn).

http://dx.doi.org/10.1016/j.advengsoft.2017.02.014 0965-9978/© 2017 Elsevier Ltd. All rights reserved.

efforts to understand the link of the mechanical response between the small-scale and large-scale levels of UD composites. As part of such efforts, to reflect the microstructure effects, many researchers have developed finite element (FE) model generators [3,4] to consider the statistical distribution of the equivalent elastic properties of the lamina and investigated the effect of artificially generated microstructures. They thus succeeded in a more accurate prediction higher than the FE solutions of simple hexagonal, square fiber arrays or analytic solutions using the Mori-Tanaka method with the assumption of transverse isotropy [5]. However, in a practical sense, manufacturing composite materials that have statistical fiber distributions is not a straightforward process because the microstructure configurations and material properties of composites are quite affected by the manufacturing conditions such as the curing cycle, tool-part interaction during the curing process, level of applied vacuum pressure, temperature control and so forth [5]. Therefore, artificially generated fibers in numerical modeling are somewhat different from real fibers not to mention that their shapes are not purely circular and their diameters are not constant. To accurately investigate the mechanical behavior of UD composites, more realistic FE models that reflect the

2

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

microscope images of SEM or TEM are required as presented by many researchers [6–13]. To accommodate a real microscope image directly as it is, voxel and voxel-like methods have been introduced [6]. Because these methods treat one pixel to be one square FE or four triangle FEs, a large number of nodes have to be used to obtain accurate results for complex objects. It is thus quite difficult to cover a large-size domain due to limited computer resources. In addition, generally, they generate jagged edges that can cause an undesirable stress concentration in contrast to smooth boundaries in real materials. Furthermore, continuing efforts to generate FE meshes that accommodate microstructures through image processing have been reported [7–13]. On the other hand, as a conventional method, it is possible to create a FE model directly when the CAD file of the microstructure is already arranged. However, creating CAD files of complex objects could still be a cumbersome task even though there have been many dramatic enhancements in commercial CAD software, and a few successes in converting images to CAD files have been reported [14–16]. From the viewpoint of computational mechanics, to consider the detailed configuration of microstructures efficiently, many multi-material modeling schemes have been proposed by several researchers using the extended finite element method (XFEM) [17,18], mesh-free methods [19,20], generalized finite element method (G-FEM) [21], and the element carving/trimming scheme [22]. In all of them, basically, a background mesh with square elements is used to cover the entire domain of problems. To identify each material phase on the background mesh, the level-set method [23] has been frequently used. In the level-set method, object boundaries are represented implicitly using the minimum signed distance function (SDF) values as nodal variables. Especially, when the level-set method is combined with the aforementioned multi-material modeling scheme such as the XFEM, it provides an efficient modeling tool for dealing with various singularities including cracks, multi-material interfaces, dislocations and so forth without remeshing. It also allows one element to cover multi-material constituents with a subdomain (distinguished by the signed distance)-wise Gauss integration. In contrast, in the element carving/trimming scheme, the elements crossing multi-material interfaces are locally trimmed so that each element cover only one material. In other words, this scheme splits the square elements in the background mesh to several polygonshaped pieces. In particular, when the trimmed square elements, which are polygonal-shaped, is further divided into three-node triangle elements in the two-dimensional (2D) domain (six-node prism elements in the three-dimensional (3D) domain), FE models obtained by this scheme can be directly used in FE analysis by means of any commercial FE software such as ABAQUS. The goal of this work was to propose a simple and accurate multi-material modeling scheme using the SDF based boundary smoothing and element trimming techniques. The remainder of this paper is organized as follows. In Section 2, we explain how to construct trimmed FE meshes to fit multi-material interfaces represented by the SDF. This is followed by an element merging scheme to enhance the mesh quality. Some distorted and small-area elements, which can be generated during the element trimming, should be merged into large elements to avoid the ill-conditioning or nearly zero values of the stiffness matrix. This manipulation can thus provide accurate solutions comparable to those of conventional FE meshes with CAD files. Next, in Section 3, to verify the efficiency and accuracy of the proposed scheme compared with other numerical schemes, we solve benchmark problems regarding plates containing a hole or an inclusion. In Section 4, we show the effectiveness of the proposed scheme in the FE modeling of UD composites with complex fiber configurations, which are acquired from image processing. Finally, we finish the paper with concluding remarks in Section 5.

2. Multi-material representation using the SDF based boundary smoothing and element trimming 2.1. SDF based boundary smoothing algorithm Consider an analysis domain covered by square-background elements. Fig. 1 shows two material regions, 1 and 2 , and their interface, ∂  = 1 ∩2 . If a few objects are located in the domain, the material interface is represented using the minimum SDF. For any material point x in the domain at time t, the minimum SDF φ (x,t) is defined as the minimum distance between the nodes of the square elements and the material interfaces which is a scalar variable independently defined at the nodes. The zero SDF value (φ (x,t) = 0) corresponds to the interface of the two materials. For the case of multi-materials (n-materials), the independent variables, for which the number is (n − 1), have to be considered at the nodes. With this concept, the X-FEM has reported remarkable successes for many singularity problems [17,18] without additional local remeshing because it describes the material interfaces implicitly using the SDF. However, in our approach, we explicitly split the square elements crossing multi-material boundaries into several polygon-shaped pieces (see Fig. 2(a)) which are represented by the combination of simple three-node triangle elements in the 2D domain (six-node prism elements in the 3D domain) shown in Fig. 2(b). As a result, using element trimming based on the SDF values, all of the material boundaries have a smoothed representation. Before calculating the SDF values, a high-quality image should be first prepared. To obtain an image with local contrast enhancement and noise reduction, we conduct a preprocessing using two well-known image filters; a median filter taking the median value to remove speckles/dots on an image, and a Gaussian high-pass filter to sharpen an image. The two filters are supported by MATLAB image processing toolbox [24]. Subsequently, for a given black-andwhite patterned image passed through the filters, we consider the following two techniques to define the minimum SDF for detecting material boundaries. The first way of obtaining SDF values from raw images is to solve iteratively a level-set equation given by Eq. (1) on squaretype background meshes, until the solution of the level-set equation is converged to constant values within a tolerance limit [23,25]:

∂φ (x, t )  ∇φ = 0, +ν ∂t

(1)

φ (x, t ) < 0 for x ∈ 1 ,

(2)

φ (x, t ) > 0 for x ∈ 2 ,

(3)

φ (x, t ) = 0 for x ∈ ∂ .

(4)

where the minimum SDF φ should satisfy the following properties;

 is taken, Eq. (1) beWhen the normal component ν N of velocity ν comes the following partial differential equation:

∂φ (x, y ) + νN |∇φ| = 0. ∂t

(5)

Considering the curvature-based level-set evolution, or mean curvature flow, with a curvature-based force that smoothes the curve, ν N = −bκ , Eq. (5) can be rewritten as [25]

∂φ (x, y ) = bκ|∇φ|, ∂t

(6)

where κ is the curvature, and b is a weighting parameter for the curvature-based force. Eq. (6) can be discretized using central differencing with a uniform Cartesian grid. In this work, we use an

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

3

Fig. 1. A multi-material domain represented by the minimum SDF.

shown in Fig. 3(a) and (b) which has been widely used in the image processing field and has exhibited excellent performance [28]. In this work, we use this approach to extract material interfaces. In the case of a bulk object that has a non-zero area, the output after the edge-detection is simply a set of nodes describing a closedloop polygon. Its line segments are generally not smooth but somewhat jagged due to the pixel resolution shown in Fig. 3(c). To increase the smoothness of the polygon, we apply a localized cubic spline fitting with a set of neighboring nodes of the polygon and then find the magnitudes of the minimum distances between the FE nodes and the fitting curve. To assign the sign of the distance values, the crossing number algorithm (also termed the ray casting algorithm [29]) is used. We can thus define the SDF values as having more accurate and smooth interface information. Subsequently, we explicitly define a position X(x, y) of the intersection point of the material boundary and the element edge using the minimum SDF values, φ 1 and φ 2 , at two nodes associated with the element edge as follows:



X=

 |φ1 |x2 + |φ2 |x1 |φ1 |y2 + |φ1 |y1 , , |φ1 | + |φ2 | |φ1 | + |φ2 |

(7)

where φ 1 × φ 2 < 0 is satisfied. Finally, based on such an intersection point X(x, y), a square background element is trimmed, and the trimmed elements are further split into several triangle elements. 2.2. Treatment of bad-shaped elements Fig. 2. Element trimming of multi-materials: (a) a few subcases of element trimming; and (b) reconstruction of a polygonal region by triangle elements.

open MATLAB toolbox, which Sumengen [26] developed following the Osher and Fedkiw’s book [25], to obtain the minimum SDF values. However, allowing a tolerance in solving Eq. (6) may lead to the loss of geometric accuracy. Generally, it also requires considerable computational cost for large-size raw images, although the cost may be reduced by the resampling technique of the raw image [25]. Besides this, it is well known that a traditional form of solving the level-set equation based on the finite difference method cannot detect sharp corners and cannot preserve the volume well although several solutions have been reported [27]. In solid mechanics problems, such geometric deviation is not acceptable because sharp corners function as sources of stress singularities. The second approach is based on the canny algorithm to detect object or material boundaries of black-and-white patterned images

As mentioned previously, to accommodate multi-materials explicitly, we split one square element into several triangle elements according to the element-crossing cases shown in Fig. 2. However, this element subdivision could generate tiny or skewed elements that lead to an ill-conditioned stiffness matrix or a nearly-zero stiffness term. In order to simply avoid such a situation, we merge such small elements to their neighboring large elements, shown in Fig. 4, keeping the nodes with zero SDF values untouched. Note that the material boundaries are not changed even though the configurations of some elements are modified. This modification assists in slightly improving the accuracy of FE solutions, which will be discussed through numerical tests in Section 3. 3. Verification of the proposed multi-material modeling scheme To demonstrate the performance of the proposed scheme, we solved two benchmark problems under tensile conditions: an infi-

4

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

Fig. 3. Image processing of microstructures dispersed in a UD composite: (a) raw image; (b) raw image with detected edges by the canny algorithm; (c) a smoothed interface by spline interpolation; and (d) computed SDF values.

Fig. 4. Two examples of element merging strategies.

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

5

alize the infinite body state into a finite domain, we calculated the exact nodal forces by integrating the exact stress distributions for the infinite body, given by Eq. (9), with the 12th order Gauss integration along the boundary, and imposed these forces on the outer boundary [18]:

 

 3 a4 a2 3 σ11 (r, θ ) = σ0 1 − 2 cos 2θ + cos 4θ + cos 4 θ , 2 r4 r 2   2

 3 a4 a 1 σ22 (r, θ ) = σ0 − 2 cos 2θ − cos 4θ − cos 4θ , 2 r4 r 2  2

  3 a4 a 1 σ12 (r, θ ) = σ0 − 2 sin 2θ + sin 4θ + sin 4θ . (9) 4 r

2

2r

where u is a displacement vector; ε is a strain matrix and C is an elastic modulus matrix. Eq. (8) is calculated by the summation of the integral over an individual element domain i , where the number of elements is ne . The superscripts “exact” and “fe” indicate the exact solution and the approximate FE solution, respectively.

To handle this problem, four FE solution methods were prepared: (a) conventional conforming meshes based on CAD files; (b) trimmed FE meshes by the proposed scheme shown in Fig. 6; (c) X-FEM solution in [18] and (d) voxel FE meshes. Except for the XFEM solution, we obtained results using the commercial software ABAQUS [30]. When the voxel meshes were used, the geometric centers of the elements were used to assign the material properties of every element. To check the convergence rate of the solution by the mesh trimming, we prepared four square-type background meshes with the following resolutions: 10 × 10, 20 × 20, 40 × 40, and 80 × 80 elements. The relative errors in strain energy norm were plotted in Fig. 7 against the mesh size h which indicates the average characteristic length of the used elements. The overall convergence rate using the trimmed FE meshes by the proposed scheme (1.0) is similar to those of the conventional conforming mesh solution (0.96) and X-FEM solution (0.96), not to mention that it is much higher than that of the voxel FE mesh solution (0.81). From these results, it can be deduced that the proposed scheme provides accurate FE solutions with a comparable convergence rate to the conventional FE scheme.

3.1. An infinite plate containing a hole

3.2. A circular plate containing a circular inclusion

Consider an infinite plate with a hole subjected to a tensile load, σ 0 = 1 Pa, at the far-field boundary in the X1 -direction shown in Fig. 5. Under the plane strain condition, we dealt with a block domain for which the dimension is 2 m × 2 m, containing a hole with a radius of 0.3 m. Then, we only used a quarter model by imposing the proper symmetric boundary conditions. The Young’s modulus and Poisson’s ratio are 106 Pa and 0.3, respectively. To re-

Consider a circular plate containing a circular inclusion subjected to radial tension shown in Fig. 8. A material discontinuity exists across the interface 2 ( r = a). The plate under the plane strain condition is subjected to uniform displacements on the boundary 2 ( r = b), given by u1 = x1 and u2 = x2 (i.e., ur = r and uθ = 0). The Young’s modulus and Poisson’s ratio of the plate are 10 Pa and 0.30, respectively. Those of the circular inclusion are set

Fig. 5. A square domain containing a circular hole under horizontal tension.

nite plate containing a hole, and a circular plate containing a circular inclusion. As an error measure for the benchmark problems, the relative error in strain energy norm was used as follows:

⎛

Ee =

⎝

ne  

⎞1/2

(ε f e − ε exact ) C(ε f e − ε exact )di ⎠ T

i  i

⎛ ×⎝

ne  

/

⎞1/2 (ε exact ) C(ε exact )di ⎠ T

,

(8)

i  i

Fig. 6. Trimmed FE mesh construction of the problem of the infinite plate containing a hole.

6

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

Fig. 7. Convergence rate in terms of the relative strain energy error norm, with respect to different numerical schemes applied to the problem in Fig. 5.

Fig. 9. Convergence rate in terms of the relative strain energy error norm, with respect to different numerical schemes applied to the problem in Fig. 8.

In Eqs. (11) and (12), the parameter α is defined as

α=

Fig. 8. A circular plate containing a circular inclusion under radial extension.

as 1.0 Pa and 0.25, respectively. In this example, the radii of circular inclusion and plate are set as a = 0.4 m and b = 2.0 m, respectively. Similar to the previous example, we considered a quarter model of the 2 m × 2 m block containing the circular inclusion and applied the appropriate nodal forces obtained by integrating the exact stress distributions to the boundaries of the square domain. The exact stresses are given as follows:

σrr (r, θ ) = 2μεrr + λ(εrr + εθ θ ), σθ θ (r, θ ) = 2μεθ θ + λ(εrr + εθ θ ), σrθ (r, θ ) = 0,

(10)

where the Lamé constants λ and μ are assigned appropriately for the plate and inclusion. The exact strains are



εrr (r, θ ) =

1−



b2 a2



α+

b2 , a2

0≤r≤a



,

α − br2 , a ≤ r ≤ b   2 2 1 − ab2 α + ab2 , 0 ≤ r ≤ a εθ θ (r, θ ) = ,  2 2 1 − br2 α + br2 , a ≤ r ≤ b 1+

b2 r2

2

εrθ (r, θ ) = 0,

(11)

and the exact displacements are



ur (r, θ ) =

uθ (r, θ ) = 0.

1−

r−

b2 a2

b2 r





α+

α+

b2 a2

b2 r





r, 0 ≤ r ≤ a

, a≤r≤b

, (12)

(λ1 + μ1 + μ2 ) , (λ2 + μ2 )a2 + (λ1 + μ1 )(b2 − a2 ) + μ2 b2

(13)

where the subscripts, “1” and “2” indicate the inside and outside domains, respectively. In this example as well, we obtained FE results using the four solution methods described in the previous example. The relative errors in strain energy norm were evaluated varying the mesh size h, and their convergence rates were plotted in Fig. 9. The overall convergence rates by the proposed scheme were 1.103 and 1.106, which were obtained using the trimmed meshes and the trimmed meshes with bad-shaped element merging, respectively. These values are slightly higher than that of conventional conforming meshes (1.04), as well as that of the X-FEM (0.90) [18], not to mention that they are much higher when compared to using the voxel (0.38). The major reason that the proposed scheme has a better solution accuracy might be that we put more actual nodes directly along the material interfaces, which can represent the material discontinuity well. To examine the accuracy of the stress solution in a more strict sense, we compared the radial stress and displacement distributions (σ rr and ur ) along the material interface. First, we obtained all stress and displacement components in the Cartesian coordinates along the material interface and then, transformed those to the radial components. Fig. 10 shows the radial stress and displacement distributions along the material interface and von Mises stress distributions over the analysis domain, varying the mesh size h. The maximum discrepancies in the radial displacements between the analytic solution and the FE solutions, which were obtained using the conventional conforming meshes, the trimmed meshes, and the trimmed meshes with bad-shaped element merging, are almost identical regardless of the mesh types and element sizes, and they are also negligibly small. In terms of the stresses, the trimmed meshes show a little fluctuation, about 1.5% in the case of h = 0.05, which however decreases gradually as the mesh is refined. From the results in Fig. 10, it was confirmed that the stress obtained by the proposed scheme shows good agreement with the analytic solution as well as the conventional conforming mesh solution within an acceptable error level. These slight errors in the stress components might be caused by a stress recovery algorithm based on the element configurations in ABAQUS because the trimmed meshes do not have regular meshes along

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

7

Fig. 10. Comparison of radial stress and displacement distributions along the inclusion boundary and von Mises stress distributions over the analysis domain: (a) h = 0.05; (b) h = 0.025; and (c) h = 0.0125.

8

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

Fig. 10. Continued

the bi-material interface in contrast to the conventional conforming meshes. The nodal stresses are obtained averaging the values extrapolated from Gauss integration points of several elements that share a node; to do this, we use a stress output option of “averaged at nodes” supported in ABAQUS. In addition, the effect of the badshaped element merging does not seem to be significant in this example; however, tiny or skewed elements should be eliminated because they may cause ill-conditioning problems in the stiffness matrix and may sometimes make FE analysis impossible. 4. Modeling of UD composites with image processing In this section, we deal with FE modeling of UD composites, as practical applications of the scheme verified in the previous sections. In particular, we perform image processing, using the MATLAB image processing toolbox [24], to accommodate the irregular configurations of microstructures in UD composites. First, we created artificial microstructure images with the assumption that perfectly circular fibers with a constant radius maintain their proper distance from each other and then, constructed FE meshes using the proposed scheme. Second, from a raw image obtained by optical microscope, we modeled real fiber arrangement without any assumptions on the shape, size, and distance of the fibers. Finally, through computational homogenization, we discuss the results obtained from the artificial and real images. In this section, all numerical results by the proposed scheme were obtained with ABAQUS. 4.1. Random fiber models from artificially created images To investigate the relationship of the microscopic constituents and macroscopic responses of composites and to demonstrate the effectiveness of the proposed scheme, we first dealt with random fiber arrangements that were artificially generated. Here, blackand-white images of 100 representative volume elements (RVEs)

were prepared with random fiber arrangements. It has been widely accepted that such random fiber models are more realistic and accurate in evaluating elastic properties by computational homogenization than simple periodic models containing one or two fibers [5]. To determine the fiber positions in the RVEs, random fiber generation was iteratively done with the constraint that overlapping of fibers was not allowed as in [3,4]. Each RVE image included 50 circular fibers in 2D domains, which is necessary for statistically evaluating the equivalent elastic properties of composites as verified in [31]. To construct conforming meshes in an efficient manner, we calculated the minimum SDF for all artificially created random fiber images and defined material interfaces with the aforementioned scheme in Section 2.1 implemented with the MATLAB image processing toolbox. Then, we trimmed the square background meshes across the material interfaces. For example, as shown in Fig. 11(a), a 2D background mesh with square elements was prepared, and the square elements crossing the material interfaces were then split into several triangle elements so that 2D conforming meshes were obtained. Subsequently, by dragging the 2D meshes along the fiber longitudinal direction, 3D FE models for UD composites were constructed with six-node prism elements and eight-node hexahedral elements shown in Fig. 11(b). There is no configuration change in the fiber longitudinal or thickness direction. Therefore, as long as the number of element layers in the thickness direction is two or more in order to impose the periodic boundary conditions (PBCs) on all three pairs of RVE boundaries, thickness does not affect the evaluation of equivalent elastic moduli. For an efficient analysis, we considered only two layers in the thickness direction, which corresponds to the 1-direction in Section 4. In addition, the element length in the 1-direction was chosen to be the same as the lengths of a 2D background element in the 2- and 3-directions. In this example, because periodic fiber arrangements in the RVEs were considered, all FE models of the RVEs were represented

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

9

Fig. 11. An artificially created random fiber image and its corresponding trimmed FE mesh: (a) 2D black-and-white image and trimmed FE mesh; and (b) magnified view of a 3D trimmed FE mesh.

by periodic meshes that have compatible nodes and element-faces for each pair of exterior boundary surfaces. Therefore, the PBCs can be directly imposed at the nodes on the boundary surfaces of the RVE. For 100 RVEs, 200 × 200 square elements in the 2D domains were used as a background mesh for element trimming, and the fiber-volume fraction was controlled to be 68.33%. By applying the PBCs on the RVEs, we evaluated the equivalent elastic properties of 100 RVE samples which will be discussed in Section 4.3. 4.2. Random fiber models from real images In contrast to the artificially created fibers in the previous subsection, real fibers are not perfectly circular in 2D domains, and their sizes are constant. In addition, fibers can be located directly contacting neighboring fibers, without a proper gap between the fibers. Because the real configurations of fibers can affect the overall response of composite structures, direct FE modeling from image processing is necessary. As seen in Fig. 12, we prepared a microscope image of M55J/M18 UD composites taken with an Axioimager 2 [32] using the proper surface treatment which consisted of 740 × 758 pixels. After converting the colored images to black-

and-white patterned images with local contrast enhancement and noise removing, the signed minimum distances were defined at the four nodes of each square background element. Then, the background elements that covered the entire domain were trimmed along the material boundaries with zero SDF values. Finally, the conforming meshes, which can represent the material interfaces well, were obtained without any assumptions on the shape and size of the fibers. In this example, because extremely tiny or skewed elements are unavoidably generated between the fibers very close to each other, the bad-shaped elements should be eliminated as described in Section 2.2. It is noted that in this example, the bad-shaped elements having a nearly zero stiffness term lead to the ill-conditioning problem, and thus, FE analysis itself fails without removing the bad-shaped elements. Here, we took two samples, Sample #1 (A-B-C-D) and Sample #2 (A-E-F-G), from the microscope image shown in Fig. 13. Using the SDF based boundary smoothing and element trimming schemes, we obtained two trimmed FE meshes shown in Fig. 13 which have similar fiber-volume fractions (Vf ), 68.38% and 68.33%, respectively. By dragging the 2D trimmed meshes along the fiber longitudinal direction with two element layers,

10

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

Fig. 12. Conversion procedure of real microscope images to black-and-white patterned images.

Fig. 13. Real microscope image and trimmed FE meshes over two sample regions.

Fig. 14. Three simple regular RVEs: (a) a square RVE; (b) a hexagonal RVE; and (c) a diamond RVE.

the 3D trimmed meshes were obtained. For Sample #1, based on 253 × 215 background elements in the 2D domain, 186,168 nodes were used in the 3D model. In the case of the Sample #2 based on 301 × 294 background elements in the 2D domain, the number of nodes in the 3D model is 302,385.

If we used conventional 3D voxel discretization which regards a square pixel in image processing as a square element, the number of nodes for Sample #2 should be 1687,257 (741 × 759 × 3), which is about five times more than that of the Sample #2 mesh obtained by the proposed scheme. That is, the proposed scheme can effi-

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

11

Fig. 15. Deformed configurations of RVEs obtained by imposing PBCs: (a) a hexagonal RVE; and (b) an RVE with real microstructures.

ciently cover much larger domains than that of the conventional voxel discretization. 4.3. Results and discussion in terms of equivalent elastic properties Imposing the displacement-controlled PBCs to the meshes created from the artificial and real images in the previous subsections, we computed the equivalent elastic properties of M55J/M18 composites. The displacement relationships, corresponding to each of the six macroscopic strains (<ε 11 >, <ε 22 >, <ε 33 >, <ε 23 >, <ε 13 >,

<ε 12 >) for 3D models of UD composites, were first prescribed at the counterpart nodes on each pair of RVE boundaries. The average stresses (<σ 11 >, <σ 22 >, <σ 33 >, <σ 23 >, <σ 13 >, <σ 12 >) over the RVE were then calculated in the FEM framework, which correspond to the macroscopic stresses of the UD composites. As a result, the equivalent elastic moduli can be extracted from the macroscopic stress-strain (<σ > − <ε>) relationship [33]. The fiber and matrix elastic properties are presented in Table 1. The fiber properties were obtained by inverse analysis with the Mori-Tanaka method

12

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

Fig. 16. Comparison of elastic properties obtained using various RVE models: (a) E11 ; (b) E22 ; (c) G12 ; (d) G23 ; and (e) ν 12 .

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14 Table 1 Elastic properties of the M55J fiber and M18 matrix. Elastic property E11 (GPa) E22 (GPa) G12 (GPa) G23 (GPa)

ν 12

Reference value M55J (Fiber)a

M18 (Matrix)b

496.52 6.38 17.92 2.78 0.25

4.2 4.2 1.5 1.5 0.4

a Predicted values using the scheme proposed by Rupnowski et al. [34] with M55J/M18 elastic properties [35]. b Acquired values from the Hexcel corporation [36].

[34,37], in which the known matrix and lamina properties and the fiber-volume fraction values were used. In the case of RVEs with the periodic fiber arrays as in Section 4.1, node-to-node periodic pairs between two opposite RVE faces/edges are naturally generated, and thus, it is straightforward and easy to apply the PBCs. In contrast, the RVEs that reproduce real microstructures as in Section 4.2 fail to make matching node pairs for all boundary elements of the RVEs because the real fiber arrays are not periodic at all. In this work, to enforce the periodic constraints on the RVEs containing nonperiodic microstructures, only the boundary elements with nonmatching node pairs were locally triangulated by adding nodes so that all matching nodeto-node pairs were created between two opposite RVE boundaries. Consequently, we can directly impose the periodic constraints for the RVEs that come from both the artificial and real images. For comparison purpose, we additionally prepared three simple RVEs with hexagonal, square and diamond arrays shown in Fig. 14. Through the computational homogenization with PBCs [33], we obtained the equivalent elastic properties of the M55J/M18 composites for the simple RVEs, as well as for the RVEs created from the artificial and real images. All RVEs have a similar fiber-volume fraction at 68.33%. The deformed configurations for the prescribed macroscopic strains with the PBCs are shown in Fig. 15, representatively, for the hexagonal RVE and the real microstructure RVE. Moreover, the obtained values of the equivalent elastic properties are plotted in Fig. 16. The results of the 100 RVEs containing the artificially created microstructures are shown in the form of histograms and their corresponding normal distribution curves with the mean (m) and standard deviation (s) values, and the other results are indicated by the point symbols. There is mostly good agreement among the results for the different RVEs, and we can thus conclude that the proposed scheme is a promising tool for efficient microstructure modeling. Moreover, from a relatively large variation of the elastic properties except for the longitudinal Young’s modulus E11 , it was found that the equivalent properties are affected by the configurations of the microstructures. Hence, it is important to represent the detailed microstructures in numerical models for accurate and realistic prediction of material behaviors. Meanwhile, the values of G12 obtained using the RVEs containing real microstructures show a deviation 15% higher than those obtained using the 100 RVEs with artificial microstructures. The deviation might be attributed to enforcing the PBCs, instead of the kinematic uniform boundary conditions, on the real microstructure RVEs which have no periodic configurations of microstructures, as reported in [38–40]. 5. Conclusion In this work, we proposed a simple and efficient microstructure modeling scheme using the SDF based boundary smoothing algorithm and element trimming technique, with emphasis on its application to UD composites. The performance of the proposed

13

scheme was verified by two benchmark problems. From the results, it was found that the proposed scheme provides solutions as accurate as those when conventional FE meshes are used based on CAD software. Furthermore, as a practical application to UD composites, we constructed FE models based on image processing. The proposed scheme can reproduce detailed configurations of artificial and real microstructures in an accurate and efficient manner. In particular, in the case of the RVEs based on real images, no assumption was made about the shape, size, and distance of the fibers. Although in this work, we limited our attention to the evaluation of the elastic properties of UD composites, the proposed scheme can be applied to the evaluation of other multiphase solids or other mechanical properties, such as compressive strength, residual thermal stress, and so forth. The scheme can be also extended to 3D modeling of randomly oriented microstructures in RVEs. The testing of further applications are currently in progress and will be reported when completed. Acknowledgment This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A1A2058549 continued with 2017R1A1A2A10 0 0 0716). This was also supported by research funds for newly appointed professors of Chonbuk National University in 2016. References [1] Moore JA, Ma R, Domel AG, Liu WK. An efficient multiscale model of damping properties for filled elastomers with complex microstructures. Compos Part B Eng 2014;62:262–70. [2] McVeigh C, Liu WK. Linking microstructure and properties through a predictive multiresolution continuum. Comput Methods Appl Mech Eng 2008;197(41):3268–90. [3] Melro A, Camanho P, Pinho S. Generation of random distribution of fibres in long-fibre reinforced composites. Compos Sci Technol 2008;68(9):2092–102. [4] Wongsto A, Li S. Micromechanical finite element analysis of unidirectional fibre-reinforced composites with fibres distributed at random over the transverse cross-section. Compos Part A Appl Sci Manuf 2005;36(9):1246–66. [5] Trias D, Costa J, Mayugo J, Hurtado J. Random models versus periodic models for fibre reinforced composites. Comput Mater Sci 2006;38(2):316–24. [6] Terada K, Miura T, Kikuchi N. Digital image-based modeling applied to the homogenization analysis of composite materials. Comput Mech 1997;20(4):331–46. [7] Reid AC, Langer SA, Lua RC, Coffman VR, Haan S-I, García RE. Image-based finite element mesh construction for material microstructures. Comput Mater Sci 2008;43(4):989–99. [8] Sharma R, Mahajan P, Mittal RK. Elastic modulus of 3D carbon/carbon composite using image-based finite element simulations and experiments. Compos Struct 2013;98:69–78. [9] Tarleton E, Charalambides M, Leppard C. Image-based modelling of binary composites. Comput Mater Sci 2012;64:183–6. [10] Lian W-D, Legrain G, Cartraud P. Image-based computational homogenization and localization: Comparison between X-FEM/levelset and voxel-based approaches. Comput Mech 2013;51(3):279–93. [11] Legrain G, Cartraud P, Perreard I, Moës N. An X-FEM and level set computational approach for image-based modelling: application to homogenization. Int J Numer Methods Eng 2011;86(7):915–34. [12] El Moumen A, Imad A, Kanit T, Hilali E, El Minor H. A multiscale approach and microstructure design of the elastic composite behavior reinforced with natural particles. Compos Part B Eng 2014;66:247–54. [13] Kim JH, Lee M-G, Wagoner R. A boundary smoothing algorithm for image-based modeling and its application to micromechanical analysis of multi-phase materials. Comput Mater Sci 2010;47(3):785–95. [14] Duarte A, Silva B, Silvestre N, de Brito J, Júlio E. Mechanical characterization of rubberized concrete using an image-processing/XFEM coupled procedure. Compos Part B Eng 2015;78:214–26. [15] Scan2CAD (accessed 17.01.17). [16] Img2CAD (accessed 17.01.17). [17] Fries TP, Belytschko T. The intrinsic XFEM: a method for arbitrary discontinuities without additional unknowns. Int J Numer Methods Eng 2006;68(13):1358–85. [18] Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite-element method. Comput Methods Appl Mech Eng 2001;190(46):6183–200.

14

J.H. Lim et al. / Advances in Engineering Software 109 (2017) 1–14

[19] Wu C, Guo Y, Askari E. Numerical modeling of composite solids using an immersed meshfree Galerkin method. Compos Part B Eng 2013;45(1):1397–413. [20] Tian R, To AC, Liu WK. Conforming local meshfree method. Int J Numer Methods Eng 2011;86(3):335–57. [21] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling. Model Simul Mater Sci Eng 20 09;17(4):0430 01. [22] Sohn D, Han J, Cho Y-S, Im S. A finite element scheme with the aid of a new carving technique combined with smoothed integration. Comput Methods Appl Mech Eng 2013;254:42–60. [23] Sethian JA. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge University Press; 1999. [24] Matlab User’s Guide (R2014a), MathWorks, 2014. [25] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. Springer Science & Business Media; 2006. [26] B. Sumengen, A matlab toolbox implementing level set methods, http:// barissumengen.com/level_set_methods (accessed 01.05.16). [27] Enright D, Fedkiw R, Ferziger J, Mitchell I. A hybrid particle level set method for improved interface capturing. J Comput Phys 2002;183(1):83–116. [28] Canny J. A computational approach to edge detection. IEEE Trans Pattern Anal Mach Intell 1986;8(6):679–98. [29] Roth SD. Ray casting for modeling solids. Comput Graph Image Process 1982;18(2):109–44. [30] Abaqus Analysis User’s Manual (6.14), Dassault Systèmes, 2014.

[31] Swaminathan S, Ghosh S. Statistically equivalent representative volume elements for unidirectional composite microstructures: part II-with interfacial debonding. J Compos Mater 2006;40(7):605–21. [32] Axio-imagers 2 (accessed 17.01.17). [33] Lim JH, Henry M, Hwang D-S, Sohn D. Numerical prediction of fiber mechanical properties considering random microstructures using inverse analysis with quasi-analytical gradients. Appl Math Comput 2016;273:201–16. [34] Rupnowski P, Gentz M, Sutter J, Kumosa M. An evaluation of the elastic properties and thermal expansion coefficients of medium and high modulus graphite fibers. Compos Part A Appl Sci Manuf 2005;36(3):327–38. [35] Romeo G, Frulla G. Analytical and experimental results of the coefficient of thermal expansion of high-modulus graphite-epoxy materials. J Compos Mater 1995;29(6):751–65. [36] Hexcel Corporation (accessed 17.01.17). [37] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica 1973;21(5):571–4. [38] Zhu H, Hobdell J, Windle A. Effects of cell irregularity on the elastic properties of 2D Voronoi honeycombs. J Mech Phys Solids 2001;49(4):857–70. [39] Silva MJ, Hayes WC, Gibson LJ. The effects of non-periodic microstructure on the elastic properties of two-dimensional cellular solids. Int J Mech Sci 1995;37(11):1161–77. [40] Silva MJ, Gibson LJ. The effects of non-periodic microstructure and defects on the compressive strength of two-dimensional cellular solids. Int J Mech Sci 1997;39(5):549–63.