## Abstract

Motivated by a model that qualitatively captured the response of vertically
aligned carbon nanotube (VACNT) pillars in uniaxial compression, we consider the
uniaxial tensile response of a class of compressible elastic-viscoplastic
solids. In Hutchens et al. [“Analysis of Uniaxial Compression of Vertically
Aligned Carbon Nanotubes,” J. Mech. Phys. Solids, **59**, pp. 2227–2237
(2011), Erratum **60**, 1753–1756 (2012)] an elastic viscoplastic
constitutive relation with plastic compressibility, plastic non-normality, and a
hardening-softening-hardening hardness function was used to model experimentally
obtained uniaxial compression data of cylindrical VACNT micropillars. Complex
deformation modes were found in uniaxial compression, which include a sequential
buckling-like collapse of the type seen in experiments. These complex
deformation modes led to the overall stress-strain signature of the pillar not
being of the same form as the input material hardness function. A fundamental
question that motivates exploring the deformation of this class of
materials—both experimentally and theoretically—is how to extract the intrinsic
material response from simple tests. In this study we explore the relation
between the input material response and the overall stress strain behavior in
uniaxial tension using the constitutive framework of Hutchens et al. A simple
one-dimensional analysis reveals the types of instability modes to be expected.
Dynamic, finite deformation finite element calculations are carried out to
explore the dependence of diffuse necking, localized necking, and propagating
band deformation modes on characteristics of the hardness function. Attention is
devoted to uncovering implications for obtaining intrinsic material properties
of complex hierarchical structures; for example, vertically aligned carbon
nanotubes (VACNTs), from uniaxial tension experiments.

## Introduction

Vertically aligned carbon nanotubes (VACNTs), also referred to as carbon nanotube (CNT) forests, turfs, brushes, and mats exhibit a variety of deformation behaviors. For example, under compression and flat-punch indentation, some VACNTs display essentially a fully elastic or viscoelastic mechanical response [1–8] and full recoverability under cyclic loading [4,5,9,10], while others deformed permanently [9,11–17]. In uniaxial compression, VACNTs have been observed to deform in a progressive buckling-like deformation mode [1,11–14,17–21]. A phenomenological elastic-viscoplastic constitutive relation, which incorporated plastic compressibility, plastic non-normality, and a hardening-softening-hardening hardness function [22,23] qualitatively reproduced the deformation mode and overall stress-strain behavior seen in the experiments in [13,24]. Additional calculations in [22] revealed a variety of deformation modes depending on material parameters and enabled the identification of the parameter regime that gives rise to progressive buckling type of deformation. The calculations also showed that gradients in properties that may arise from the influence of the VACNT growth process on density, tortuosity, and intertubular friction, can play a significant role in determining the overall response [13,24].

The existing experimental and computational reports on the deformation of the entangled arrays of carbon nanotubes have largely been limited to compressive loading—either compression or indentation. VACNTs have potential uses in a variety of applications, which include electromechanical devices, switches, and actuators [12,25,26], in which they are subjected to extensive thermomechanical stresses and strains. Hence, it is useful to develop a predictive framework for the mechanical behavior of VACNTs under a wide range of loadings, which requires the knowledge of their intrinsic material properties. A key finding in [22,23] was that the response in compression is a structural response, i.e., the observed overall compressive stress-strain relation does not directly reflect the intrinsic material properties. Extraction of material properties from the overall stress-strain response is complex and not unique. Hence, it is of interest to explore whether other deformation modes could provide a more direct link between the overall structural and material response and the inherent material properties, either by itself or in conjunction with the compressive response. A natural choice in this regard is uniaxial tension.

Here we do not attempt to quantitatively model the response of VACNTs in uniaxial tension. Rather, we focus on assessing the relation between material parameters characterizing the hardening-softening-hardening hardness function in the constitutive formulation of Refs. [22,23] and the qualitative nature of the predicted response in uniaxial tension. We first use a one-dimensional analysis to set the stage and then carry out finite deformation axisymmetric finite element calculations using the constitutive framework in [22,23] for compressible elastic-viscoplastic solids with plastic non-normality. We restrict consideration to the variation of three constitutive parameters and explore their influence on the deformation mode and the overall stress-strain response.

## 1D Model

Insight into the possible nonuniform deformation modes that can be expected in uniaxial tension for the type of constitutive relation used here can be gained from a simple one-dimensional model. There is a significant body of literature on using one-dimensional analyses to gain insight into tensile instabilities, e.g., [27–33]. The presentation here follows that in [34]. The possible nonuniform deformation modes considered are: (i) diffuse necking, (ii) localized deformation, and (iii) a propagating neck.

where $J=V/V0=HA/H0A0$ ($V$ and $V0$ are the current and reference volumes, respectively) and $F=H/H0$.

When phrased this way, the maximum load condition for a compressible solid has the same form as the Considére [27] condition for an incompressible solid. Within this one-dimensional context, attainment of a force maximum (or, equivalently, a maximum nominal stress) is regarded as a criterion for the onset of diffuse necking in a tensile specimen.

where $[\u2003]$ denotes a jump.

so that in this one-dimensional context the conditions for the onset of diffuse necking and of deformation localization coincide.

In three dimensions and for rate independent response, the localization condition involves a bifurcation under all around displacement boundary conditions. The one-dimensional analog of this can be regarded as $A\xb7=0$ throughout the deformation history. Then, with $A0$ uniform, Eq. (11) also implies $\sigma \xb7=0$.

and if $\nu \u2227<1/2$, $\sigma \xb7$ can be negative with positive $\tau \xb7$.

Another type of instability is a propagating deformation band, which arises in a variety of contexts and with a variety of underlying physical mechanisms, e.g., [22,31,35,36]. In the constitutive framework here, it can occur for a hardening-softening-hardening solid and the formulation given follows that in [31,34]. The conditions that must be satisfied across a neck front propagating with velocity $Wb$ are balance of linear momentum and compatibility.

where $u$ is the displacement, and $\rho $ is the density in the reference configuration.

where $\Phi e$ is the stored elastic energy and $\Phi d$ is the plastic dissipation. If $\tau =g$, $\Phi =\u222bgd\u025b$.

where $save$ is the average of the values of $s$ on each side of the traveling front and a common factor $Wb$ has been canceled from each term in Eq. (20).

where $sM$ is termed the Maxwell stress. Equation (21) expresses the “equal area” construction to determine the Maxwell stress. For a hardening-softening-hardening relation, $A$ and $B$ denote points on the two hardening branches and $\Phi A$ is the area under the maximum and $\Phi B$ is the area above the minimum. With material rate independence and $\tau =g$, $\Phi =\u222bg(\u025b)d\u025b$.

The cylindrical samples considered here are constrained by a substrate, so the initial deformation is not uniform. The material is rate dependent so the simple expressions for rate independent material response do not apply. The effect of the three-dimensional pillar geometry (axisymmetric in the analyses here) and material rate dependence can differ for the three potential instability modes. The finite aspect ratio of the pillars delays the onset of diffuse necking to beyond attaining a maximum nominal stress [37], and material rate dependence delays the onset of both the diffuse and the localized necking [29,33].

The main focus of the simulations here is to quantify the dependence of the mechanical response of the pillars on several parameters that characterize the function $g$ in Eq. (26). The pillar geometry, the constraint-induced deformation inhomogeneity and material rate dependence can all play a role in determining the instability mode and the postinstability response. The relation between the overall stress-strain response of the pillars and the material stress-strain relationship is also of interest since it determines the extent to which a uniaxial tension test can be used to determine material properties.

## Problem Formulation

The formulation of the initial/boundary problem and the numerical solution procedure are given in [22,23] and are briefly outlined here. Further details and references are given in [22,23]. A Lagrangian finite deformation formulation is used, and axisymmetric finite calculations are carried out (cylindrical coordinates $r,z,\theta $ in the reference configuration with all field quantities independent of $\theta $) using the dynamic principle of virtual work.

Initially, the stress free configuration of the pillar occupies $0\u2264r\u2264R0,0\u2264z\u2264H0$. Except for a short ramp up time, a constant velocity $W$ is imposed at $z=H0$. The surface at $z=0$ is fixed to a substrate so that all velocity components vanish there. The remaining boundary conditions are zero traction conditions.

The finite element mesh is comprised of uniform square elements, each consisting of four “crossed” triangles. Time integration is carried out by the explicit Newmark $\beta $ method with a lumped mass matrix. The constitutive update is carried out via a rate tangent method.

The constitutive relation is that of a compressible elastic-viscoplastic solid, as described in [22,23], and is formulated in terms of the Kirchhoff stress $\tau =J\sigma $, where $J$ is the determinant of the deformation gradient and $\sigma $ is the Cauchy (or true) stress. There is no fundamental reason for choosing to phrase the constitutive relation in terms of Kirchhoff rather than Cauchy stress.

where $L$ is the tensor of elastic moduli, $E$ is Young's modulus, $\nu $ is Poisson's ratio, $tr(\xb7)$ denotes the trace, $I$ is the identity tensor, $A:B=AijBji$, and $\tau \u2227$ is the Jaumann rate of Kirchhoff stress.

In Eq. (24)$\u025b\xb70$ is a reference strain rate and $m$ is the rate hardening exponent.

When $h1>0$, $h2<0$, $h3>0$, the parameters $h1$, $h2$, and $h3$ give the sort of hardening-softening-hardening relation that can qualitatively describe the progressive buckling-type deformation mode of VACNTs (see Ref. [22]). The value of $\u025b1$ specifies the plastic strain at which the transition from hardening to softening takes place and $\u025b2$ specifies the strain at which the transition back to hardening occurs. Here we explore a wider range of possibilities including; for example, cases where $h2$ is positive and other cases where $h2=h3$ in which case the value of $\u025b2$ is irrelevant.

Dynamic finite element calculations are carried out because of the instabilities that occur. In some cases, the postinstability response can be difficult, if not impossible, to calculate quasi-statically while including inertia with a material rate provides a physically based regularization.

For a quasi-static response, the values of parameters such as Young's modulus $E$, reference stress $\sigma 0$, pillar length $H$, and radius $R$ do not separately affect the response. The response depends on the values of appropriate ratios and, since the focus here is on quasi-static response, these ratios are chosen to be reported. The fixed parameters (and parameter ratios) are $E/\sigma 0=100$, $\alpha p=0.2$, $\beta p=0.28$, $m=0.02$, and the pillar aspect ratio was fixed at $H0/R0=3$. The imposed velocity $W$ was specified such that $W/H0\u025b\xb70=0.004$.

Unless specified otherwise, the parameters $h1$ and $\u025b1$ are fixed at $h1=24$ and $\u025b1=0.085$, and the only parameters varied in the calculations are $h2$, $h3$, and $\u025b2$ in Eq. (26).

## Numerical Results

### Diffuse and Localized Necking Modes.

All calculations of diffuse and localized necking modes used a $30\xd790$ quadrilateral mesh which gives a uniform mesh of square elements in the reference configuration.

Figure 1 shows the hardness function $g(\u025bp)$ for two materials: material A with $h2=5.0$, $h3=15.0$, and $\u025b2=0.6$ and, for comparison purposes, material B with $h2=h3=5.0$. In terms of the variation of $g$ with plastic strain, for $\u025bp\u2265\u025b2$ material A exhibits hardening-increased hardening, whereas material B has a constant hardening rate $h2=5.0$. The computed overall stress-strain curves for materials A and B are shown in Fig. 2.

Here $U$ is the displacement on the loaded surface at $z=H0$, $F$ is the deformation gradient, $(\u2003)-T$ denotes the inverse transpose, and $\u2207$ is the gradient. Also, $\u025bt=ln(H/H0)$ with $H=H0+U$.

The curves of $\sigma t$ versus $\u025bt$ are nearly identical for materials A and B. Also, the values of the nominal stress $sn$ are essentially identical for both materials up to $\u025bt=0.55$ with a maximum in $sn$ occurring at $\u025bt\u22480.43$. However, for $\u025bt>0.55$, $sn$ increases for material A and decreases for material B. As a consequence, the evolution of the deformation distributions within the two specimens differs significantly.

Figure 3(a) shows the distribution of normalized plastic strain rate $e\xb7=\u025bp/(W/H0)$, for material A at $\u025bt=0.50$, which is a bit past the strain at maximum nominal stress. A very shallow neck has formed at the end $z=H0$ and the maximum value of $e\xb7$ is $0.76$. Subsequently, the plastic strain locally exceeds $0.6$ in the incipient neck and the material there stiffens since $h3=3h2$. The deformation, and therefore the stiffening, propagates down the bar and at the last stage of deformation, $\u025bt=0.76$ shown in Fig. 3(b), most of the bar is uniformly straining with $e\xb7\u22480.4$. The constraint at $z=0$ leads to a strongly inhomogeneous strain rate distribution near that end with a maximum value of $e\xb7$ of about $1$ at $(R0/2)$.

Distributions of normalized plastic strain rate are shown in Fig. 4 for material B at the same two values of $\u025bt$ as in Fig. 3. The plastic strain rate distribution at $\u025bt=0.50$ in Fig. 4(a) is essentially the same as that for material A. However, without the increased hardening of material A, the diffuse neck continues to develop as seen in Fig. 4(b), where the maximum plastic strain rate $e\xb7\u22489.5$ is in the neck. Thus, the additional hardening acts to suppress diffuse necking.

For both materials A and B the true stress continues to increase strongly over the strain range shown while the nominal stress attains a shallow maximum, Fig. 2. We now consider hardening relations for which the true stress attains a maximum. Figure 5 shows the hardness functions $g(\u025bp)$ for material C $h2=h3=1.0$ and material D $h2=h3=0.5$. For comparison purposes $g(\u025bp)$ for material B is also shown.

The overall stress-strain curves for materials C and D are shown in Fig. 6 with the curve for material B from Fig. 2 included for comparison. For materials C and D the overall true stress $\sigma t$, as well the overall nominal stress $sn$, reaches a maximum. For material C a sharp drop in the stress-strain curves occurs at a relatively small strain, $\u025bt\u22480.14$, and the values of $\sigma t$ and $sn$ nearly coincide. For material D, the sharp drop in the overall stress values occurs at a much larger value of $\u025bt$ and there is a clear difference between the values of $\sigma t$ and $sn$ until the sharp drop occurs.

Distributions of normalized plastic strain rate $e\xb7=\u025bp/(W/H0)$ for materials C and D are shown in Fig. 7. For material C considerable necking has taken place before the strain rate has localized into a band, whereas for material D the strain concentration due to the constraint at $z=0$ precipitated localization into a band before any significant necking has taken place. Although dynamic calculations are carried out, material inertia plays a rather small role in the overall response in the calculations for materials A, B, C, and D. For example, for material D the latter stages of straining are characterized by the rapidly localized deformations, the kinetic energy is of the order of $10-4$ times the work input. For materials A and B this ratio is an order of magnitude smaller.

With plastic normality and rate independent material response, localization requires strongly negative hardening while plastic non-normality, as in the constitutive relation here, promotes localization [38]. Rate dependence, as also included in the constitutive relation here, is a stabilizing factor [29] but for slight rate sensitivity the main qualitative features of the rate independent constitutive response are preserved.

From the simple one-dimensional model relation Eq. (15), and for linear hardening, the value of the true stress eventually decreases and can become negative at a sufficiently large strain. The possibilities then include the true stress decreasing to essentially zero or the eventual emergence of a localized deformation mode. The degree of rate dependence and the extent of the deviation from normality play essential roles in determining which of these occurs. The increased hardening at large strains as for material A acts to delay both of these outcomes.

### Propagating Band Mode.

Figure 8 shows the hardness function $g(\u025bp)$ for three materials: the one with $h2=0.5$ and $\u025b2=0.6$ is the same as material D up to $\u025bp=0.6$ and then $h3=15.0$ as for material A. The other two materials have $h2=-3.90$ and $h3=15.0$ but different values of $\u025b2$. With $h2=-3.90$, $g(\u025bp)$ would vanish at $\u025bp\u22480.77$. The minimum value of $g(\u025bp)$ is limited to $0.1\sigma 0$ so that for this material $g(\u025bp)=0.1\sigma 0$ from $\u22480.77$ to $5.0$.

The calculations for a propagating band used a $40\xd760$ quadrilateral mesh. The reason for using this mesh stems from the observation of Ballarin et al. [39] that mesh induced oscillations in the overall stress-strain response occur for a propagating band deformation mode. The amplitude of the oscillations depends on the orientation of the mesh boundaries and it was found that this mesh led to a noticeable reduction in the oscillation amplitude as compared with a square quadrilateral mesh. Also, it was reported in [23] that calculating the overall stress-strain response from the quasi-static principle of virtual work (Eq. (27)) rather than the dynamic principle of virtual work reduced the amplitude of the oscillations but did not significantly affect the mean value of the overall stress.

Figure 9 shows the computed stress-strain response for material E ($h2=-3.90$, $h3=15.0$, $\u025b2=0.6$). There is an initial peak followed by oscillations of both $\sigma t$ and $sn$ about a plateau with a mean value of about $1.4\sigma 0$ for $sn$ and $1.8\sigma 0$ for $\sigma t$ until $\u025bt\u22480.37$ at which point there is a much larger increase in $\sigma t$ than in $sn$. The oscillations continue until $\u025bt\u22480.56$ and vanish at higher strains. Using Eq. (21) with $\tau =g$ gives a Maxwell stress of $1.26\sigma 0$, which is a bit smaller than the mean oscillating value of $sn/\sigma 0$. A uniaxial compression calculation was carried out for material E (not shown here). The deformation mode involved progressive bucklinglike deformations as in [22] but with a different amplitude and wavelength.

The oscillations in the overall stress-strain response are associated with the deformation mode being the propagation of a band of increased strain rate. Figure 10 shows contours of the normalized strain rate $e\xb7=\u025bp/(W/H0)$ at three values of overall strain $\u025bt$. At $\u025bt=0.20$ the band of high strain rate has propagated about $2/3$ the length of the pillar. The strain rate contours at $\u025bt=0.40$ reveal that the reason for the sharp increase in the value of the overall true stress $\sigma t$ is mainly due to the decrease in area that occurs when the high strain rate band has reached at $z=H0$. At $\u025bt=0.50$ the band has propagated down to near the pillar base and the magnitude of the strain rate in the band is smaller than in the earlier stage of deformation at $\u025bt=0.20$. The oscillations in the overall stress-strain curves disappear when band propagation dies out.

Insight into the constitutive features contributing to propagating band propagation and to the source of the oscillations in the overall stress-strain response is given by the response of materials F and G. Figure 11 shows the normalized overall nominal stress $sn$ versus overall strain $\u025bt$ for these two materials. The stress-strain curve for material F exhibits no oscillations, while that for material G exhibits large amplitude oscillations that have minima near zero. The response for material F can be compared with that for material D in Fig. 2. The strong hardening for $\u025bp>0.6$ gives rise to a very different overall stress-strain response. For material F there is a plateau in $sn$ until $\u025bt\u22480.46$ and then a drop associated with the arrival of the deformation band at $z=H0$.

The normalized strain rate distributions for material F (Fig. 12) show a wide plastic strain rate band (compared to the element size) that propagates along the length of the pillar during the strain range that corresponds to the plateau in Fig. 11. The value of $sn/\sigma 0$ on the plateau is smaller than the value of the Maxwell stress obtained from Eq. (21) with $\tau =g$ which is $sM=2.22\sigma 0$. This is consistent with the deformation mode not involving a sharply defined traveling front as assumed for the Maxwell stress.

The plastic strain rate band in Fig. 13 for material G is the width of a deformed element. Each stress oscillation corresponds to the propagation of the strain rate band through a row of elements. The stress increases in a row of elements and then decreases to almost zero (the minimum is kept above zero by requiring that $g(\u025bp)\u22650.1$) as the strains in that element row become extremely large. Then due to the hardening at very large strains the stress increases and deformation localizes in the next row of elements and so on. There is a burst of kinetic energy with the ratio of kinetic energy to work input being of the order of $10-3$ (an order of magnitude greater than for the calculations with localized necking), as the overall stress goes through a minimum and then increases as the deformation band crosses a mesh boundary. This process, in less extreme form, is what gives rise to the oscillations seen in Fig. 9.

A calculation was also carried out for the extreme case of $\u025b2=50$. The overall stress-strain response was the same as that for material D in Fig. 6 except for dynamic ringing when the overall stress approached zero. The dynamic ringing with a decreasing amplitude continued to large strains and the deformation remained concentrated in the elements near $z=0$.

It is worth noting that in the one-dimensional analysis of band propagation in Sec. 2 material inertia plays a key role. The Maxwell stress emerges as an approximate limiting case for quasi-static propagation. The derivation of the Maxwell stress depends on canceling the band propagation speed $Wb$ from both sides of Eq. (20) and so the quasi-static limit ($Wb=0$) is approached but never reached. For material E, as for material G, the ratio of kinetic energy to work input spikes where this ratio is of the order of $10-3$. On the other hand, for material F, which has a broad, slowly traveling strain rate concentration, this ratio remains an order of magnitude smaller with no spikes.

## Concluding Remarks

We have carried out finite element analyses of uniaxial tension of cylinders using a constitutive relation for strain rate-dependent, plastically compressible solids with plastic non-normality. The flow rule was expressed in terms of the effective Kirchhoff stress and the hardness function had three linear regimes. The slope and the extent of the middle linear regime were varied. Depending on the parameters characterizing this regime, the deformation modes that emerge are diffuse necking, localized necking, and/or a propagating band deformation mode.

A diffuse neck forms when the slope of the middle regime $h2$ is sufficiently positive but still small enough for a maximum nominal stress to be attained. When the subsequent high hardening regime is reached, the deformation becomes more uniform (the calculation with $h2=h3$ can be viewed as one in which $\u025b2$ is essentially infinite). In both cases in Fig. 2, the nominal stress (which is what is directly measured in experiments) can be related to the intrinsic material response if the geometry change due to necking is accounted for.

For a sufficiently large negative slope $h2$, the deformations can localize into a narrow band at relatively small strains. The condition for this involves an interplay between the slope of the hardness function, plastic non-normality (destabilizing) and material rate dependence (stabilizing). Nevertheless, when such a localization occurs in uniaxial tension, the strain at which it occurs can give at least a qualitative indication of the form of the hardening response.

The situation when band propagation occurs is more complex. However, at least within the constitutive framework considered here, band propagation is a direct consequence of the hardening-softening-hardening shape of the hardness function. Furthermore, the results for materials E and F indicate that the width of the propagating band is related to the magnitude of the negative slope $h2$. Thus, at least in principle, measurement of the overall stress strain response together with the in situ observation of the propagating band would give considerable insight into the plastic properties of the material. This, of course, presumes that the constitutive framework used in the calculations is appropriate for characterizing the mechanical response.

A numerical issue, pointed out in [39] and explored here, is that when band propagation occurs the stress-strain response computed via a standard finite element method as used here contains mesh induced oscillations so that the value of the propagation stress needs to be estimated from an oscillating response.

Developing a straightforward relation between the measured stress-strain response and material properties of VACNTs may be further hindered by the presence of a density gradient arising from their synthesis [24], which induces a gradient in material properties along the sample axes. This, in turn, leads to an increased apparent hardening as shown in [22]. Preliminary calculations (not shown here) carried out with a gradient in $E$ and $\sigma 0$ for uniaxial tension also showed an increase in apparent hardening due to such a gradient.

The analyses here may not be directly applicable to VACNTs. For example, as noted above, the calculations do not account for the effects of a density gradient. Also, the constitutive formulation of Refs. [22,23] is isotropic while VACNTs are expected to be mechanically anisotropic (different properties along the vertical direction from those perpendicular to it). Nevertheless, based on the results in [22] it seems likely that a hardening-softening-hardening characterization is appropriate for VACNTs and the results here indicate that qualitatively different responses are obtained depending on the character of the hardness function. At present, no experimental data on uniaxial tension of entangled VACNTs is available. Although directly relating uniaxial tension response to intrinsic material properties is difficult, the results here indicate that insight into such properties could be obtained from the uniaxial tensile response of VACNTs (the level, slope, and extent of the stress propagation plateau, and the band width and propagation speed). This provides a motivation for such experiments, which would also serve as a test of the predictability of the constitutive framework in [22,23] as well as providing information for developing an improved constitutive characterization of VACNTs.

## Acknowledgment

The authors gratefully acknowledge the financial support of the Army Research Office through the Institute for Collaborative Biotechnologies (ICB) at Caltech (ARO Award number UCSB.ICB4b).