This paper derives and demonstrates a one-dimensional acoustic metamaterial homogenization method. The homogenization method uses a multiple-scales approximation with Hamilton's principle, a weak-form representation of the dynamic equation. While the multiple-scales approximation makes the predicted effective material properties of this method inexact, the method is shown to be highly versatile. Analytical and numerical examples are given showing the ability of the homogenization method to account for viscosity and finite-amplitude effects.

## I. INTRODUCTION

Metamaterial research, or the effort to create synthetic composite materials with specific, often exotic, effective properties, relies upon homogenization methods. Homogenization methods are ways to replace inhomogeneous systems with an effective homogeneous system for analysis purposes. Many homogenization methods for acoustic and elastic metamaterials have been proposed. These methods range from analytical schemes such as scattering-based methods^{1,2} and statistical mechanics (see, for example, Ref. 3) and numerical methods such as self-consistent models,^{4} to experimental methods.^{5,6} The majority of the theoretical methods describe micromechanical behavior using a strong-form, or differential, representation rather than a weak-form, or integral, representation. There are many reasons for using a strong-form representation. The physical interpretation of strong forms is relatively easy to understand, and most existing analytical descriptions of fluid and solid mechanics are presented in differential form. On the other hand, weak-form representations, such as Hamilton's principle, have traditionally been used in high level derivations such as in the Lagrangian and Hamiltonian formulation of mechanics and, more recently, in finite element-based numerical methods.

The fact that integration inherently considers finite-sized domains makes using a weak-form representation of the microstructure attractive from a homogenization standpoint. In fact, a statistical mechanics analysis of a perfect gas proceeds from Hamiltonian mechanics. A more recent example is the work by Andia *et al.*,^{7} where a Lagrangian-mechanical approach was developed with the intention to incorporate molecular dynamics into their homogenization scheme. Their work focused primarily on periodic systems, and obtains effective properties by taking the variation of the microscopic Lagrangian with respect to the averaged (macroscopic) displacement.

The purpose of this paper is to present a new, multi-scale homogenization method based on Hamilton's principle for mechanical media. Hamilton's principle is, as mentioned above, a weak-form representation and is equivalent to Newtonian mechanics for continuous media. An important feature of Hamilton's principle is the ease with which constraints, such as constitutive equations and complicated geometries, may be imposed.

The outline of this paper proceeds as follows. Section II presents a basic introduction to Hamilton's principle and then derives the new homogenization method. Examples that illustrate aspects of the new homogenization method are then given in Sec. III. The examples include accounting for viscous damping and considering homogenization of material properties associated with finite-amplitude propagation. Finally, Sec. IV summarizes the paper.

## II. HAMILTON'S PRINCIPLE AND ENERGETIC HOMOGENIZATION

### A. A brief introduction to Hamilton's principle

Newton's second law for one-dimensional acoustic systems may be cast in weak form by Hamilton's principle, which takes the form^{8,9}

where *δ* denotes the variation, *S* is the cross-sectional area (assumed to be constant throughout the analyses below), $Ekin$ is the kinetic energy density of the system, $Epot$ is the potential energy density of the system, *F* is the sum of all body force densities (which may be a function of space, time, and acoustic fields), and *u* is the particle displacement. The difference of the kinetic and potential energy densities is also called the Lagrangian density. Hamilton's principle is an optimization problem, and so additional constraints may be added via a Lagrange multiplier. For example, incorporating the constraint *C* = 0, where *C* is some function of the acoustic variables, may be done by adding the term $S\mu C$ to the Lagrangian density, yielding a modified equation,

Here, the yet-to-be determined function *μ* acts as a Lagrange multiplier.

For homogeneous and linear acoustic systems, the kinetic and potential energy densities may take the form, respectively, of

where *ρ* is the mass density, *β* is the linear compressibility, $v=\u2202u/\u2202t$ is the particle velocity, and *p* is the acoustic pressure. A relevant constraint is the constitutive equation

Using these definitions and assuming no external body forces (i.e., *F* = 0; this restraint will be relaxed below), Eq. (2) becomes

Carrying out the variation then yields

The variations of the particle velocity and pressure are assumed to be smooth functions of space. Then integrating by parts to isolate the variations of the acoustic pressure and particle velocity yields

Since variations are arbitrary the integrand of the left-hand side of the above equation yields three independent equations that must be simultaneously satisfied in the bulk medium:

The last of these equations is just the constraint *C* = 0, and the first two may be used to eliminate the Lagrange multiplier *μ*, yielding

which is just Newton's second law recovered.

The initial and final terms on the right-hand side of Eq. (8) may be neglected by assuming the pressure is zero at the possibly remote times *t*_{0} and *t*_{1}. The boundary term may be re-written as

As is shown below, this form of the boundary term shows that *u* and *p* must be independently continuous at interfaces. Consider an integral over $x\u2208(x0,x1)$. This integral may be arbitrarily broken into two integrals, one over $x\u2208(x0,x\u2032)$ and one over $x\u2208(x\u2032,x1)$, where $x0<x\u2032<x1$. Each of these integrals lead to boundary conditions of the same form as Eq. (11), leading to the corresponding equation

where the subscript < denotes the quantity approaching $x\u2032$ with $x<x\u2032$ and the subscript > denotes the quantity approaching $x\u2032$ with $x>x\u2032$. Since the final term on the right is identical to the boundary term without splitting the integral, we find that the remaining terms on the right must be equal to zero. This can only be the case if $p<\delta u<=p>\delta u>$. Since the variation is arbitrary, we may therefore conclude that $p<=p>$ and $\delta u<=\delta u>$, or, since the variation is also assumed to be smooth, $u<=u>$. The two integrals may be treated as separated domains that might have had independent material properties. Thus, the above result implies that the acoustic pressure and the particle displacement are continuous at arbitrary interfaces.

Hamilton's principle as described above is a global description of the system, or every part of the system must be accounted for simultaneously. However, the Lagrangian density and the constraint equation are written in terms of local variables. Thus, Hamilton's principle provides a way to transition from locally valid approximations to a globally valid approximation. Section II B presents a multiple scales-based analysis providing locally valid approximations of the acoustic variables, and Sec. II C explains the globally valid homogenization.

### B. Locally valid multiple scales approximation

Consider a static, layered, inhomogeneous, one-dimensional, lossless acoustic system with a representative volume element length of *L*. Each layer is assumed to be a homogeneous acoustic/elastic material. By defining *L* we implicitly assume that the statistics of the inhomogeneities are stationary and ergodic, or that in any sub-domain of length *L* of the system the relevant statistics are, approximately, the same. (Note that choosing *L* for a given system is a non-trivial problem; see Ref. 10 for a discussion.) Locally, the behavior of the system may be represented by the constitutive equation and the equation of conservation of momentum, which may be written (in Lagrangian coordinates), respectively, as

where

is the strain and $B$ is a coefficient of quadratic nonlinearity in the constitutive relation. The ellipsis at the end of the constitutive equation denotes that higher-order terms in the strain may also appear. To determine the relative importance of each of these terms, define the characteristic amplitudes of the pressure and velocity *p*_{0} and *v*_{0}, respectively, and the smallest characteristic time scale of the acoustic signal *T*. These characteristic quantities are defined such that $\Pi \u2261p/p0$ and $\Lambda \u2261v/v0$ are *O*(1) and derivatives of acoustic variables with respect to time are *O*(*T*). Defining $\sigma \u2261x/L$ and $\tau \u2261t/T$, Eqs. (13)–(15) may then be written as

where $\u2202\nu /\u2202\tau \u2261\Lambda $. Let the characteristic acoustic Mach number be $\u03f5\u2261v0/cmin,$ where $cmin$ is the smallest small-signal sound speed present in the system. Assuming weak nonlinearity such that terms $O(\u03f52)$ relative to the largest order term may be neglected, we may write

where $\eta \u2261L/cminT,\u2009M\u2261\beta p0cmin/v0,\u2009R\u2261\rho cminv0/p0,N\u22611+B$.

For $\eta \u226a1$ and $\eta \u22122\u03f5=O(1)$ (long-wavelength limit with relatively weak nonlinearity) a multiple-scales analysis may be used to perform an asymptotic homogenization. Define the inhomogeneity-scale position $\sigma 0\u2261\sigma $ and the *n*th macroscopic position $\sigma n\u2261\sigma /\eta n$ such that

Assume that any acoustic field $X$ (e.g., Λ or Π) may be written as a convergent power series in *η* as $X=X0+\eta X1+\u2009\eta 2X2+\cdots $. Substituting these series forms into Eqs. (19) and (20) and collecting powers of *η* yields

The coefficients of each order of *η* may be treated as independent terms leading to separate equations. The $O(\eta 0)$ equations are

and the $O(\eta 1)$ equations are

Note that the nonlinear terms do not appear in the $O(\eta 1)$ equations due to the $O(\eta 0)$ equations.

Since $\Pi 0$ and $\Lambda 0$ are independent of *σ*_{0} the $O(\eta 1)$ acoustic variables will contain secular terms (i.e., unbounded in *σ*_{0}) if the right-hand sides of Eqs. (25) are not identically equal to zero, yielding

These equations may be combined yielding the linear wave equations accurate to $O(\eta 2)$ (see Ref. 11):

Note that while the effective wave speed is a function of *σ*_{0} and is inhomogeneous, the wave equation is homogeneous in *σ*_{1}. Any solution which is smooth in *τ* (assumed to be true since the smallest characteristic time scale is *T* > 0) is also smooth in *σ*_{1}, and so a convergent Taylor series may be written as

Recalling that $\sigma 1=\sigma /\eta $ then shows that the first term is $O(\eta 0)$, the second term is $O(\eta 1)$, the third term is $O(\eta 2)$, and so on. However, the solution $X=X0$ as derived above is only valid to $O(\eta 2)$, and so we may write

Thus the acoustic variables may be considered completely independent of the microstructure to $O(\eta 1)$, and vary only linearly with respect to the position within a domain of length *L* to $O(\eta 2)$. Written in this form it becomes clear that $X0$ may be interpreted as the macroscopic approximation of $X$.

Returning to dimensional variables yields the results in the form

where $x=X+\chi $, *X* is the macroscopic position, $\chi =O(L)$ is the position relative to the macroscopic position, and *V* and *P* are the macroscopic particle velocity and acoustic pressure. See Fig. 1 for a summary schematic. These results are valid even for weakly nonlinear waves.

It is important to note that the analysis presented above assumed that the inhomogeneities only consist of layers of material with constant mass density and compressibility. Thus, a modified analysis is necessary to account for layers with dynamic material properties, additional physics (e.g., piezoelectricity), or hidden degrees of freedom (e.g., side-wall resonators).

### C. Globally valid homogenization

Hamilton's principle is global in nature due in large part to the fact that it is represented as an integral equation, and so naturally the integral itself is important to global homogenization. Consider the spatial integral over some function *f*(*x*),

Any integral may be divided into a summation of sub-integrals that span the same domain as the original integral. Let $I$ be divided into segments of length *L*, or into *M* parts such that $M=(x1\u2212x0)/L$:

See Fig. 2 for a schematic with explanations of the relevant variables. If we identify $x0+(m\u22121/2)L$ with the macroscopic position *X _{m}* then the integral may be written as

where $\chi =x\u2212Xm=O(L)$. The quantity within the square brackets is then identified as the local *χ*-average of *f* centered at *X _{m}*, and may be written as

where $\u27e8\u2009\xb7\u2009\u27e9$ is the local *χ* average, defined as

As a continuum approximation (which is itself a form of a multiple-scales approximation), the element length *L* may be approximated as being infinitesimal, and the summation may be approximated by an integral over *X*:

This approximation is valid if *L* is smaller than the smallest characteristic length scale of the macroscopic behavior, i.e., $cminT$. Thus, the continuum approximation given here is appropriate for $\eta \u226a1$. Extending this result to larger values of *η* requires special care (see Ref. 12).

Using the results from Eqs. (30) and (36), Hamilton's principle in Eq. (2) with inhomogeneous material properties may be written approximately as

where

Since the Lagrange multiplier *μ* takes whatever form necessary to enforce the constraint, it may be assumed to not be a function of *χ* and therefore need not be averaged over. The average of the constraint term, however, requires some care. Since the original constraint is just *C* = 0 for all *x*, then $\alpha C=0$ could also have been used, where *α* is any known function of *x*. However, *α* then appears in the constraint average as $\u27e8\alpha C\u27e9$, and is therefore interpreted as a weight. This observation leads to the question, what weight should be chosen? A useful tool to answer this question is the covariance, defined as

A particularly useful identity of the covariance is that $cov(a,b)\u2264var(a)var(b)$, where $var(a)=cov(a,a)$ is the variance of *a*. Thus, if the variance of *a* is zero then so is the covariance of *a* and *b*, regardless of what *b* is. With this in mind, consider the constraint of the constitutive equation. The average of this constraint then becomes

The covariances on the right represent the errors of the continuum approximation, and so it is desirable for them to be as small as possible. From the multiple-scales analysis above we know that the variance of *p* and *v* are $O(\eta 4)$ [since *p* and *v* are approximated up to $O(\eta 2)$ and $var(p)$ is $O(p2)$], and so $cov(\alpha \beta ,p)=O(\eta 2)$ as long as $var(\alpha \beta )=O(1)$. However, there is no guarantee that the variance of $\u2202v/\u2202x$ is small. In fact, throughout the RVE the velocity gradient is proportional to the compressibility times the near-constant pressure. By choosing *α* = 1 (or any other constant) the final covariance may be set identically to zero, leaving the accuracy of the continuum approximation $O(\eta 2)$. From this example a useful rule is to minimize the variance of any function multiplying the gradient of an acoustic variable.

The relevant variations for Eq. (37) may be written to $O(\eta 2)$ as

Derivatives of the variations of *P* and *V* may be moved via integration by parts (boundary terms will be considered below). Thus the integrand of Eq. (37) may be written as

Eliminating *μ* from the resulting *δV* and *δP* equations yields the homogenized dynamic equation

Together with Eq. (40) and $\u27e8C\u27e9=0$, the homogenized dynamic equation describes the effective behavior of the system to $O(\eta 2)$, and yields the wave equation

where $\rho eff$ is the effective mass density, $\beta eff$ is the effective compressibility, and $\psi eff$ is a material property that couples the traditional momentum balance equation with the traditional pressure−strain relationship, and has been called the Willis coupling.^{6,13,14} As has been noted in previous research,^{15,16} Willis coupling appears only for intrinsically asymmetric microstructures. The effective Willis coupling does not appear in the effective wave equation, as has been noted in previous research,^{17} and is only manifested at interfaces, making it interesting for applications such as impedance matching and surface waves. It should be noted that Willis coupling is only guaranteed to be absent in one-dimensional effective wave equations; no such guarantee occurs for higher dimensional systems, even if only plane waves are considered.^{18}

Now consider the boundary terms obtained by integrating Eqs. (43)–(45) by parts. Neglecting temporal boundaries, the boundary terms may be written as

where *U* is the macroscopic particle displacement. From Eq. (46) we see that $SP+\u2202\mu /\u2202t=O(\eta )$, and since $\u27e8\chi \beta \u27e9=O(\eta )$ the coefficients of *δP* are $O(\eta 2)$ and may be neglected. Therefore, the boundary terms may be written as

The term $\u2202\mu /\u2202t$ may be substituted to $O(\eta 2)$ with $\u2212P+(\u27e8\chi \beta \u27e9/\u27e8\beta \u27e9)\u2202V/\u2202t$. Following the same logic as found around Eq. (12) we therefore conclude that the macroscopic displacement and the pressure-related quantity

are continuous at interfaces to $O(\eta 2)$.

It is interesting to note that the macroscopic continuity and dynamic equations may be written in terms of $P\u02c7,$ respectively, as

These equations lead to a slightly different set of effective material properties, the difference being that the effective Willis coupling $\psi eff=\u27e8\rho \u27e9\u27e8\chi \beta \u27e9/\u27e8\beta \u27e9$ is replaced by $\psi eff=\u27e8\chi \rho \u27e9$.

The homogenization results derived above are similar to the results of many other methods. However, the present method has the benefit of versatility beyond other methods, as shown in Sec. III.

## III. EXAMPLES

For the sake of simplicity, all analyses presented below will neglect terms $O(\eta )$. As a consequence of this simplification, Willis coupling is also neglected.

### A. Viscous media

The homogenization results derived above neglected body forces, but they may be readily accounted for. Consider the case of a one-dimensional heterogeneous acoustic system with weak viscous drag, such that the body force in Eq. (2) may be expressed as^{19}

where $\mu visc$ is the coefficient of viscosity. Note that this viscous term is actually just the first term of a constitutive equation, and so is only valid up to $O(\mu visc)$. The locally valid approximation given above may still be used as long as the viscosity term is negligible over the length scale *L*, which is the case for $\eta \u226b(\mu visc/\rho cL)max,$ where the subscript “max” denotes the largest value in the system. In this case, Eq. (37) becomes

where $\u27e8Ekin\u27e9,\u2009\u27e8Epot\u27e9$, and $\u27e8C\u27e9$ are given as

The homogenization of the left-hand side proceeds as above. The right-hand side is written in terms of gradients of the velocity and the variation of the displacement. Using the local dynamic and constraint equations and integrating-by-parts in time the right-hand side of Eq. (56) may be written in the form

where no gradients of acoustic variables are present. Then the homogenization procedure given above may be utilized to write

Thus Eq. (56) becomes

which leads to the equations

with $\u2212\u2202\mu /\u2202t=P$ and *U* being continuous at interfaces. Combining these three equations leads to the effective dynamic equation

Since this equation is only valid up to $O(\mu visc)$, the quantity $\u22022V/\u2202t2$ on the right may be substituted with $(1/\u27e8\beta \u27e9\u27e8\rho \u27e9)\u22022V/\u2202X2$ without changing the order of accuracy, leading to

where the gradient of the compressibility is assumed to be negligibly small. Thus, the effective homogeneous coefficient of viscosity may be defined as

### B. Finite-amplitude systems

The homogenization methods described above are valid for weakly nonlinear systems, such that $O(\u03f5)$ terms may be kept without changing the homogenization procedure. This example considers a heterogeneous one-dimensional system without losses that has amplitudes sufficiently high that first-order nonlinear terms must be maintained for accuracy.

The form of the potential energy defined in Eq. (4) is only valid for linear propagation, and is a special case of a more general expression,

with the pressure and strain being related through Eq. (13), which may be inverted as

Since the pressure is continuous at interfaces and the strain is not, the most accurate homogenization may be obtained by writing the potential energy entirely in terms of pressure and using Eq. (70) as a constraint equation. In this case and neglecting body forces, Hamilton's principle yields

Since the strain is isolated from anything except the Lagrange multiplier, this form of Hamilton's principle is appropriate for homogenization, which yields

where

is the macroscopic strain. In this case further analysis is performed more naturally in terms of the macroscopic particle displacement *U* rather than the velocity *V*, and so we write Hamilton's principle as

Taking the variation, integrating by parts, collecting like variations and setting the associated coefficients to zero then yields the equations

and $\mu (1+\u2202U/\u2202X)\delta U$ is continuous at interfaces. These equation may be combined to yield the nonlinear wave equation for the macroscopic displacement:

The coefficient of the nonlinear term on the right may be identified as the effective coefficient of nonlinearity $\beta nl$.^{20} For homogeneous systems, $B$ may be written in terms of $\beta nl$ as

Then, the effective coefficient of nonlinearity may be written as

It is interesting to compare this result with the effective coefficient of nonlinearity obtained by Everbach *et al.*,^{21} which is just $\beta nl,eff=\u27e8\beta nl\beta 2\u27e9/\u27e8\beta \u27e92$. Note that in the limit of small variance of the compressibility these two predictions coincide. Another note is that while the present analysis is cast in Lagrangian coordinates, the effective material properties do not rely on the coordinate system and are applicable for Eulerian analyses as well.

Incorporating nonlinearity in propagation usually leads to changes in the spectral amplitudes of a signal. In particular, quadratic nonlinearity without sufficient absorption leads to shock formation. Since shocks are characterized by very rapid changes in the acoustic variables, the validity of the multiple scales approximation that is used in the present homogenization scheme comes into question. Signals propagating without microscopic or macroscopic losses will eventually steepen to the point that $\eta \u226a/1$ and treating the system as an effective medium will become inappropriate. However, for propagation distances much shorter than the macroscopic shock formation distance it may be reasonable to use the effective properties given here.

As a demonstration of the homogenization method, consider a system consisting of 90% water and 10% alcohol placed in periodic layers of length *L* = 10 cm (see Fig. 3). The constituent materials were chosen due to the ready availability of their traditional and nonlinear material properties. The mass density, compressibility, sound speed, and coefficient of nonlinearity for the water, alcohol, and the effective medium are summarized in Table I. It is interesting to note that the effective mass density and sound speed are very similar to that of water, while the compressibility and, especially, the coefficient of nonlinearity are significantly larger.

Property . | Symbol . | Water . | Alcohol . | Effective . |
---|---|---|---|---|

Mass density (kg/m^{3}) | ρ | 998 | 790 | 977 |

Sound speed (m/s) | c | 1480 | 1150 | 1421 |

Compressibility ($10\u221210$ Pa^{−1}) | β | 4.6 | 9.6 | 5.1 |

Coefficient of nonlinearity | $\beta nl$ | 3.5 | 6.3 | 4.9 |

Property . | Symbol . | Water . | Alcohol . | Effective . |
---|---|---|---|---|

Mass density (kg/m^{3}) | ρ | 998 | 790 | 977 |

Sound speed (m/s) | c | 1480 | 1150 | 1421 |

Compressibility ($10\u221210$ Pa^{−1}) | β | 4.6 | 9.6 | 5.1 |

Coefficient of nonlinearity | $\beta nl$ | 3.5 | 6.3 | 4.9 |

COMSOL Multiphysics 5.4 was used to simulate a signal propagating through the above described system. The source signal is a Gaussian pulse given by

with $p0=29.11$ MPa, $t0=5L/ceff=0.35$ ms, and $ceff=1421$ m/s is the effective sound speed. The amplitude was chosen such that the effective shock formation distance would be at 8 m. It should be noted that using the Everbach prediction of the effective coefficient of nonlinearity with the above waveform and system would lead to an estimate of the shock formation distance being 8.22 m. The nonlinear propagation was modeled by the second-order Westervelt equation without losses.^{24} The maximum element size was chosen to be $L/40$ and the maximum frequency to resolve was chosen to be $10/t0=28.4$ kHz.

The COMSOL simulation using the full system and the effective system are shown as a function of distance at several times in Fig. 3. Immediately evident is the nonlinear steepening of the propagating waves, and by 5 ms a shock wave has almost completely formed around the shock formation distance of 8 m. The signals for the full system (solid lines) and the effective system (dashed lines) are very similar for all of the times shown. Of particular note, the shock fronts of the full and effective predictions occur at the same locations for a given time. However, at 4 and 5 ms the full system predictions begin to show small scale high-frequency oscillations that do not appear in the effective system predictions, especially for positions past 5 m. These oscillations are likely due to the fact that shorter wavelengths become important as the waves steepen, making the effective medium approximation become less applicable. This argument may be explored analytically, as shown below.

Since the present analysis is valid to $O(\eta )$, the quantity $\eta =L/cminT$ is the important quantity to analyze. The time scale *T* depends on the time-derivative of the acoustic pressure, which was shown by Muhlestein and Gee to evolve as^{25}

By definition $\u2202p/\u2202t=O(p0/T)$, and so Eq. (82) dictates that

In the example above, $T(0)=t0$, and so at the source $\eta =(ceff/cmin)/5=0.25$. The approximation is only valid for $\eta <1$; this threshold is met at *x* = 5.17 m, which correlates well with the appearance of the high-frequency oscillations that appear in Fig. 3.

While further quantitative analyses of this example could be given, the qualitative discussion given here is sufficient to demonstrate the efficacy of the present homogenization method.

## IV. SUMMARY

A one-dimensional acoustic metamaterial homogenization method has been derived and demonstrated. A formal multiple-scales method was used to obtain a locally valid approximations of the acoustic variables. These approximate variables were then substituted into Hamilton's principle, which is global in nature, and then averaged to obtain estimates of the effective material properties of heterogeneous systems. While this homogenization method is approximate in nature, it is also highly versatile. This versatility was demonstrated with two examples: accounting for viscosity, and accounting for finite-amplitude effects.

This method may be considered a framework for additional development. Three possible improvements are mentioned here. First, the present analysis is restricted to one dimension, but there is no fundamental reason that it may not be extended to two and three dimensions. Second, hidden degrees of freedom, such as side-wall resonators for plane wave tubes and mass-in-mass systems, may also be incorporated. Third, additional physics, such as thermal and piezoelectric phenomena, may be accounted for using the method of Lagrange multipliers.

## ACKNOWLEDGMENTS

The author would like to thank Carl Hart and Kyle Dunn for excellent discussions of the manuscript at many stages, and the anonymous reviewers for their helpful suggestions and guidance. Permission to publish was granted by the Director, Cold Regions Research and Engineering Laboratory.