An Atomistic-Based Continuum Modeling for Evaluation of Effective Elastic Properties of Single-Walled Carbon Nanotubes

The mechanical behavior of SWCNTs is characterized using an atomistic-based continuum method. At nanoscale, interatomic energy among carbon atoms and the corresponding force constants are defined. Subsequently, we used an atomistic finite element analysis to calculate the energy stored in the SWCNT model, which forms a basis for calculating effective elastic moduli. In the finite element model, the force interaction among carbon atoms in a SWCNT is modeled using load-carrying structural beams. At macroscale, the SWCNT is taken as cylindrical continuum solid with transversely isotropic mechanical properties. Equivalence of energies of both models establishes a framework to calculate effective elastic moduli of armchair and zigzag nanotubes. This is achieved by solving five boundary value problems under distinct essential-controlled boundary conditions, which generates a prescribed uniform strain field in both models. Elastic constants are extracted from the calculated elastic moduli. While results of Young’s modulus obtained in this study generally concur with the published theoretical and numerical predictions, values of Poisson’s ratio are on the high side.


Introduction
Extensive research work done by researchers from science and engineering background in composite materials opens new prospects for future short and long term technologies, which will reshape the practical application of modern composites.Currently, the research themes on nanocomposites and/or composites with nanoreinforcements face the challenges of characterization, fabrication, and application.Significant amount of experimental and numerical research work is done to characterize the nanoreinforcement.But further research is needed to bring these to the level of practical application.These nanocomposites are becoming favorable candidates for materials with a bright future in a wide variety of industries such as transport, defense, electronics, and biomedicine, to name a few.Hence, it is important that the mechanical properties of these composite constituents, particularly the carbon nanotubes (CNTs), be predicted accurately.Further, the potential use of carbon nanotubes (CNTs) as a reinforcing material in nanocomposites and light weight composite structures has triggered a need to explore their mechanical properties and assess their deformation under mechanical loading.The unique structure and geometric configuration of CNTs along with their high stiffness, low density, and large aspect ratio have propelled an increasing demand in furthering the research to quantify their elastic properties as well as to explore possible applications in different fields.
Various experimental and theoretical approaches have been developed or used to characterize the elastic behavior of SWCNTs.Several investigators [1,2] have conducted experimental studies to investigate the mechanical properties of carbon nanotubes.These experiments were mainly based on atomic force microscopy (AFM) and transmission electron microscopy (TEM) and were able to confirm that CNTs possess superior mechanical properties.However, the experimental error bars are too large to state exact characteristics of CNTs of different configurations, sizes, and structures.The wide scatter in the experimentally reported values of the elastic constants of the CNT can be attributed to the lack of proper direct measuring techniques at nanometer scale, difficulties in test specimen preparation, and the dissimilarities in the method of nanotubes manufacture [3].Such high complexity in the experimental characterization has prompted many researches to pursue a variety of theoretical studies on determining the effective mechanical properties of nanotubes.
The theoretical approaches found in literature can be divided into three main categories: the atomistic methods, the continuum mechanics modeling, and the equivalent continuum modeling using finite element method.The atomistic approaches include classical molecular dynamics (MD) [4], tight-binding molecular dynamics, and density functional theory (DFT) [5].Phenomenological interatomic potential energy functions are used in these approaches to model the nanoscale systems in order to determine the force applied by the carbon atoms.Therefore, the more realistic and accurate these potentials are, the more closely the results match the experimental data and the better they reflect the actual properties of real system.Although these atomistic methods can simulate any problem associated with molecular or atomic motions, their huge computational tasks bound their application to problems with small number of molecules or atoms.
The continuum mechanics based approaches employ theories of shells, trusses, and beams [6] to model CNTs.One advantage of the continuum shell modeling is that it can efficiently calculate both static and dynamics properties of CNTs.However, traditional continuum models cannot accurately describe the mechanical properties of CNT structures because they lack the atomistic representation and appropriate constitutive relations that govern material behavior at nanoscale.Hence, there is a need for the development of new modeling techniques to accurately capture the mechanical behavior of CNTs.Thus, nanomechanics continuum theories that integrate continuum mechanics theories with interatomic potentials of atomic and molecular structure were developed.Equivalent continuum modeling (ECM) approach is one of the major developments of continuum method.It has been regarded as a very efficient method, especially with nanostructures modeled at large scale.Molecular mechanics combined with finite element method (FEM) involving shell, beam, spring, rod, and combination of these models form the framework of ECM approach.Over the past years, many ECM models were presented in literature.Atomistic-based continuum multiscale modeling techniques were used to predict the mechanical behavior of CNTs considering the interatomic interactions at nanoscale [7].Li and Guo [8] proposed an equivalent continuum beam model that is capable of modeling interatomic forces between carbon atoms to compute effective elastic constant of CNTs.The elastic constants of beam elements in a finite element model were determined using a linkage between molecular and structural mechanics.An improved beam element, which includes the bond inversion energy, is proposed by Lu and Hu [9] to evaluate mechanical properties of graphene and SWCNTs based on molecular mechanics.A finite element approach based on molecular mechanics was proposed by Sun and Zhao [10].The chemical bond was simulated with a two-node elastic rod element and an elastic joint at each end.A tensile modulus of 0.4 TPa was found.Shokrieh and Rafiee [11] investigated Young's modulus of CNTs based on a nanoscale continuum modeling by employing frame elements to simulate C-C bonds.Another continuum model that allows calculation of Young's and shears moduli based on structural mechanics combined with FEA was developed by Muc [12].Young's moduli are derived from the natural frequencies of CNT structures.These models assume that the material is transversally isotropic.
Moreover, different kinds of atomistic finite elements including rods, trusses, beams, and springs have been used to model carbon-carbon (C-C) link in CNTs [6].Giannopoulos et al. [13] constructed a computational FE model to simulate the SWCNT using linear interatomic potentials for C-C bonds.Meo and Rossi [14] developed a finite model including both nonlinear elastic and linear torsional spring elements to represent the modified Morse potential when simulating SWCNTs.A full nonlinear finite element model was developed using spring elements accounting for both C-C bond stretching and C-C-C bond angle variations to investigate the effect of chirality and the diameter on Young's modulus of SWCNT [15].The elastic constants, Young's and shear modulus, of SWCNTs were numerically computed via finite element method incorporating Poisson's effect in the estimation of Young's and shear modulus of SWCNT [16].The numerically calculated values of Young's and shear modulus were approximately 1.046 TPa and 0.424 TPa, respectively, for a SWCNT having a thickness of 0.34 nm.
In spite of the variety of theoretical studies, there still remain differences in approaches and calculated values regarding the effect of geometric structure of CNTs on elastic constants, as evident by the wide scatter among elastic constants reported in the literature.The objective of this paper is to study the elastic behavior of single-walled carbon nanotubes (SWCNTs) using a multiscale modeling approach.At nanoscale, interatomic interactions among carbon atoms are modeled by a structural beam in atomistic finite element approach in conjunction with molecular structural mechanics.At macroscale, an equivalent continuum modeling method is proposed to compute directly the effective elastic moduli of SWCNTs having different chirality and configuration.The elastic constants, that is, Young's modulus and Poisson's ratio, are calculated from the elastic moduli.

Atomic Structure of SWCNTs
Several methods to construct CNTs are reported in literature.In a common procedure, CNT is built by rolling a graphene sheet in two directions: a specific rolling direction and the circumference of the tube cross section, as shown in Figure 1.In general, two chiral indexes (, ) are used to define the configuration of carbon nanotubes, where () normally is greater than ().The nanotubes with (, ) are typically labeled as armchair, while the structure with (, 0) is usually labeled as zigzag.The translation vector, T, is parallel to the tube axis and perpendicular to the tube chiral vector C h .The unit vectors of graphene sheet lie along two sequent "zigzag" lines and are represented by a 1 and a 2 .The vector a 2 has a different magnitude than a 1 and when added together, they equal the chiral vector C h .A complete graphical description of the variables and the relations governing the geometry of SWCNT can be found in [17].The following formulae are used to define the diameter,  CNT , and the coordinate of 3-D SWCNT:

T Tube axis
where  c-c is the bond length between carbon atoms, , ,  are the generalized coordinates of an arbitrary point in nanotube,  is the radius of nanotube (i.e.,  =  CNT /2), and ,  are the coordinates of graphene sheet.The properties of CNTs are dependent on its structural configurations, the number of the concentric layers, and their structural defect and impurities.The present paper focuses on characterizing the effect of geometric configuration on mechanical properties of SWCNTs.

Molecular Mechanics Theory
Molecular mechanics forms the basis for the developed method to characterize the effective mechanical properties of SWCNTs at macroscale.In a molecular structure of CNT, the nuclei of carbon atoms in hexagonal arrangement are treated as material points.The mechanics of the nuclei motion are controlled via the force field generated by electron-nucleus interactions and nucleus-nucleus interactions.This force field is derived from the steric potential energy, which is a function of the relative position of the nuclei in the molecule.However, due to a very large difference in mass between electron and nucleus of the carbon atom, the electronic-nucleus interactions are neglected in atomistic analysis according to Born-Oppenheimer approximation.Moreover, all thermal, electromagnetic, and quantum effects are also neglected.The interatomic potential depends on the relative positions of the carbon nuclei in a deformed state and is defined as the sum of all energies due to the interatomic interactions [18].It can be expressed as where   is energy due to bond stretching,   is energy due to the bond angle variation,   is energy due to dihedral torsion, and   is energy for an improper torsion (out-ofplane).The first four terms of (2) are known as bonded energies.The last two terms are associated with van der Waals and electrostatic interactions, respectively, and are known as nonbond energies.Generally, for a covalent system, bonded interactions have a dominant contribution in the steric potential energy, while the effect of nonbonded terms are negligible, which is also the case considered in current study.Assuming small deformation of SWCNT, the motion of material points can be approximated as a simple harmonic motion.This is used to identify the component of energy needed to define total steric potential energy in order to develop atomistic-based continuum mechanics model for determining elastic moduli of SWCNT.Accordingly, the bond stretching and bending energies are defined as while the sum of the torsion related energies,   , is defined as where   ,   , and   are the bond stretching, bond angle variation, and bond-torsional force constants, respectively, and Δ  , Δ  , and Δ  represent the change in the bond length, bond angle, and twisting bond angle, respectively.It is important to note that force constants are a direct measure of the resistant forces that are obtained either from MD simulation or from atomistic computation [18].The total atomistic energy in a SWCNT is found by summing the total energies of individual bonds among all carbon atoms in the deformed state and expressed as where superscript "" refers to th bond between C-C atoms and   is the total number of bonds in the SWCNT.Each bond is represented by a structural beam element as explained in the following section.

Atomistic Finite Element Modeling of SWCNTs
In this study, each carbon atom is considered as a material point, which is modeled as a node in finite element model.Hence, all carbon atoms are represented as nodes and all C-C bonds are represented as beam elements, which form the basic building blocks, as shown in Figure 2. The classical structural 3D beam finite element, BEAM188, in ANSYS © , a commercial finite element software, is used to model the interatomic forces of the covalent bonds present between carbon atoms.This element is defined by two end nodes with each node having six degrees of freedom (i.e., 3 translations and 3 rotations).

Geometrical and Material Properties.
The coordinate of all nodes in the atomistic finite element model is calculated using (1b).The beam length,  c-c , is equal to the bond length between carbon atoms in the hexagonal lattice in their lowest energy state.The diameter of the beam is assumed to be equal to the wall thickness of the nanotube.However, different values of the wall thickness have been assumed in the literature in order to find out the nearest values of the elastic properties computed by other simulation methods.In the present study, the material properties of the beam with circular cross section are calculated using the force field constants that are obtained for molecular dynamic simulations, as will be explained later.A MATLAB © code is developed to generate different configurations of SWCNT that are used in this study.The generated script output files are then imported to ANSYS © to construct a working finite element model.
It is important to note that the required geometrical and material properties of the beam element are determined by equating the molecular potential energies, described by (5), to strain energies of equivalent structural beam in the proposed finite element model.No coupling is assumed among the distinct types of stored energies due to stretching, bending, and torsional deformations.As a result, equivalency of the corresponding terms in two systems can be established.For pure bending, the total strain energy in a uniform beam of length  c-c subjected to a uniaxial force  is given by where    is the strain energy of the th beam element due to axial stretch and , , and  c-c are Young's modulus, crosssectional area, and length of the beam.This energy represents the oscillation between two carbon atoms at atomistic scale.Comparing ( 6) with (3a) gives equating the strain energy stored in the th beam due to pure bending deformation, which is the energy generated from an angle change between two connected beams to atomistic level energy (see (3b)).This can be expressed as Following similar procedure for pure torsional deformation yields Algebraic manipulation of ( 7), (8), and ( 9) yields geometrical and material properties of beam in terms of resistive force constants as given below:

Loading and Boundary Conditions.
Once the finite element is developed with complete details of geometry and material, a set of boundary conditions are needed to obtain the solution.The solution in this case will depend on the type of continuum model considered at macroscale, which is discussed in Section 5. Due to the hexagonal symmetry on the cylindrical surface of the SWCNT, it is treated as transversely isotropic material, with the plane of isotropy perpendicular to the axial axis of the nanotube.A transversely isotropic material requires five independent material parameters to determine the entire set of elastic constants.Each one of these material parameters is determined from a single set of boundary conditions applied to SWCNT such that they yield unique value of independent elastic constants.This requires determination of atomistic energy in SWCNT under single boundary conditions relevant to elastic constant under consideration.
Five sets of boundary conditions considered in this study are summarized in Figure 10.The steps used for finite element analysis of SWCNT to calculate the total energy are summarized in Table 1.The finite element model is run separately for each set of boundary conditions.The SENE label on the ETABLE command, in ANSYS © , is used to Represent carbon atom as a node and C-C bond as beam element in finite element model.
(3) MATLAB5 software is used to generate the nodal coordinates with appropriate element connectivity for armchair and zigzag SWCNT finite element model using (1a) and (1b).( 4) Import the generated script files to ANSYS software to construct the finite element model.(5) Determine geometrical and material properties of beam finite element using (10) and then feed them to the FE model.(6) Apply loading and boundary conditions for each loading case described in Figure 10.For every simulation run: (6.1) Set the total energy to zero   = 0. (6.2) Solve the FE model for the selected loadings and boundary conditions.(6.3) Calculate the energy for every th bond (i.e.,    ).(6.4) Add this energy to the total energy for the selected loading case.(7) Use the total strain energy,  atomistic , to calculate the elastic constants,   , for each SWCNT configuration.

Equivalent Continuum Modeling of SWCNTs
As described in the previous section, the SWCNT is modeled at atomistic level to determine the total strain energy in the system subjected to prescribed loading/boundary conditions.
In this section we focus on developing an equivalent continuum model of SWCNT under the same prescribed loading/boundary conditions.We assume that, at macroscale, the SWCNT can be modeled as a hollow cylinder with unknown effective elastic material moduli   , where   depends upon type of material symmetry.The hollow cylinder is taken as an equivalent continuum model of SWCNT at macroscale having an inner diameter of " CNT ," thickness "," and length " CNT ," in cylindrical coordinate system (, , ), as shown in Figure 3.A similar approach was adopted to evaluate effective mechanical properties of graphene sheets [30].
Assuming transverse isotropic material behavior, linear elastic deformation, the generalized state of stress,   , in the equivalent continuum model under given generalized strain,   , can be characterized using the following constitutive law: where   and   are second-order tensors of the stress and strain, respectively, and   is the fourth-order elasticity tensor.Incorporating the transversely isotropic symmetry, the stress-strain relationship which involves five independent moduli is reduced to where  11 ,  33 ,  13 ,  23 , and  66 are five independent moduli of equivalent continuum model and   = 1/( 33 +  23 ).It can be also expressed in terms of the engineering constant as where The strain energy for the continuum is given by Substituting (13a) and (13b) in ( 16) and integrating over volume give where  continuum is the strain energy and   is the initial volume of the equivalent continuum nanotube shown in Figure 3.The calculation of elastic moduli ( 11 ,  33 ,  13 ,  23 , and  66 ) requires solving five boundary value problems under a prescribed uniform (constant) infinitesimal strain filed inside the SWCNT, which corresponds to predefined essential-controlled boundary conditions on SWCNT boundaries according to in which    is a prescribed uniform strain filed in the interior domain,   is the corresponding displacement on boundaries, and   is the coordinate of the boundaries.
Accordingly, five different distinct sets of loading cases are required to solve the boundary value problems and each load case is used to determine a corresponding elastic modulus at a time.These loadings can be classified based on the type of the deformation to uniaxial, biaxial, and pure shear loading.Thus, each load case will lead to a special form of strain energy given in (17).In the first run (loading case # 1),   = 0.001 is specified and all other strains are assigned a zero value.This translates a uniaxial deformation in the radial direction prescribed on the tube boundaries, as shown in Figure 10, while both tube ends are restricted from deforming in the axial direction.According to (17), the stored energy in the continuum model is reduced to from which  11 is calculated directly.Similarly, in the second run (loading case # 2),   = 0.001 is specified and all other strains are assigned a zero value from which  33 is calculated.
In the third run (loading case # 3), a uniform biaxial displacement in the radial and circumferential direction is applied to the tube.The values of   and   are fixed to 0.001, while all other strains are assigned a zero value.The stored energy in the continuum model for this case is reduced to from which  13 is calculated directly after substituting obtained values of  11 and  33 .Likewise, in the fourth run (loading case # 4), a uniform biaxial displacement in the circumferential and axial direction is applied to the tube and all other strains are assigned to a zero value from which  23 is calculated.
In the fifth run (loading case # 5), a pure shear deformation is applied to the tube using   = 0.002.All other strains are set to zero.In other words, an angular displacement is prescribed on the tube circumference at one end, as shown in Figure 10, while the other end is restricted from deforming in the axial direction.The stored energy is then reduced to from which  66 is calculated directly.In summary the energy of atomistic structure, for every run, is obtained by the atomistic finite element analysis, as described in the previous section, and equated to the corresponding strain energy of the continuum solid as given by the special forms.Once the five components of the elastic moduli are determined, the elastic constants of the SWCNT can be calculated using (15).

Results and Discussion
The aforementioned methodology is applied to characterize the elastic deformation of various configurations of armchair and zigzag single-walled carbon nanotubes with a diameter that varies from 0.8 nm to 3.1 nm and an aspect ratio ranging from 5 to 30. Figure 4 shows atomistic FE models of (10, 10) armchair SWCNTs with different aspect ratios.Essential boundary conditions as given in (18) are prescribed on exterior edges of the carbon nanotubes to determine the corresponding elastic modulus.Five sets of loading cases are used to calculate the effective elastic moduli ( 11 ,  33 ,  13 ,  23 , and  66 ), from which all elastic constants (  ,   , ]  , ]  , and   ) are determined.Table 2 shows elastic moduli for both armchair and zigzag nanotubes having an aspect ratio equal to 20.The values of the axial shear modulus,   , are also derived from the elastic relation as is also shown in Table 2.It can be observed that, at low aspect ratio, the elastic properties are strongly dependent on size and chirality.Nevertheless, it asymptotically converges to a constant value for large diameter nanotubes.Also, elastic constants are slightly larger for armchair nanotubes compared to zigzag nanotubes, although the difference is minimal.
Figure 5 depicts the variation in radial Young's modulus versus nanotube diameter at different aspect ratios for both armchair and zigzag geometries, respectively.For the armchair nanotube, as the diameter increases radial Young's modulus converges to a value equal to 1.274 TPa.Also, at small diameters, a large variability in Young's modulus exists particularly at low aspect ratios (AR = 5-10).This variability vanishes at large diameters.The zigzag nanotubes follow similar behavior with an averaged converged value of radial Young's modulus equal to 1.252 TPa.The difference between radial Young's modulus for both geometries is approximately 2% where Young's modulus of armchair is slightly higher than the zigzag configuration.The figures also show that, for diameters greater than 2.5 nm, there is no effect of aspect ratio on   .Also, nanotubes with diameter less than 0.8 nm offer no advantage over conventional metallic materials.However between 0.8 nm and 2.5 nm diameter nanotube, one has to be carefully selecting right configuration of nanotube, diameter, and length, in order to exploit its strength potential.Different behaviors of both configurations (armchair and zigzag) in the axial direction exist, as shown in Figure 6.For armchair, the axial Young's modulus is higher at small diameters and decreases as diameter increases until it converges to a value of 1.261 TPa.The zigzag nanotube follows similar behavior with an asymptotic axial Young's modulus equal to 1.284 TPa.Contrary to the previous case, there is a small impact of aspect ratio on   for both geometries at high diameters.Similar argument of careful selection of diameter and length, as discussed before, is also true for this case (0.8 ≤  CNT ≤ 2.5 nm).
It is worthwhile to mention that the difference between calculated values of radial and axial Young's moduli in either case of armchair or zigzag nanotubes is negligible.This observation confirms the isotropic behavior of SWCNTs concerning stiffness as most researches have assumed with respect to both geometric configurations.
Figure 7 shows the variation of out-of-plane Poisson's ratio, ]  , with respect to carbon nanotube diameter and   aspect ratio for both armchair and zigzag geometries, respectively.As depicted in figures, armchair and zigzag nanotubes share a similar behavior; that is, ]  is higher at low diameters and decreases with increasing diameter till it converges to values of 0.362 and 0.339, respectively.The decrease in Poisson's ratio is steep for diameter greater than 0.8 nm and less than 2.5 nm for aspect ratio below 10, while the decrease is gradual as aspect ratio increases for same range of nanotube diameters.The in-plane Poisson's ratio, ]  , as shown in Figure 8 behaves differently as it has lower values at small diameters and increases with an increase in diameter till it converges to 0.352 and 0.323 at high aspect ratio for armchair and zigzag, respectively.One can notice a striking similarity between out-of-plane and in-plane Poisson's ratio (]  and ]  ) for both armchair and zigzag configurations.The two cases are a mirror image of each other at the whole range of aspect ratio showing a large variation in Poisson's ratio especially in the range between 0.8 nm and 2.0 nm, as shown in Figures 6 and 7.
Figure 9 shows the variation of shear modulus,   , with respect to carbon nanotube diameter and aspect ratio for both armchair and zigzag geometries, respectively.The figures show that both armchair and zigzag nanotubes behave similarly.Values of   are slightly higher at low diameters and decrease till converging to 0.442 TPa and 0.467 TPa for armchair and zigzag carbon nanotube, respectively.Also, the radial shear modulus, determined by loading # 5, matches the calculated axial shear modulus,   , shown in Table 2.This observation confirms the isotropic behavior of SWCNTs.
In order to examine the capability of current modeling method to predict the elastic behavior of SWCNTs with a reasonable accuracy, Table 3   theoretical predictions reported in the literature.As can be observed, Young's modulus and shear modulus are in a good agreement with the published results.On the other hand, the calculated values of Poisson's ratio are on the high side when compared to those calculated by other researchers which show a wide range of variation.
versus radial) of SWCNT.Instead, isotropic behavior was confirmed based on output results.The presented results demonstrate that the current multiscale-based model provides a valuable tool for prediction of elastic moduli of SWC-NTs.Therefore, it will be extended in future to modeling of thermal and mechanical behavior of MWCNT and SWCNT based nanocomposites.

Figure 1 :
Figure 1: Chiral vector defining the configurations and geometric parameters of carbon nanotubes [3].

Table 1 :
Steps involved in the development of atomistic finite element of SWCNT.Construct an atomistic finite element model of a SWCNT by rolling the graphene sheet of size (, ) considering all required parameters (e.g., physical and geometrical parameters).
The values of the elastic constants are derived by inversion of   , which can be deduced to  =  11  23 +  11  33 − 2 2 ,   are Young's moduli in the radial and axial directions of the carbon nanotube, respectively, ]  and ]  are Poisson's ratio in - and - plane, respectively, and   is the transverse shear modulus.The axial shear modulus,   , is derived from the elastic relation   =   /2(1 + ]  ).
illustrates a comparison of current values of Young's modulus, shear modulus, and Poisson's ratio with several corresponding experimental and this study, an atomistic-based continuum modeling approach is developed to calculate the effective elastic moduli of SWCNTs based on energy equivalence between atomistic model and continuum model.At low aspect ratio, Young's modulus values, as extracted from elastic moduli, show strong dependence on diameter for both armchair and zigzag configurations.Similar dependence of Poisson's ratio on geometric parameters of SWCNTs is confirmed.This effect diminishes as aspect ratio increases where both elastic constants become almost constant at high aspect ratios.Unlike other works reported in the literature, current modeling is capable of calculating elastic moduli ( 11 ,  33 ,  13 ,  23 , and  66 ) directly, from which elastic constants are extracted.Results of Young's modulus obtained in this study are in good agreement with several published theoretical and experimental predictions; however, values of Poisson's ratio are higher than those obtained by other researchers.A big discrepancy among estimated values of Poisson's ratio as reported in literature exists.Thus, further investigation is still needed.Furthermore, contrary to several previous studies, current modeling approach avoids assuming isotropic behavior (axial