Fracture model for the prediction of the electrical percolation threshold in CNTs/Polymer composites

Abstract

In this paper, we propose a 3D stochastic model to predict the percolation threshold and the effective electric conductivity of CNTs/Polymer composites. We consider the tunneling effect in our model so that the unrealistic interpenetration can be avoided in the identification of the conductive paths between the CNTs inside the polymer. The results are shown to be in good agreement with reported experimental data.

Keywords

electrical percolation / CNTs/Polymer composites / fracture model / electric conductivity / tunnelling effects

Introduction

Nowadays, composite materials are widely used in the aerospace industries due to their excellent mechanical properties and light weight. As the idea of structure-function combining design has become more and more popular, it is always preferred that the composite materials should have certain functional properties, such as anti-corrosion, heat insulation, electromagnetic shielding, electric conduction and so on. During flight, lightning strikes are always the big threats to the safety of airplanes. Due to the poor electric conductivity of traditional composite materials, a large amount of charge and heat will be accumulated at the fuselage of the plane, which can lead to severe structure damage. Therefore, it has been recognized that for flight safety composite materials used in aircarft should possess considerable electric conduction capability. In the past, carbon black and copper wires were added into polymer matrix to produce conductive composite materials. However, a relatively high concentration of these fillers are required in order to make practical electrical conductivities higher than 10⁻⁶ S/m [¹]. On the other hand, such high filler concentration will reduce the mechanical properties of the material due to the agglomeration fillers, which requires very strict and well controlled processing of materials.

Carbon Nanotubes (CNTs) are ideal candidates as fillers in composites due to their remarkable mechanical, thermal and electrical properties. At room temperature, intact CNTs can have a maximum electrical conductivity of 10⁸ S/m [²], which is even higher than the electric conductivity of copper. In addition, CNTs usually have very large length-to-diameter aspect ratio, ranging from several hundreds to thousands. A significant reduction in the percolation threshold by increasing aspect ratios of the filler particle has been confirmed by researchers [¹,³,⁴]. Hence, the high surface volume ratio of CNTs is very advantageous in making a percolation network. At the nano-scale, very high aspect ratio and high conductivity of CNTs make it possible to produce conductive composites at very low filling concentration and with much lower inhomogeneity as compared with larger particles used as fillers [⁵].

When sufficient numbers of fillers are added, a three-dimensional CNTs conductive network should be formed, and the original insulative composite material begin to behave like a conductor. The critical concentration of CNTs that characterizes a drastic increase in electrical conductivity is commonly termed as the electrical percolation threshold (EPT). This behavior is attributed to the formation of multiple, continuous electron paths, or conductive paths, in the polymer composites with the filler content at or above the percolation threshold. Over the past decades, great efforts by researchers have been made to explore how the percolation threshold is influenced by different factors such as the shape, size, aspect ratio, orientation, synthesis method and dispersion method of the inclusions. The electrical conductivity of single-walled carbon nanotubes (SWNTs) [⁶,⁷] and multi-walled carbon nanotubes (MWNTs) [³] were experimentally investigated. Some experimental literature reports ultralow value of volume fraction [¹,³,⁵,⁸], while others have used some different measures ranging from less than 1 to over 15% volume fraction [⁴,⁶]. A table of comprehensive collection of published experiments data of percolation threshold (using weight fraction wt%) from different researchers is outlined by Bauhofer et al [⁴]. When the original data are given in vol%, Bauhofer uses the conversion relation vol%=wt% for single wall nanotubes and vol%=2 wt% for multi-wall nanotubes regardless of the polymer matrix. As shown in their tables, for epoxy matrix and SWCNT, the highest measure is 1.0 wt% while the lowest is 0.005 wt%. For Epoxy matrix and MWCNT, the highest is 5.0 wt% while the lowest becomes 0.0021 wt%. It is also noticed that for each combination of different matrix and CNT, not more than two papers will give a report of percolation threshold higher than 0.2 wt%. So they provide a very reasonable speculation that with optimized dispersion methods a percolation threshold of 0.1 wt% can be achieved by nearly all CNT/Polymer systems. From the large variation of the percolation threshold and electrical conductivity in that table, we can see that percolation threshold phenomenon seems to be quite stochastic. This gives us insight to use stochastic method such as Monte-Carlo method to study it. Monte-Carlo method is one type of stochastic and undetermined simulation method first proposed by Stanislaw Ulam and John von Neumann. Its idea is very simple, namely using a large number of repeated numerical experiments to test the built-in properties of a system. But it is very powerful and can even be used to study determined problem.

To better understand and predict the overall response of nanocomposites, a lot of computational methods have been exploited. Griebel et al verify the rules-of-mixture by performing Molecular Dynamics Simulation (MD) of CNT/Polymer composites with infinite and finite CNT inside [⁹–¹¹]. Arash et al develop a Coarse-Grained (CG) model to predict the mechanical properties of CNTs/Polymer composites with high accuracy and efficiency [¹²]. MD has incomparable advantages in exploring material properties at nano scale and provides a lot of detailed and important information which is useful to help us understand some specific material properties and phenomenons occur at nano and micro scale. However, due to limited computational resources, MD alone is not sufficient and feasible for solving real engineering problems at macro scale. Homogenization techniques based on continuum mechanics and Representative Volume Element (RVE) [¹³] bridge the gap between nano and macro scale and help to extract effective properties from fine scale to coarse scale. Mortazavi et al study the effective thermal conductivity and elastic modulus of two phases composites with different cylindrical, spherical and platelet inclusions by Finite Element Method, Mori-Tanaka and Strong-Contrast Method [¹⁴,¹⁵]. While most contributions are focused on deterministic modeling of composites, there are also a few papers accounting for uncertainties [¹⁷–²⁰] including their propagation over several length scales [²²–²³], see also the software framework [²⁵,²⁶] for uncertainty analysis. The contributions [²⁷,²⁸] are among the contributions aiming to predict fracture related parameters of composites.

In applications in Aerespace Engineering, the macroscopic properties such as the electric conductivity of the composite is of interest. However, the macroscopic properties are governed by features at the finer length scales. For instance, the growth and interaction of micro-cracks during operation might influence local deformation, interrupt the conductive paths and change the macroscopic electric conductivity and percolation threshold of the composite. In order to account for such multi-scale phenomena, multiscale methods have been developed. They can be categorized into hierarchical, semi-concurrent and concurrent methods. In hierarchical multiscale methods, information is transferred from the fine-scale to the coarse-scale only. Semi-concurrent and concurrent multiscale approaches transfer also information back from the coarse-scale to the fine-scale. Typical concurrent multiscale methods include the Arlequin method [³⁰,³¹] which aimed to couple two continuum models through a handshake coupling, their extension to the bridging domain method [³²] coupling atomistic with continuum models, the bridging scale method [³³] or quasi-continuum approach [³⁴,³⁵], just to mention a few. They have been extended to fracture [³⁶–⁴⁴] exploiting the extended finite element method [⁴⁵,⁴⁶] which allows modeling crack growth without remeshing. Note that XFEM has been applied to other interesting applications including topology optimization [⁴⁸–⁴⁹] and inverse analysis [¹³⁰] to name a few. Other competetive methods for fracture include the phantom node method [⁵³], virtual node XFEM [⁵⁵], SXFEM [⁵⁶–⁵⁸], efficient remeshing techqniues [⁶¹–⁷⁰], meshfree methods [⁷¹–⁸⁶,¹²²–¹²⁴], dual-horizon peridynamics [⁸⁷], cracking particles methods [⁸⁸–⁹¹] and extended IGA [⁹³–¹⁰⁸] among others [¹¹¹,¹¹³].

Many semi-concurrent methods are based on computational homogenization, see e.g. the seminal work of [¹¹⁶]. Hence, they are not suitable to model fracture since the separation of length scales, one key assumption in computational homogenization, is violated. Hierarchical multiscale methods which are based on computational homogenization have the same problems. They are more efficient when the material response is linear but get increasingly demanding for non-linear material responses.

In this paper, we propose a 3D stochastic model to predict the percolation threshold and effective electric conductivity in CNT based polymer matrix composites with randomly dispersed CNTs. Besides, we also account for tunneling effect to explore how the percolation threshold is affected by the length of the CNTs and the tunneling effect domain width. This model can be very promising in future applications due to the possibility of using CNT/Polymer composites as strain and damage sensors for nondestructive testing. This is quite beneficial for structural health monitoring to identify the current state of a system and detect damage and degradation that could ultimately lead to a system failure in real time. In our model, percolation threshold is defined as the critical vol% of CNTS at or above which a conductive network will be formed with one hundred percent possibility. It is worth noting that this value does not mean that a conductive network have no chance to be formed below the critical value. Because of the stochastic nature of percolation threshold phenomenon, a conductive network can still be produced in a lower concentration system with elaborate processing methods. But the smaller the volume fraction of CNTs inside matrix is, the harder for a conducting path to be formed and the more human intervention is needed. Our model is completely 3D and avoids using the overlapping or interpenetration criteria to detect conductive paths [¹¹⁷]. Instead, we consider a certain width of tunneling effect domains around the hard-core CNTs and check the conduction between CNTs by the interference of tunneling effect zone. In the next section, we will explain our model for predicting the percolation threshold in CNT based poylmer-matrix composites. Subsequently, we address computational homogenization before we present our results in Section 3. The manuscript finishes with conclusions and future research perspectives.

Electro-mechanical model for polymer-matrix composites

In our model, we consider the CNTs as straight hard-core particles without any intersections. Around the CNT, there exist a capsular region with a certain thickness W representing the distance that the tunneling effect can happen between two CNTs. Fig. 1 shows a Carbon Nanotube and its surrounding tunneling effect zone. If the tunneling effect zones of two neighbor CNTs have overlapping, we think they can form a conductive path for electricity (Fig. 2).

In our simulation, we assume the CNTs as straight solid cylinders with a radius of 2 nm. Their length follows a gaussian distribution, with a mean value of 100 nm, 120 nm and a standard deviation of 0 nm, 20 nm, 40 nm respectively. The RVEs we generate have size of 300 nm. Figs. 3 and 4 shows two samples with 400 and 6000 CNTs respectively. For each volume fraction value of CNTs, we take 1000 random samples and get the possibility by checking whether a percolation network can be formed or not in three principal direction in each sample. In details, we check whether a conductive path can form in all of three principal directions. If any of these three tests fail, we regard the sample as a failure test and it can not form a overall conductive path. If all the three tests passed, we count the sample as a success. We use these short CNTs because they require smaller size of Representative Volume Element (RVE), which makes it easier to generate a microstructure of specific volume fraction. We implement this model by developing a Fortran code to generate geometry information and write the data to a file, then we use Python scripts to read the microstructure data and make it visible in ABAQUS\CAE. Fig. 5 displays a conductive path connecting the two opposite surface in X direction. As can be seen from the figure, the conductive path is usually not a simple straight line. In most cases, it will be a very complicated curve path.

Fracture modeling based on XFEM

Cracks can be modelled by several techniques in computational methods. One way is to align the discretization to the cracks. In case of crack propagation, the mesh needs to be adjusted at least in the vicinity of the crack tip [⁶¹,⁶⁶–⁶⁹,¹¹³]. Meshfree methods such as the element-free Galerkin method [¹¹⁸–¹²¹] or partition of unity enriched meshfree methods [⁸⁸–⁹¹,¹²²–¹²⁴] also offer the opportunity to capture the crack topology. However, modeling the complex micro-structure remains a challenge in meshfree method due to their lack of geometry description. One of the most popular method to treat fracture is the extended finite element method [⁴⁵,⁴⁶] which allows for modeling crack propagation without remeshing. Though the goal of this research is focused on predicting electric properties for stationary cracks, XFEM still has a substantial advantage. Ensuring a good mesh quality for such a complex micro-structure including a large number of cracks is very difficult. Therefore, we decided to model the cracks through a step-enriched XFEM-type approach which is similar to the so-called phantom node method [¹²⁵], a special version of the extended finite element method. The phantom node method does not require a crack tip enrichment which complicates the implementation and integration. Recent approaches [⁶⁶] also allow the crack tip to be located within an element which allows the use of relatively coarse meshes. In the step-enriched XFEM, the crack tip is also avoided. There are a few electro-mechanical XFEM formulations available [¹²⁶–¹³⁰]. However, the approaches [¹²⁶–¹²⁹] employ a crack tip enrichment and deal only with a few cracks.

The crack is considered as impermeable. Hence the approximation of the electrical field is given by

(1){\phi }^h\left(\mathrm{X}\right)={\displaystyle \sum _{I\in N}{{\displaystyle N}}_{\mathrm{1}}}\left(\mathrm{X}\right)\phi I+{\displaystyle \sum _{I\in N_c}{{\displaystyle N}}_I^D}\left(\mathrm{X}\right){{\displaystyle b}}_I,

where N and N_c denote the set of all nodes and enriched nodes in the discretization, {{\displaystyle \phi }}_I and {{\displaystyle b}}_I are the standard and enriched degrees of freedom associated electrical potential, respectively; NI (X) denote the standard finite element shape functions while {{\displaystyle N}}_I^D\left(X\right) are the enriched or discontinuous shape functions which are computed by the product of the standard finite element shape function and an enrichment function. In order to obtain a discontinuity in the electrical potential, the shifted Heaviside function {{\displaystyle H}}_I\left(f\left(\mathrm{X}\right)\right)=H\left(f\left(\mathrm{X}\right)\right)-H\left(f\left({{\displaystyle \mathrm{X}}}_I\right)\right) is chosen with

(2)

where X_c is a point on a crack segment and n is the normal to the associated crack segment. Hence, the enriched shape function is given by

(3){{\displaystyle N}}_I^D\left(\mathrm{X}\right)=H\left(f\left(\mathrm{X}\right)\right){{\displaystyle N}}_I\left(\mathrm{X}\right).

We plan to generate random initial micro cracks in polymer and analyze its influence on the effective electric conductivity. However, we believe close micro cracks have minor effect on the effective electric conductivity since both void and polymer have almost zero electric conductivity. When this kind of conductive composites are at service, they will also be exposed to different mechanical loads, so we intend to add various external load to our RVE model. At this stage, these initial micro cracks may have significant influence on the overall electric conductivity because cracks will change the distribution of local stress and strain field. Due to some local deformation, the connectivity of the tunneling effect zone will be broken. Due to Joule Effect, electricity will produce heat inside the composite, which will affect the electrical properties of constituents and the mechanical field inside the composites. So in the long term, we plan to include thermal field and perform multifield analysis to solve this coupled problem.

Computational homogenization

The ability to make artificial materials with tailored behavior and properties is the key to the success of modern engineering. A general method to obtain macroscopically desired responses is to enhance a base material by the addition of microscopic and nanoscopic particles. Usually this kind of composite materials have highly complex microstructures. The macroscopic characteristics of the composites are the aggregate and integral response of the assemblage of different phases. However, it is not feasible to perform direct simulations of macroscopic engineering structures composed of all microscale details and heterogeneity because very fine spatial discretization mesh is needed to get local field distribution accurately and capture the effects of micro heterogeneities. This detailed model using FEM will lead to system of equations containing billions of unknowns, which is far exceed the memory and capability of most modern computers [¹³]. Due to these facts, the utilization of homogenized material models and aforementioned multiscale techniques is very important to solve emerging problems in all the areas of modern science and engineering.

With the dramatic increase of computational power, computational homogenization method are more and more widely used in the analysis of heterogeneous materials. They can overcome a lot of of inherent limitation of the the analytical homogenization theory, such as the incapability to handle complicated microstuctures and the interaction between inclusions. The basic strategy of computational homogenization is the simulation of behavior of a Representative Volume Element (RVE) of microstructure of material. The RVE should contain enough information about the statistical description of microstructure and heterogeneity, but still be very small compared with the size of structural components at macroscale and can be regarded as a material point. Fig. 6 [¹³] illustrates the relation between the size of structural component (L₁), the size of RVE (L₂) and the size of inclusions (L₃).

The central aim of computation homogenization is to connect the effective properties of composites with the properties and arrangement of its constituents. This is accomplished by computing the volume average of internal fields which is the solution to a series of Boundary Value Problem. For computing the effective electric conductivity, we use the following formula:

\overline{\sigma }={\displaystyle \frac{\overline{J}}{E}}={\displaystyle \frac{\mathrm{1}}{V}}\times {\displaystyle \frac{{\displaystyle {\int }_{V}J\mathrm{d}V}}{\nabla U}}.

Among this equation, U is the potential difference between two opposite surfaces of a RVE, E is the electric field intensity, J and \overline{J} are the local and average electric current density respectively, and \overline{\sigma } is the effective electric conductivity. We set constant potential difference between two opposite surfaces of a RVE sample, then we calculate the electric current density distribution among the RVE. Next we compute the volume average of the electric current density and get the corresponding effective electric conductivity in that direction by using the aforementioned formula. Since each RVE has three pairs of opposite surfaces or three perpendicular direction, we need to repeat the previous process three times for each RVE sample.

Since we use solid cylinder to represent hollow CNT with a certain thickness, we need to do some reduction of CNT’s electric conductivity due to the change of cross section. The equivalence is based on the assumption that the number of electrons that pass through the cross section per second should be the same since each CNT has a fixed ability to provide electrons. This assumption leads to a conclusion that the equivalent electric conductivity should be inversely proportional to the cross section area. Gupta et al. [¹³³] report that the wall thickness of single-walled CNT will stay 0.137 nm after the radius exceeds 1.5 nm. So we take the cross section of CNTs as a annulus with average radius of 2 nm and thickness of 0.137 nm, shown in Fig. 7. If we take the electric conductivity of pure CNTs as 10⁷ S/m, we can get a equivalent electric conductivity of 1.37×10⁶ S/m. For polymer, we assign a electric conductivity of 10⁻¹⁴ S/m. For Tunneling effect zone, since there is no report about its electric conductivity, we assume it to be the geometric average, which is 1.17 × 10⁻⁴ S/m, of the electric conductivity of CNT and Polymer.

Fig. 8 depicts the finite element discretization of a RVE model. This RVE model has an edge of 0.3 μm and contains 300 identical CNTs. Each CNT has a radius of 0.002 μm and a length of 0.1 μm. The width of tunneling effect zone is set to be equal to the radius, 0.002 μm. Fig. 9 shows the distribution of electric potential when voltage difference is applied in X, Y and Z direction. respectively. Table 1 shows the homogenized electric conductivity of this RVE model. However, the results indicate some anisotropy in different directions. In the following work, we will take more and bigger samples to give better prediction. We will also verify the isotropy of our model by statistical analysis.

Results

By percolation possibility we mean the chance of one hundred percent to form a continuous conductive path to connect two opposite surface of a RVE in three directions. Some researches use the volume fraction corresponding to 0.5 percolation possibility as the percolation threshold. But it is also possible to use the volume fraction that corresponds to a percolation possibility slightly large than zero, because that value should have some connections with the ultra low percolation thresholds reported. However, in this work, we use the definition we give in previous section.

Fig. 10 shows the curve of percolation possibility with respect to volume fraction for different tunneling effect zone width (W). Figs. 11 and 12 are the percolation possibility curves for different standard deviation and different average length of CNTs (or aspect ratio).

Conclusion

From the simulation, we can draw the following conclusions:

1) From Fig. 12, we can see that the average length of CNTs seems to have mirror effects on the percolation threshold value since they all appear around 2.8 vol%. However, the start point of climbing in each curve shifts left with longer average length. This means that with longer average length or larger aspect ratio, it is more possible to produce materials with low volume fraction of inclusions but still contain conductive paths accidentally. This also explains why for CNTS with such a high aspect ration around 1000, we can get a ultralow percolation threshold 0.005 wt% by accident.

2) From Fig. 11, we find that if we increase the standard deviation of the length of CNTs, the start point of climbing in each curve seems to be the same and close to 1.5 vol%, but the percolation threshold of each curve increases and shifts to the right. This shows that with a fixed average length of CNTs, if the uniformity of the length of CNTs is not good, it will bring some bad influence and impede us to make conductive composite materials with low percolation threshold.

3) From Fig. 10, we can see clearly that larger width of tunneling effect zone reduces the percolation threshold significantly. This is quite reasonable because if we imagine that the width tends to infinite, no mater how many CNTs the material contain and how they are dispersed and oriented, the CNTs can always form conductive paths easily. So the corresponding percolation threshold tends to zero. This is quite consistent with the overall left shift of the curves with respect to increasing width of tunneling effect zone. For width about three times of the radius of CNT, the percolation threshold drops to around 0.8 vol%.

In the future work, we will continue to improve our model by further study. First of all, due to lack of experimental data and limitation of experimental tools at nanoscale, we are not clear about how to choose the width of tunneling effect zone, namely how close two CNTs should be to make tunneling effect possible. Besides, the shape of real tunneling effect zone may differs significantly from the hollow capsule shape we propose. So in this paper, we just propose a simple model to make some reasonable approximations and get some useful conclusions which is consistent with some experimental reports about percolation threshold in CNTs/Polymer nanocomposites, but still a lot of further work needs to be done. In addition, as we mentioned previously, these CNTs/Polymer-matrix nanocomposites work as both mechanical and functional materials in industrial applications, so they will always be under some external load. The initial micro defects in the nanocomposites can develop under the applied load and the polymer matrix may also suffers ageing at the same time. All these factors can influence the overall mechanical properties and change local deformation state inside. As a chain effect, the original conductive paths may be cut by cracks or closed by relative pull-apart of CNTs. This will for sure reduce the electric conductivity of the material. So next we need to develop a mutiscale model to include both fracture and tunneling effect.

References

[1] Sandler J, Shaffer M S P, Prasse T, Bauhofer W, Schulte K, Windle A H. Development of a dispersion process for carbon nanotubes in an epoxy matrix and the resulting electrical properties. Polymer, 1999, 40(21): 5967–5971

[2] Michael G H, Tejas N. Radiofrequency interaction with conductive colloids: permittivity and electrical conductivity of single-wall carbon nanotubes in sallne. Bioelectromagnetics, 2010, 31(8): 582–588

[3] Martin C A, Sandler J K W, Shaffer M S P, Schwarz M K, Bauhofer W, Schulte K, Windle A H. Formation of percolating networks in multi-wall carbon-nanotube–epoxy composites. Composites Science and Technology, 2004, 64(15): 2309–2316

[4] Bauhofer W, Kovacs J Z. A review and analysis of electrical percolation in carbon nanotube polymer composites. Composites Science and Technology, 2009, 69(10): 1486–1498

[5] Bryning M B, Islam M F, Kikkawa J M, Yodh A G. Very low conductivity threshold in bulk isotropic single-walled carbon nanotube–epoxy composites. Advanced Materials, 2005, 17(9): 1186–1191

[6] Ounaies Z, Park C, Wise K E, Siochi E J, Harrison J S. Electrical properties of single wall carbon nanotube reinforced polyimide composites. Composites Science and Technology, 2003, 63(11): 1637–1646

[7] Kymakis E, Amaratunga G A J. Electrical properties of single-wall carbon nanotube-polymer composite films. Journal of Applied Physics, 2006, 99(8): 56

[8] Ramasubramaniam R, Chen J, Liu H. Homogeneous carbon nanotube/polymer composites for electrical applications. Applied Physics Letters, 2003, 83(14): 2928–2930

[9] Griebel M, Hamaekers J. Molecular dynamics simulations of the elastic moduli of polymer–carbon nanotube composites. Computer Methods in Applied Mechanics and Engineering, 2004, 193(17): 1773–1788

[10] Frankland S J V, Caglar A, Brenner D W, Griebel M. Molecular simulation of the influence of chemical cross-links on the shear strength of carbon nanotube-polymer interfaces. Journal of Physical Chemistry B, 2002, 106(12): 3046–3048

[11] Zhu R, Pan E, Roy A K. Molecular dynamics study of the stress–strain behavior of carbon-nanotube reinforced epon 862 composites. Materials Science and Engineering A, 2007, 447(1): 51–57

[12] Arash B, Park H S, Rabczuk T. Mechanical properties of carbon nanotube reinforced polymer nanocomposites: A coarse-grained model. Composites. Part B, Engineering, 2015, 80: 92–100

[13] Quayum M S, Zhuang X, Rabczuk T. Computational model generation and rve design of self-healing concrete. Journal of Contemporary Physics, 2015, 50(4): 383–396

[14] Mortazavi B, Baniassadi M, Bardon J, Ahzi S. Modeling of two-phase random composite materials by finite element, mori–tanaka and strong contrast methods. Composites. Part B, Engineering, 2013, 45(1): 1117–1125

[15] Mortazavi B, Bardon J, Ahzi S. Interphase effect on the elastic and thermal conductivity response of polymer nanocomposite materials: 3d finite element study. Computational Materials Science, 2013, 69: 100–106

[16] Hamdia K M, Msekh M A, Silani M, Vu-Bac N, Zhuang X, Nguyen-Thoi T, Rabczuk T. Uncertainty quantification of the fracture properties of polymeric nanocomposites based on phase field modeling. Composite Structures, 2015, 133: 1177–1190

[17] Silani M, Talebi H, Ziaei-Rad S, Kerfriden P, Bordas S P A, Rabczuk T. Stochastic modelling of clay/epoxy nanocomposites. Composite Structures, 2014, 118: 241–249

[18] Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T. Stochastic predictions of interfacial characteristic of polymeric nanocomposites (pncs). Composites. Part B, Engineering, 2014, 59: 80–95

[19] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimization of fiber distribution in fiber reinforced composite by using nurbs functions. Computational Materials Science, 2014, 83: 463–473

[20] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[21] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimal fiber content and distribution in fiber-reinforced solids using a reliability and nurbs based sequential optimization approach. Structural and Multidisciplinary Optimization, 2015, 51(1): 99–112

[22] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[23] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96(PB):520–535

[24] Msekh M A, Silani M, Jamshidian M, Areias P, Zhuang X, Zi G, He P, Rabczuk T. Predictions of j integral and tensile strength of clay/epoxy nanocomposites material using phase field model. Composites. Part B, Engineering, 2016, 93: 97–114

[25] Hamdia K M, Zhuang X, He P, Rabczuk T. Fracture toughness of polymeric particle nanocomposites: Evaluation of models performance using bayesian method. Composites Science and Technology, 2016, 126: 122–129

[26] Ben Dhia H. Multiscale mechanical problems: the arlequin method. Comptes Rendus de l’Academie des Sciences Series IIB Mechanics Physics Astronomy, 1998, 326(12): 899–904

[27] Dhia H B, Rateau G. The arlequin method as a flexible engineering design tool. International Journal for Numerical Methods in Engineering, 2005, 62(11): 1442–1462

[28] Xiao S P, Belytschko T. A bridging domain method for coupling continua with molecular dynamics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(17): 1645–1669

[29] Wagner G J, Liu W K. Coupling of atomistic and continuum simulations using a bridging scale decomposition. Journal of Computational Physics, 2003, 190(1): 249–274

[30] Tadmor E B, Ortiz M, Phillips R. Quasicontinuum analysis of defects in solids. Philosophical Magazine A, 1996, 73(6): 1529–1563

[31] Shenoy V B, Miller R, Tadmor E B, Rodney D, Phillips R, Ortiz M. An adaptive finite element approach to atomic-scale mechanicsthe quasicontinuum method. Journal of the Mechanics and Physics of Solids, 1999, 47(3): 611–642

[32] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[33] Silani M, Talebi H, Hamouda A M, Rabczuk T. Nonlocal damage modelling in clay/epoxy nanocomposites using a multiscale approach. Journal of Computational Science, 2016, 15:18-23

[34] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[35] Silani M, Ziaei-Rad S, Talebi H, Rabczuk T. A semi-concurrent multiscale approach for modeling damage in nanocomposites. Theoretical and Applied Fracture Mechanics, 2014, 74(1): 30–38

[36] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[37] Budarapu P R, Gracie R, Bordas S P A, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[38] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. Molecular dynamics/XFEM coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture. International Journal for Multiscale Computational Engineering, 2013, 11(6): 527–541

[39] Belytschko T, Moës N, Usui S, Parimi C. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 2001, 50(4): 993–1013

[40] Sukumar N, Moës N. B Moran, and T Belytschko. Extended finite element method for three-dimensional crack modelling. International Journal for Numerical Methods in Engineering, 2000, 48(11): 1549–1570

[41] Ghasemi H, Park H S, Rabczuk T. A level-set based iga formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[42] Nanthakumar S S, Valizadeh N, Park H S, Rabczuk T. Surface effects on shape and topology optimization of nanostructures. Computational Mechanics, 2015, 56(1): 97–112

[43] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92-93: 242–246

[44] Kumar S, Singh I V, Mishra B K, Rabczuk T. Modeling and simulation of kinked cracks by virtual node xfem. Computer Methods in Applied Mechanics and Engineering, 2015, 283: 1425–1466

[45] Chen L, Rabczuk T, Bordas S P A, Liu G R, Zeng K Y, Kerfriden P. Extended finite element method with edge-based strain smoothing (ESM-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209-212: 250–265

[46] Bordas S P A, Natarajan S, Kerfriden P, Augarde C E, Mahapatra D R, Rabczuk T, Pont S D. On the performance of strain smoothing for quadratic and enriched finite element approximations (xfem/gfem/pufem). International Journal for Numerical Methods in Engineering, 2011, 86(4-5): 637–666

[47] Areias P, Msekh M A, Rabczuk T. Damage and fracture algorithm using the screened poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[48] Rabizadeh E, Saboor Bagherzadeh A, Rabczuk T. Goal-oriented error estimation and adaptive mesh refinement in dynamic coupled thermoelasticity. Computers & Structures, 2016, 173: 187–211

[49] Areias P, Rabczuk T, Msekh M A. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 322–350

[50] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh refinement. Computational Mechanics, 2016, 58(6): 1003–1018

[51] Areias P, Reinoso J, Camanho P, Rabczuk T. A constitutive-based element-by-element crack propagation algorithm with local mesh refinement. Computational Mechanics, 2015, 56(2): 291–315

[52] Areias P, Rabczuk T, Camanho P P. Finite strain fracture of 2d problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72(1): 50–63

[53] Areias P, Dias-Da-Costa D, Sargado J M, Rabczuk T. Element-wise algorithm for modeling ductile fracture with the rousselier yield function. Computational Mechanics, 2013, 52(6): 1429–1443

[54] Areias P, Rabczuk T, Dias-da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[55] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[56] Nguyen-Xuan H, Liu G R, Bordas S, Natarajan S, Rabczuk T. An adaptive singular es-fem for mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 252–273

[57] Rabczuk T, Belytschko T, Xiao S P. Stable particle methods based on lagrangian kernels. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12-14): 1035–1063

[58] Rabczuk T, Areias P. A meshfree thin shell for arbitrary evolving cracks based on an extrinsic basis. Computer Modeling in Engineering & Sciences, 2006, 16(2): 115–130

[59] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367–382

[60] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[61] Rabczuk T, Areias P M A, Belytschko T. A meshfree thin shell method for non-linear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[62] Amiri F, Millán D, Arroyo M, Silani M, Rabczuk T. Fourth order phase-field model for local max-ENT approximants applied to crack propagation. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 254–275

[63] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[64] Amiri F, Anitescu C, Arroyo M, Bordas S P A, Rabczuk T. Xlme interpolants, a seamless bridge between xfem and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[65] Talebi H, Samaniego C, Samaniego E, Rabczuk T. On the numerical stability and mass-lumping schemes for explicit enriched meshfree methods. International Journal for Numerical Methods in Engineering, 2012, 89(8): 1009–1027

[66] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[67] Nguyen V P, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[68] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically nonlinear three-dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[69] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[70] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[71] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37-40): 2437–2455

[72] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29-30): 2777–2799

[73] Nguyen-Thanh N, Valizadeh N, Nguyen M N, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on kirchhoff-love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[74] Ghorashi S S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based xiga for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[75] Thai T Q, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[76] Chan C L, Anitescu C, Rabczuk T. Volumetric parametrization from a level set boundary representation with pht-splines. CAD Computer Aided Design, 2017, 82: 29–41

[77] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[78] Thai C H, Ferreira A J M, Bordas S P A, Rabczuk T, Nguyen-Xuan H. Isogeometric analysis of laminated composite and sandwich plates using a new inverse trigonometric shear deformation theory. European Journal of Mechanics. A, Solids, 2014, 43: 89–108

[79] Thai C H, Nguyen-Xuan H, Nguyen-Thanh N, Le T H, Nguyen-Thoi T, Rabczuk T. Static, free vibration, and buckling analysis of laminated composite reissner-mindlin plates using nurbs-based isogeometric approach. International Journal for Numerical Methods in Engineering, 2012, 91(6): 571–603

[80] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wchner R, Bletzinger K U, Bazilevs Y, Rabczuk T. Rotation free isogeometric thin shell analysis using pht-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47-48): 3410–3424

[81] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S P A, Rabczuk T. Isogeometric analysis using polynomial splines over hierarchical t-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21-22): 1892–1908

[82] Zhuang X, Huang R, Liang C, Rabczuk T. A coupled thermohydro-mechanical model of jointed hard rock for compressed air energy storage. Mathematical Problems in Engineering, 2014, 2014

[83] Areias P, Rabczuk T, Camanho P P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[84] Frédéric Feyel, Jean-Louis Chaboche. Fe 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials. Computer Methods in Applied Mechanics and Engineering, 2000, 183(3): 309–330

[85] Zeng X, Xu X, Prathamesh M. Shenai, Eugene Kovalev, Charles Baudot, Nripan Mathews, and Yang Zhao. Characteristics of the electrical percolation in carbon nanotubes/polymer nanocomposites. Journal of Physical Chemistry C, 2011, 115(44): 21685–21690

[86] Belytschko T, Yun Y L, Gu L. Element-free galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[87] Belytschko T, Lu Y Y, Gu L, Tabbara M. Element-free galerkin methods for static and dynamic fracture. International Journal of Solids and Structures, 1995, 32(17): 2547–2570

[88] Zhuang X, Augarde C E, Mathisen K M. Fracture modeling using meshless methods and level sets in 3d: framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[89] Zhuang X, Zhu H, Augarde C. An improved mesh-less shepard and least squares method possessing the delta property and requiring no singular weight function. Computational Mechanics, 2014, 53(2): 343–357

[90] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[91] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[92] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23): 1391–1411

[93] Song J H, Areias P, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[94] Enderlein M, Ricoeur A, Kuna M. Finite element techniques for dynamic crack analysis in piezoelectrics. International Journal of Fracture, 2005, 134(3-4): 191–208

[95] Kuna M. Finite element analyses of crack problems in piezoelectric structures. Computational Materials Science, 1998, 13(1): 67–80

[96] Shang F, Kuna M, Abendroth M. Meinhard Kuna, and Martin Abendroth. Finite element analyses of three-dimensional crack problems in piezoelectric structures. Engineering Fracture Mechanics, 2003, 70(2): 143–160

[97] Béchet E, Scherzer M, Kuna M. Application of the x-fem to the fracture of piezoelectric materials. International Journal for Numerical Methods in Engineering, 2009, 77(11): 1535–1565

[98] Nanthakumar S S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 1: 153–176

[99] Gupta S S, Bosco F G, Batra R C. Wall thickness and elastic moduli of single-walled carbon nanotubes from frequencies of axial, torsional and inextensional modes of vibration. Computational Materials Science, 2010, 47(4): 1049–1059

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg