Effect of Surface Elasticity on the Elastic Response of Nanoporous Gold

: In this work, the effect of surface elasticity on the effective elastic properties of nanoporous gold is studied. To this end, a theoretical framework for surface elasticity effects in submicron-sized solids is implemented as a user-defined finite element subroutine. This allows the use of the theory in large-scale engineering problems. The theory suggests a zero-thickness surface accommodating unique energetic properties and surface tension. For the example of ball-and-stick diamond cubic unit cell structures for nanoporous gold, it is shown that incorporation of surface excess elasticity and surface tension allows prediction of the size effect associated with the change of the surface area-to-volume ratio by capturing, e.g., the increase in the effective Young ’ s modulus and decrease in the effective Poisson ’ s ratio with decreasing ligament diameter, a phenomenon that is not accessible to classical continuum elasticity approaches. DOI: 10.1061/(ASCE) NM.2153-5477.0000126. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/.


Introduction
Cellular solids are formed of interconnected networks of solid struts or plates that make up cell edges or walls, respectively. A class of nanoporous metals can be categorized as open cellular solids, which exhibit a network structure of ligaments with characteristic diameters steplessly adjustable within a 5-500 nm interval; see, e.g., Huber et al. (2014) and references cited therein. Even smaller ligament sizes in the range of 1-2 nm have been reported in the literature (Jin et al. 2008). This results in a highly porous microstructure with up to 30% solid fraction, and a very high surface-area-to-volume ratio; see, e.g., the scanning electron microscope (SEM) image of nanoporous gold with a ligament size of approximately 30 nm in Fig. 1 (Wilmers et al. 2017). Because properties exhibited by the surfaces of the bodies are different from those associated with their interiors, for such nanosized samples accommodating pores, the stored energy in the surfaces becomes comparable to that of the bulk. A classical example is the dramatic increase of the Young's modulus of nanoporous metal samples with reduction of the average ligament diameter (Mathur and Erlebacher 2007;Biener et al. 2006). This kind of size effect is not considered in conventional continuum mechanics estimates.
Following Gibbs's abstraction of material interfaces as a mathematical surface, thus treating surfaces as bidimensional geometrical boundaries of zero thickness, a rigorous mathematical framework to study the mechanical behavior of material surfaces was developed by Gurtin and Murdoch (1975b, 1975a, 1978, which was later generalized to account for interfaces (Gurtin et al. 1998) and incorporation of, e.g., surface curvature; see, e.g., Ogden (1999, 1997) and Chhapadia et al. (2011). In this approach, physical quantities between the surface and bulk are discontinuous. Alternatively, the surface can be treated as a layer with small but finite uniform thickness where the physical quantities undergo a smooth transition, an approach that dates back to van der Waals (Guggenheim 1940). Analytical models, e.g., that adapt Gurtin-Murdoch theory to study static and dynamic behavior of various nanoscale structures such as nanofilms under bending (Lim and He 2004;Lu et al. 2006) and nanobeams under bending (Liu et al. 2011), as well as effective elastic properties of nanoporous gold (Nazarenko et al. 2015), although useful, are limited to simple geometries and boundary conditions. There have been attempts to develop structural elements, e.g., beams, incorporating surface elasticity effects (Liu et al. 2011;He and Lilley 2009). Nevertheless, these involve various simplifications, e.g., the stress state being plane stress, which might not be realistic for various cellular materials, or geometrically linear analysis, in which the effect of residual stresses over the stiffness and frequency of the structure with configurational changes after deformation is disregarded, as noted in Shi et al. (2012).
Combined experimental-numerical studies on the elastic and plastic Poisson's ratio were presented in Lührs et al. (2016) without considering surface elasticity effects. Contributing to the ongoing research on nanoporous metals (Saane et al. 2014;Bargmann et al. 2016;Lührs et al. 2016;Husser et al. 2017), the present work aims at elucidating the important physical consequences of incorporation of surface energy for nanoporous gold. To this end, the framework proposed in Javili et al. (2014) is implemented as an ABAQUS version 6.12 user-defined subroutine where geometrical nonlinearity is consistently incorporated. This makes the theory applicable not only to academic examples but also to more realistic nanostructured materials and large-scale engineering problems. As an example, the authors investigate the effect of surface elasticity on the elastic behavior of nanoporous gold, devising an idealized 3D diamond cubic unit cell structure. On this basis, the size effect associated with an increasing surface-area-to-volume ratio as a result of a decreasing ligament size is investigated.

A Word on Notation
For the repeated indices, lowercase Latin letters i; j; : : : imply summation over the range f1; 2; 3g, whereas for lowercase Greek letters, i.e., α; β; : : : , the range is f1; 2g. Assuming a and b as two vectors and A, B, and C as three second-order tensors, a · b ¼ ½a i ½b i , where ½a i denotes the contravariant components of a and ½b i denotes the covariant components of b. The compo- where E, F , and G represent fourth-order tensors. Following the notation introduced in Javili et al. (2014), any quantity associated with the surface is distinguished from its bulk counterpart ð⋄Þ by the use of a hat, i.e., ð⋄Þ. In this line, Div and Grad represent bulk material divergence and gradient operators and d Div and d Grad their surface counterparts. The expressions ð⋄Þ ⊤ and ð⋄Þ −1 denote the transpose and inverse of ð⋄Þ, respectively. log denotes the natural logarithm.

Differential Geometric and Kinematic Fundamentals
Let B 0 denote the material configuration taken by the continuum body at time t 0 with associated material particles X. Let S 0 ¼ ∂B 0 denote the surface of the body at material configuration with associated material particlesX. The nonlinear deformation maps φ and φ, respectively, transform B 0 and S 0 to the spatial configurations B t and S t ¼ ∂B t at time t > t 0 at which the particle placements are denoted by x for the bulk andx for the surface, which is firmly attached to the bulk; i.e., φ ¼φ on ∂B 0 throughout the motion. The authors define dX ∈ TB 0 (tangent to B 0 ) and dX ∈ TS 0 (tangent to S 0 ) as two material line elements in the bulk and on the surface, respectively. The bulk deformation gradient F, as a linear deformation map associated with φ, and surface deformation gradientF, as a linear deformation map associated withφ, respectively, map dX and dX to the corresponding spatial line elements dx ∈ TB t and dx ∈ TS t with dx ¼ F · dX and dx ¼F · dX. Inverse mappings are possible; viz. dX ¼ f · dx and dX ¼f · dx.
Using the notion of convective coordinates, let fg 1 ; g 2 ; g 3 g denote a covariant basis in the spatial configuration associated with the bulk. The dual basis fg 1 ; g 2 ; g 3 g then satisfies g i · g j ¼ δ j i , where δ j i denotes the Kronecker delta. Analogously, let G i and G i denote covariant and contravariant base vectors in the material configuration for which G i · G j ¼ δ j i . Transformations between covariant and contravariant coordinates in spatial and material configurations are possible via where For the orthonormal Cartesian basis fe 1 ; e 2 ; e 3 g, on the other hand, one has e i ¼ e i to give e i · e j ¼ δ ij . Furthermore, transformation between Cartesian and spatial or material curvilinear coordinates devises Similarly to the bulk, let fĝ 1 ;ĝ 2 g denote a covariant basis in the spatial configuration associated with the surface. The dual basis fĝ 1 ;ĝ 2 g then satisfiesĝ α ·ĝ β ¼ δ β α , where δ β α . Analogously, let G α andĜ α denote covariant and contravariant base vectors in the material configuration for whichĜ α ·Ĝ β ¼ δ β α . Transformations between covariant and contravariant coordinates in spatial and material configurations are possible viâ whereĝ αβ ¼ĝ α ·ĝ β ,ĝ αβ ¼ĝ α ·ĝ β , andĜ αβ ¼Ĝ α ·Ĝ β , and G αβ ¼Ĝ α ·Ĝ β with ½ĝ αβ ¼ ½ĝ αβ −1 and ½Ĝ αβ ¼ ½Ĝ αβ −1 . Reference and current configurations of the bulk with B 0 and B t and the surface with S 0 and S t , with associated covariant bases and metric tensors, are given in Fig. 2. These allow the definition of F ¼ g i ⊗ G i as a full-rank tensor andF ¼ĝ α ⊗Ĝ α as a rank-deficient tensor. Hence, the inverse mappings are defined as f ¼ conventional second-order identity tensors with 3 × 3 matrix representation having unit diagonal entries; hence i ≡ I, andf ·F ¼Î andF ·f ¼î withî ¼ δ α βĝα ⊗ĝ β ¼ i −ĝ 3 ⊗ĝ 3 . Moreover, J ¼ dv=dV andĴ ¼ da=dA with dA, dV representing differential area and volume elements in the reference configuration, whereas da, dv denote their spatial counterparts and Finally, denoting the curvilinear coordinates associated with the bulk and surface ξ i andξ α the following operators are defined: where, by the chain rule following identities, hold (Holzapfel 2000) gradð⋄Þ ¼ Gradð⋄Þ · f and divð⋄Þ ¼ Gradð⋄Þ∶f ⊤ and Linear Momentum Balance and Weak Form Let b 0 andb 0 denote the body force per unit reference volume and analogically the surface force per unit reference area, respectively. Representing the first Piola-Kirchhoff stress tensor with P, the balance of linear momentum is postulated (Javili et al. 2014) for the bulk and surface as N corresponds to the normal to S 0 at the reference configuration. Multiplying both parts of each equation with corresponding test functions δφ and δφ and applying the extended divergence theorems (Steinmann 2008) while assuming that the surface is firmly bonded to the bulk (i.e., φ ≡φ on ∂B 0 ) yields the associated weak form Z As a consequence of angular momentum balance, one has P · F ⊤ ¼ F · P ⊤ over the bulk andP ·F ⊤ ¼F ·P ⊤ over the surface (Javili et al. 2014).

Finite Element Discretization
This step is followed by finite element discretization over the bulk and surface of the domain, where surface elements are selected to comply with the underlying bulk (Javili and Steinmann 2010) where N B 0 and N S 0 denote the total number of bulk and surface elements to make up the total number of elements N el ¼ N B 0 þ N S 0 in the finite element mesh.
The three-dimensional and two-dimensional interpolations with corresponding shape functions N I andN I as functions of natural coordinates, ξ ¼ ðξ 1 ; ξ 2 ; ξ 3 Þ for the bulk andξ ¼ ðξ 1 ;ξ 2 Þ for the surface, make up the geometry of the bulk and surface elements, respectively, to give where N node andN node denote the total number of nodes in bulk and surface elements, respectively. Substituting spatially discrete forms into the weak form, the authors carry out the following residual vector R I e;⊎ for a local node I: where ⊎ denotes contributions from the connected bulk element n through R I e;n and the connected surface element k, fully complying with the underlying bulk element throughR I e;k , where Analogically, the local tangent stiffness for a local node couple fI; Jg, with contributions from the connected bulk element n through K IJ e;n and the connected surface element k, fully complying with the underlying bulk element throughK IJ e;k , is then computed as where ½K IJ e;n AC ∶ ¼ The fourth-order bulk and surface elasticity tensors A andÂ are computed using Reference and current configurations of the bulk with B 0 and B t and the surface with S 0 and S t , with associated covariant bases and metric tensors; the bulk Lagrangian and Eulerian covariant base vectors G i and g i are tangential to the Lagrangian and Eulerian material lines along the curvilinear coordinates Θ i and θ i ; the surface Lagrangian and Eulerian covariant base vectorsĜ α andĝ α are tangential to the Lagrangian and Eulerian material lines along the curvilinear coordinatesΘ α andθ

Hyperelastic Bulk and Surface Potentials
Restricting the theory to the case of material isotropy, in the current hyperelastic framework, P andP are derived from the bulk energy density per unit reference volume W 0 ðFÞ and surface energy density per unit reference areaŴ 0 ðFÞ; viz λ and μ,λ andμ are the Lamé constants associated with the bulk and the surface. Moreover,γ denotes the surface tension.
With these potentials, one ends up with the elastic stress definitions Finally, with the use of Eq. (16), one computes the following tangents: which concludes the theory. For further details, the reader is referred to Javili et al. (2014) and the references therein. In view of Eqs. (13) and (15), two separate user-defined element subroutines (UELs) are implemented in ABAQUS, one for the bulk element and one for the surface element. Eq. (13) then amounts to −1× RHS and Eq. (15) for AMATRX in ABAQUS terminology. This makes the theory applicable not only to academic examples but also to more realistic large-scale engineering problems.

Finite Element Model and Material Properties
In order to study surface effects on the effective elastic response of nanoporous gold, an idealized representative volume element (RVE) as a ball-and-stick assembly with a face-centered cubic diamond lattice is selected (Huber et al. 2014). The gold fraction is kept constant at φ ≃ 9.55% to eliminate its influence on the effective properties. The influence of φ was analyzed, e.g., in Bargmann et al. (2016). Here, the authors observed a positive correlation between φ and the effective elasticity modulus. The increase of macroscopic yield stress with increasing φ was also reported. Despite their inherent structural anisotropy, diamond cubic unit cell idealizations and their geometrically perturbed forms have previously been used in the literature for the purpose of modeling nanoporous gold and their composites, see, e.g., Huber et al. (2014) and Roschning and Huber (2016). This makes the model development tractable, especially as compared to modeling real microscale topologies (cf. Fig. 1), and hence allows deriving plausible inferences pertaining to force transmission along randomly distributed struts, deformation, and, as in the current study, excess surface elasticity effects in the actual microstructures. The bulk and surface of the gold skeleton in the RVE are discretized into 34,304 trilinear hexahedral solid elements and 11,520 bilinear quadrilateral surface elements; cf. Fig. 3. To mimic smooth transitions at ligament junctions, the nodes are modeled with a slightly larger node radius.
For both the bulk and the surface elasticity, the material behavior is assumed to be isotropic with corresponding potentials in Eq. (17). Material parameters together with the corresponding references are listed in Table 1. The excess surface elasticity is reported in Elsner et al. (2017). The determination of surface tension in solid materials is still a challenging task. In this respect, the authors revert to the published data for gold in Tyson and Miller (1977) and Kumikov and Khokonov (1983) which, respectively, suggest magnitudes as 1.5 J=m 2 resp. 1.41 J=m 2 . Admissible limits of the surface parameters were studied in Javili et al. (2012). Using the linearized theory of elasticity, the parameters should fulfill −1 < ν < 1=2 and −1 <ν < 1. Moreover, strong ellipticity at the bulk imposes μ > 0 and λ þ 2μ > 0, whereas for the surface, it imposeŝ μ > 0 andλ þ 2μ > 0. Pointwise stability in the bulk leads to μ > 0 and λ þ μ > 0, whereas at the surface, it gives rise toμ > 0 and λ þμ > 0. Finally, boundary-complementing conditions hold if μ > 0 andλ þ ½2λ þ 3μ=½λ þ 2μμ > 0. All parameters identified Fig. 3. FE-discretization of (a) bulk; (b) surface; generation of the surface model was done in ABAQUS using the built-in method Convert Solid to Shell and used in the parametric studies of the current study comply with these limits. In the following, the authors start by investigating the response of the selected diamond cubic unit cell in the absence of surface elasticity effects. Then, the influence of the ligament diameter D on the effective elastic response is investigated considering the interval D ∈ ½1; 50 nm. For the stationary volume fraction φ, the ligament diameter controls the ratio of surface area to bulk volume. This is clearly exemplified in Fig. 4 where unit cell combinations of different ligament diameters are embedded in identically sized cubes. As the RVE size gets smaller, the ligament dimensions and therefore the ligament diameter D get smaller, and the surface-area-to-volume ratio increases exponentially. Particular ratios under investigation are in the interval from A=V ≃ 65.65 nm −1 to A=V ≃ 3282.27 nm −1 . This comprises the source of the size effect associated with surface elasticity. The authors also independently evaluate the influences of surface excess elasticity and surface tension on the effective elastic material properties such as Young's modulus and Poisson's ratio. Forλ,μ, andγ, each parametric variation is realized considering the interval ½0; 15 J=m 2 . This interval encloses the experimentally obtained magnitudes for the parameters previously given in Table 1. For the sake of completeness, unreasonably large values are also considered, up to 300 J=m 2 . Unlike surface excess elasticity, surface tension alters the self-equilibrium state of the material.
Diamond cubic unit cells are anisotropic. Considering the point symmetry operations of rotation about an n-fold axis, reflection, inversion, rotation-reflection, and rotation-inversion, the point group of the diamond structure is the full octahedral group O h with 48 operations (Srivastava 1990). Hence, with reference to Fig. 5, the authors consider uniaxial compressive loading along crystal directions ½001, ½011, and ½111 of the diamond cubic unit cell-based microstructure under periodic boundary conditions. Because of symmetry, the associated response relates to the whole family of directions, h100i, h110i, and h111i. Corresponding effective Young's moduli are denoted as E Ã n , where n denotes one of the directions, ½001, ½011, or ½111, as the direction of uniaxial (compressive) stress. Using the notation ν Ã ðn 1 ;n 2 Þ , the Poisson's ratio is associated with two orthogonal directions, n 1 and n 2 . Here, n 1 represents the direction of uniaxial (compressive) stress, whereas n 2 represents the direction for which the Poisson's effect is observed. Each effective elastic property is computed relative to an equilibrium configuration by perturbing the state with 0.1% of strain in the loading direction under the uniaxial state of macrostress. For the simulations without surface tension, this corresponds to the undeformed configuration. With surface tension, a selfequilibrium state is yet to be attained where the structure experiences deformation that aims at minimizing the surface energyγ prior to the subsequent loading. This self-equilibrium configuration constitutes the reference configuration based on which the elastic properties are computed. The homogenization procedure used in the determination of effective elastic properties considering surface effects is given in the Appendix.

Directional Dependence of Elastic Response Without Surface Elasticity
Imposed loads on the diamond cubic unit cell are supported by means of forces transmitted through struts forming the microstructure. Because they have different orientations, in effect, even a macroscopically axial force manifests itself as axial, bending, and shear force over elements. In this subsection, the authors shall explore the directional dependence of elastic structural behavior of the diamond cubic unit cell in the absence of surface elasticity. Fig. 6(a) depicts the directional dependence of the effective Young's modulus in 3D. The direction family h111i constitutes the stiffest direction accommodating a dominating axial mode of  Table 2.
Here, the subscript "0" denotes the absence of surface elasticity. Accordingly, the effective Young's moduli for the directions ½001 and ½011 are approximately 38 and 71% of that for direction ½111, respectively giving E Ã 0;½001 < E Ã 0;½011 < E Ã 0;½111 . In Huber et al. (2014), a considerable decrease in the effective Young's modulus with increasing structural disorder was demonstrated using RVEs made up of multiple unit cells with randomly shifted connecting nodes. It is anticipated that such a disorder will also reduce the directional dependence of the elastic properties.
The interval ½0.1; 0.4 is common for an effective elastic Poisson's ratio of cellular solids such as foams. Based on the load and unload tests conducted over nanoporous gold samples with ligament diameters of D ¼ 50, D ¼ 120, and D ¼ 150 nm, the effective elastic Poisson's ratio was essentially found to be independent of D, starting at approximately 0.20 for the undeformed material with a weakly increasing trend with densification during plastic compression (Lührs et al. 2016). Fig. 6(b) depicts the directional dependence of the effective Poisson's ratio for the diamond cubic unit cell considering uniaxial stress states along crystal directions ½001, ½011, and ½111. Specific values are listed in Table 2. For loading directions ½001 and ½111, the effective Poisson's ratio is observed to be positive where the authors observe lateral expansion in response to compression. Moreover, the Poisson's ratio is invariant for rotations with respect to the loading direction for each ½001 and ½111. On the other hand, for loading direction ½011, the computed Poisson's ratios along directions perpendicular to ½011 oscillate within the interval ½−0.0556; 0.8153. The minimum of this interval complies with the finding that nearly 70% of the cubic elemental metals are reported to show a negative Poisson's ratio if stretched along the ½011 direction (Baughman et al. 1998). The thermodynamical limitations on isotropic materials require ν ∈ ½−1; 0.5. Because the current microstructure is anisotropic, this limit over the effective elastic Poisson's ratio does not apply. In other words, ν Ã 0;ðn 1 ;n 2 Þ ≃ 0.8153 > 0.5 is theoretically admissible. The difference between these experiments and the current simulation results is attributed to randomness and the connectivity of the struts in the actual material. The crucial role of disorder introduced into ligament networks formed of initially diamond cubic building blocks was reported in Roschning and Huber (2016), which makes it possible to simulate effective Poisson's ratios within the interval ½0.05; 0.45.

Surface Area-to-Volume Ratio as a Source for Size Effect in Nanoporous Materials
Uniaxial compression simulations are realized for the ligament diameter D variation spanning the interval [1,50] nm for three different loading orientations, ½001, ½011, and ½111. In each case, the surface parameters are taken as the ones listed in Table 1.
The effect of D on the effective Young's modulus is represented in Fig. 7(a). The effect is reported as a fraction of the effective Young's modulus in the absence of the surface effects denoted by E Ã 0;n in the corresponding loading direction n as given in Table 2. There is considerable agreement in the plots for different directions. Hence, the change in the ligament diameter evenly affects the effective Young's modulus in each direction. The elastic stiffness of the nanoporous gold systematically increases with decreasing ligament size, whereas with increasing ligament size, the effects rapidly saturate. It was demonstrated in Altenbach et al. (2013) that, in comparison with bulk material, the presence of surface stress leads to the increase of stiffness of nanosized specimens, whereas residual stresses may decrease or increase the effective stiffness. In accordance with the findings noted in, e.g., Altenbach et al. (2013Altenbach et al. ( , 2011, the analysis forλ ¼μ ¼ 0 and two different values ofγ withγ > 0 andγ < 0 shows that, whereas forγ > 0, an increase of the effective Young's modulus with decreasing ligament diameter is observed, forγ < 0, the effective Young's modulus decreases with decreasing ligament diameter. A significant change is observed only below D < 10 nm, where the steepness of the trend becomes especially recognized. This is in complete agreement with the analytical continuum analysis of the effective stiffness of nanowires under tensile, bending, and torsional loading, which reveals that surface excess elasticity becomes more significant at diameter sizes smaller than 10 nm (Elsner et al. 2017). At D ¼ 1 nm, the highest stiffness is observed, which is nearly 70% larger than that of the microstructure with ligament diameter D ¼ 50 nm. This corresponds to the size effect associated with the surface elasticity, which becomes dominant if the ratio of surface-area-to-volume is considerable; i.e., the porous microstructure becomes very fine. Hence, once the surface-area-to-volume ratio tends to zero, associated surface effects diminish and the overall response of the structure tends to that of the bulk elasticity.  . 6. Directional dependence of (a) effective Young's modulus E Ã 0;n ; (b) effective Poisson's ratio ν Ã 0;ðn 1 ;n 2 Þ for selected diamond cubic unit cell microstructure for D ¼ 50 nm; the minimum and maximum effective Young's moduli are along the family of directions h100i and h111i, respectively; the plot in (a) is normalized with respect to the maximum value; for uniaxial stress along ½001 and ½111, effective Poisson's ratios are invariant of the direction of observation; for uniaxial stress along ½011, a directional dependence is due where Θ n 2 ¼ 0 for n 2 ¼ ½100 and Θ n 2 ¼ π=2 for n 2 ¼ ½011 The effect of D on the effective Poisson's ratio is represented in Fig. 7(b). Unlike before, the effect is reported in terms of the difference from the Poisson's ratio in the absence of surface elasticity effects denoted by ν Ã 0;ðn 1 ;n 2 Þ ðD ¼ 1 nmÞ in the corresponding direction of uniaxial stress n 1 and direction of observation n 2 . As Fig. 7(b) depicts, except for the results for loading direction n 1 ¼ ½011 and sampling direction n 2 ¼ ½001, where an opposite trend is observed, in all loading and sampling directions, increasing D results in an increase in the Poisson's ratio. The trend is steepest at the lowest diameter D ¼ 1 nm, where a rapid saturation is observed. The observed changes in ν Ã lie below 5%, which will make corresponding experimental measurements difficult.
Considering uniaxial compression of D ¼ 5 nm samples with 0.1% axial strain along loading direction ½001 and comparing the equivalent von Mises stress σ vM distributions for the cases of surface elasticity-free response, i.e.,λ ¼μ ¼γ ¼ 0, cf. Fig. 8(a), and the response with surface excess elasticity but without surface tension, i.e.,λ ¼ −2.0 J=m 2 ,μ ¼ 3.5 J=m 2 , andγ ¼ 0, cf. Fig. 8(b), the authors observe a higher stress development over the bulk for the case without surface elasticity. Specifically, the maximum stress value reduces from ≃84 MPa to ≃77 MPa at a strain level of 0.1% with the inclusion of excess surface elasticity. This anticipated result is because of the redistribution of loads over the surface and the bulk with attained surface stiffness.
Consideration of the surface tension changes the picture drastically. Consideringγ ¼ 1.5 J=m 2 withλ ¼ −2.0 J=m 2 andμ ¼ 3.5 J=m 2 for D ¼ 5 nm samples results in a self-equilibrium state shown in Fig. 9. Macroscopic Cauchy stress at selfequilibrium configuration B e is zero; that is, M σ e ¼ 0 with hσ e i B e ¼ −hσ e i S e . The structural symmetry of the diamond cubic unit cell leads to M hσ e i ¼ −ξi. Current computations for D ¼ 5 nm and considered parameters result in ξ ¼ 56.41 MPa. Thus, a macroscopic hydrostatic compression is created over the bulk byγ > 0. Internal (von Mises) stress development reaches ≃580 MPa at the ligament centers, as demonstrated in Fig. 9(a). Close to junctions, a stress concentration region is observed where even higher magnitudes are reached. Although the mechanisms by which plastic flow is initiated in very thin structures is still a point of debate, the reason for a smaller elastic limit with increasing size is associated mainly with the presence of dislocations. In effect, as at least one dimension of the crystals structures falls below the order of a micron or less, the elastic limit increases. At the limit, crystals with sizes smaller than the Frank network dimensions should be perfect (Friedel 1964). For perfect crystals, the theoretical shear (a) (b) Fig. 7. Change of (a) effective Young's modulus E Ã n 1 for compression along n; (b) effective Poisson's ratio ν Ã n 1 ;n 2 in sampling direction n 2 for compression along n 1 as a function of D ∈ ½1; 50 nm with surface elasticity parameters given in Table 1; the size effect on the effective elastic properties associated with a high surface-area-to-volume ratio is clearly evident, especially for D < 10 nm and for all loading directions, with absf½E Ã n ð50Þ − E Ã 0;n =E Ã 0;n g < 1.5% for all n and absf½ν Ã ðn 1 ;n 2 Þ ð50Þ − ν Ã 0;ðn 1 ;n 2 Þ =ν Ã 0;ðn 1 ;n 2 Þ g < 0.15% for all n 1 and n 2 ; even for D ¼ 50 nm, the surface effects are hardly recognized because of the small surface-area-to-volume ratio Fig. 8. von Mises stress distribution over bulk elements of the nanoporous gold structure with ligament diameter D ¼ 5 nm for compression along ½001: (a) deformed configuration after 0.1% strain and for the surface-independent case, i.e.,λ ¼μ ¼γ ¼ 0; (b) deformed configuration after 0.1% strain for the caseλ ¼ −2.0 J=m 2 ,μ ¼ 3.5 J=m 2 , andγ ¼ 0; the results show a reduced stress level over the bulk when excess surface elasticity is active stress τ 0 , that is, the stress necessary to cause plastic slip, is known to be on the order of μ=10. In the current study for gold, the authors compute the theoretical elastic limit (shear stress) as τ 0 ¼ μ=10 ≃ 2743 MPa. As the size of the sample gets smaller, the observed maximum shear stress in the models gets larger. The assumption of elasticity will be valid only below the theoretical elastic limit. For nanoporous gold ligaments with 15 nm diameter, a yield stress of 1.5 GPa was measured because of the presence of very few dislocations or sources in the bulk (Volkert et al. 2006). Thin cylinders undergoing large elastic bending under stresses on the order of the theoretical elastic limit were demonstrated in, e.g., Friedel (1964) and the references therein. With surface tension, because of its inherent symmetry, the diamond cubic unit cell structure tends to shrink isotropically, with macroscopic spatial logarithmic strain of 1=2 log ÀÂ M F e · M F ⊤ e ÃÁ ≃ −0.015i, as demonstrated in Fig. 9(b), where M F e corresponds to the deformation gradient computed at the equilibrium configuration.

Influence of Surface Elasticity and Surface Tension
Considering microstructures with ligament diameter of D ¼ 10 nm, uniaxial compression simulations are realized for each parameter variation forλ,μ, andγ for three different loading orientations ½001, ½011, and ½111. In each case, the surface parameters except for the varied one are set to zero. Corresponding effective elasticity properties are computed and displayed in this section.
The effect of the selected surface parameter on the effective Young's modulus is represented in Fig. 10(a) for variations ofλ, μ, andγ spanning the interval ½0; 15 J=m 2 . The effect is reported as a fraction of the effective Young's modulus in the absence of the surface effects denoted by E Ã 0;n in the corresponding loading direction n, as given in Table 2. As the plotted curves demonstrate, for all parametric variations and loading directions, an increase of the selected parameter results in an increase of the effective Young's modulus. For the selected range of parametric variations, a linear dependence of E Ã n on the selected parameter can be judged for all n. The trends show that the increase in E Ã n under the influence ofλ ¼ 15 J=m 2 does not go beyond 4%, whereas forγ ¼ 15 J=m 2 , a nearly 17% increase is observed relative to the case without surface effects. Still, E Ã n is recorded to be most sensitive to the change of μ, where forλ ¼ 15 J=m 2 , E Ã n =E Ã 0;n reaches 1.25. The recorded change in E Ã n with a change inλ seems to be maximum for direction ½111 and the minimum for direction ½001. For the parametric Fig. 9. von Mises stress distribution for bulk elements of the nanoporous gold structure with ligament diameter D ¼ 5 nm: (a) self-equilibrium state after applied surface tension usingγ ¼ 1.5 J=m 2 ; (b) self-equilibrium state with macroscopic spatial logarithmic strain of 1=2 log ÀÂ M F e · M F ⊤ e ÃÁ ≃ −0.015i, with increased deformation scale factor of 20 for better visualization of the shrinkage behavior; compared to cases without surface tension, a drastic increase in the microscopic stress level is observed; the macroscopic Cauchy stress at the self-equilibrium configuration B e vanishes with M σ e ¼ 0 where hσ e i B e ¼ −hσ e i S e ; with the structural symmetry of the diamond cubic unit cell, the authors compute hσ e i B e ≃ −56.41i MPa (a) (b) Fig. 10. Change of (a) effective Young's modulus E Ã n 1 for compression along n; (b) effective Poisson's ratio ν Ã n 1 ;n 2 in sampling direction n 2 for compression along n 1 as a function ofλ;μ;γ ∈ ½0; 15 J=m 2 for D ¼ 10 nm; for the selected range, E Ã increases linearly with increasingλ;μ andγ for all three loading directions where the influence ofμ is the maximum, with up to an approximately 25% increase, and that ofλ the minimum, with up to an approximately 4% increase; for the effective Poisson's ratio, although generally a decreasing trend is observed with increasing surface elasticity parameters, the change does not go beyond 1.5% variations ofμ andγ, this trend is inverted, where the recorded relative change in E Ã n seems to be the maximum for direction ½001 and the minimum for direction ½111. This interesting output can be attributed to the acting loading mechanism over the struts and its association with the surface parameter. The value ofλ is a measure of surface-area compressibility, and it has a confining effect on the bulk because of surface area conservation. The constantμ is a measure of surface shear resistance. For direction ½111, an axial mode of deformation prevails over the struts, whereas for direction ½100, bending and shearing modes dominate.
The effect of the selected surface parameter on the effective Poisson's ratio, that is, the ratio of the macroscopic strain component representing effective lateral expansion to that of contraction in the direction of loading, is represented in Fig. 10(b) for variations ofλ,μ, andγ spanning the interval ½0; 15 J=m 2 . The effect is reported in terms of the difference from the effective Poisson's ratio in absence of the surface effects denoted by ν Ã 0;ðn 1 ;n 2 Þ . As plotted curves depict, the increase in surface parameters generally results in a decrease in the Poisson's ratio ν Ã ðn 1 ;n 2 Þ with a slight nonlinearity. The only exceptions occur along n 1 ¼ ½011 and n 2 ¼ ½100 for variations ofλ and along n 1 ¼ ½011 and n 2 ¼ ½011 for variations ofμ, in each of which a slight increase in ν Ã ðn 1 ;n 2 Þ is observed.
However, as compared to its influence on Young's modulus, surface elasticity affects Poisson's ratio rather marginally where introduced change is limited approximately to 1.5%. Hence, experimental measurement of this phenomenon is anticipated to be a difficult task. The overall influence on the Poisson's ratio is not straightforward to interpret because it includes the effect of lateral spread of the struts as a bulk with their own Poisson's effect and also the structural spread. Still, in agreement with the previous findings, the largest effect is due toμ andγ and the smallest due toλ.
The aforementioned proportionality of E Ã n to the surface parameters ceases to exist for larger ranges; cf. Fig. 11(a), which shows studies up to 300 J=m 2 magnitude of surface parameters for the loading direction ½001. The response concaves up forγ and down forλ andμ, but nevertheless, the monotonically increasing trend is preserved. As seen, at approximately 50 J=m 2 , the influence of surface tensionγ reaches and becomes larger than that ofμ. However, as compared to the experimentally obtained values, surface parameters as large as 300 J=m 2 seem unreasonable. Aŝ μ;γ → 50 J=m 2 , the effective Young's modulus reaches nearly twice the value of the effective Young's modulus without surface effects. This means that gained stiffness because of surface effects alone reaches that of the bulk skeleton. Again, the influence ofλ is quite small as compared to the other two parameters. Analytical derivations reported in Javili et al. (2015) concerning voided matrices with surface effects show that the overall behavior of the material becomes independent of the surface parameters at the strong-surface-elastic limit. This is attributed to the rigid inclusion-like behavior of extremely strong surfaces. Because in the current work, the microstructure corresponds to an open cell structure, the ligament surfaces continuously contribute to the load-carrying mechanism; hence, surface parameters determine the material response at the strong-surface-effect limit.
Unlike the current observations for Young's modulus, the recorded monotonic trend ceases to exist for larger ranges of the selected surface parameter. As displayed in Fig. 11(b), for ½001, the decrease in the Poisson's ratio with increasing surface parameters is only up to approximately 30 J=m 2 ofλ, 100 J=m 2 ofμ, and 180 J=m 2 ofγ. A further increase of the surface parameter results in an increase of the Poisson's ratios, where the steepest slope is observed forγ and the mildest forμ.

Conclusions
The influence of surface elasticity on the effective elastic properties of nanoporous gold is studied. For this purpose, the theoretical framework for surface elasticity presented by Javili et al. (2014) is implemented in ABAQUS as a user-defined element subroutine in addition to the user-defined element subroutine for bulk elasticity, which made the corresponding theory applicable to large-scale engineering problems.
The impact of surface excess elasticity and surface tension on the effective elastic properties of nanoporous gold was extensively analyzed considering uniaxial compressive loading along the crystal directions ½001, ½011, and ½111 of diamond cubic unit cellbased microstructure in the three-dimensional setting and under periodic boundary conditions. On the basis of the findings, the following conclusions are drawn: • The effective Young's modulus E Ã exhibits a strong dependence on the surface Lamé constantμ and on the surface tensionγ.
Less, but still substantial, influence is obtained from the surface Lamé constantλ. In all loading directions, a monotonic increase of E Ã is observed with increasing surface parameters no matter (a) (b) Fig. 11. Change of (a) effective Young's modulus E Ã ½001 for compression along ½001; (b) effective Poisson's ratio ν Ã ½001;½001 ⊥ for sampling direction ½001 ⊥ and for compression along ½001 as a function of surface parametersλ;μ;γ ∈ ½0; 300 J=m 2 for D ¼ 10 nm; the nonlinearity of the dependence of E Ã on surface parameters is obvious; a stronger influence ofμ andγ than that ofλ is observed where, for • > 50 J=m 2 , γ dominates; as for the effective Poisson's ratio, the trends become highly nonlinear with increasing surface elasticity parameters; still, the change does not go beyond 2.5% which one is selected. Also, the change of the sign ofγ inverts the observed size dependence of the effective Young's modulus for vanishing Lamé constants; • The impact of surface elasticity on the effective Poisson's ratio ν Ã is found to be more complex, but marginal. Yet, again, the influence ofμ andγ is higher as compared to that ofλ; and • A unique size effect associated with surface elasticity showing a clear surface-area-to-volume ratio dependence controlled by the ligament diameter is observed. In the presence of surface elasticity, either in terms of surface excess elasticity or surface tension, lowering the ligament diameter beyond 10 nm strongly amplifies the effective Young's modulus and suppresses the Poisson's ratio in all loading and probing directions. This work contributes to the understanding of the sizedominated processes in nanoporous materials whose microstructural design grants improved and exciting mechanical properties. The challenges posed by such advanced and complex functional materials enforces combined experimental, theoretical, and numerical studies. For instance, the analytical continuum analysis of the effective stiffness of nanowires under tensile, bending, and torsial loading in Elsner et al. (2017) reveals that surface excess elasticity becomes more significant at diameter sizes smaller than 10 nm. That conclusion is in agreement with the numerical results presented in Fig. 7. In fact, the regime around 10 nm might be understood as a rough upper limit at which surface effects takes over and hence should be taken into account.

Appendix. Application of Periodic Boundary Conditions and Determination of the Effective Elastic Properties through Homogenization
Application of periodic boundary conditions is realized via periodically located nodes A and B through the operation M F · ½X A − X B ¼ ½x A − x B by fixing the model at an arbitrary point on the x-, y-, and z-axes in order to avoid rigid body displacements. Here, M F denotes the macroscopic deformation gradient representing the deformation of the unit cell. In the absence of surface tension, compression in the z-direction is considered with ½ M F 33 ¼ 1 þ ½ M Δϵ 33 . Here, ½ M Δϵ 33 ¼ ½ū 3 =L 0 denotes the macroscopic engineering strain increment in the direction of the load based on the prescribed displacement ½ū 3 ¼ ½x 3 − ½X 3 and L 0 is the initial edge length of the unit cell. All other components of M F are computed such that corresponding first Piola-Kirchhoff stress components vanish. Surface tension results in a uniform shrinkage or expansion (depending on the sign ofγ) of the nanostructure without any macroscopic shear deformation. Hence, for computations accounting for surface tension, the authors apply, like before, a strain increment ½ M Δϵ 33 , but this time considering the deformed shape under surface tension as the reference state, with edge length L 1 . In all the demonstrated examples, the elastic material properties are computed at ½ M Δϵ 33 ¼ 0.1% with respect to the corresponding reference equilibrium configuration.
The computation of the effective Young's modulus E Ã e 3 in the direction of compression (z-direction) is realized using E Ã e 3 ¼ ½ M Δσ 33 =½ M Δϵ 33 , where M Δσ ¼ M σ t − M σ e denotes the increment of the homogenized Cauchy stress tensor as a difference of its magnitude at the current configuration M σ t from M σ e , where the subscript e denotes the reference equilibrium configuration B e . Henceforth, the authors drop t but keep e with σ as a subscript for convenience. Micro-to-macro transformations for Cauchy stresses use the following sum of average bulk and surface stresses at the current configuration: where V t = volume of the unit cell in the current configuration; and σ andσ, respectively, denote bulk and surface Cauchy stress tensors with σ ¼ ½1=JP · F andσ ¼ ½1=ĴP ·F. One can also revert to the procedure explained in Kouznetsova et al. (2001). Considering only the loading due to surface tension withγ and applying homogenization at the self-equilibrium configuration, the authors reach where V e = volume of the unit cell in the reference equilibrium configuration. Because of equilibrium, the previous computation yields M σ e ¼ hσ e i B e þ hσ e i S e ¼ 0. Furthermore, with the structural symmetry of the considered diamond cubic unit cell, the authors reach with ξ > 0 forγ > 0. That is, for the current applications, both hσi B e and hσi B e are spherical tensors with only hydrostatic stress components. Thus, the macroscopic Cauchy stress tensor vanishes, and M Δσ → M σ which can be found through homogenization of the Cauchy stresses at the current configuration. For homogenization in the reference configuration, one uses M P ¼ hPi B 0 þ hPi S 0 , where hPi B 0 ¼ ð1=V 0 Þ R B 0 PdV and hPi S 0 ¼ ð1=V 0 Þ R S 0P dA, where V 0 is the volume of the unit cell in the reference configuration; see, e.g., Javili et al. (2015).