## Abstract

Many materials under multiaxial periodic loading exhibit rate-independent internal dissipation per cycle. Constitutive modelling for such dissipation under spatially variable triaxial stresses is needed for calculating modal damping of solid bodies using computational packages. Towards a micromechanically motivated model for such dissipation, this paper begins with a frictional microcrack in a linearly elastic solid under far-field time-periodic tractions. The material is assumed to contain many such non-interacting microcracks. Single-crack simulations, in two and three dimensions, are conducted using ABAQUS. The net cyclic single-crack dissipation under arbitrary triaxial stresses is found to match, up to one fitted constant, a formula based on a pseudostatic spring-block model. That formula is used to average the energy dissipation from many randomly oriented microcracks using Monte Carlo averaging. A multivariate polynomial is fitted to the Monte Carlo results. The polynomial is used in finite-element simulation of a solid object, wherein modal analysis is followed by computation of the net cyclic energy dissipation via elementwise integration. The net dissipation yields an equivalent modal damping. In summary, starting from a known formula for a single crack, this paper develops and implements a method for computationally modelling the modal damping of arbitrarily shaped solid bodies.

## 1. Introduction

This work is motivated by an engineering design problem. Consider two possible designs (shapes and dimensions) for some engineering component, to be made of some known, lightly dissipative material. Which design has better vibration damping?

A computational approach would begin with straightforward finite-element (FE)-based modal analysis. It would also require a constitutive relation for material damping under time-periodic triaxial inhomogeneous stresses. As discussed below, such a relation is not presently available. To fix ideas, consider the solid object of figure 1*a*, modelled using an FE package, with its first mode as shown in figure 1*b*. We seek a constitutive relation that can be used, along with modal analysis, to compute the damping ratios for the first several modes of such an object.

We begin by reviewing the relevant literature on internal damping. For many solids vibrating at low frequencies, it has been observed that internal energy dissipation (or material damping, or internal damping) per unit volume and per cycle of deformation is *frequency independent* and approximately proportional to some power (*n*≥2) of the stress amplitude. We write
1.1where *D*_{m} stands for specific material damping, *σ*_{eq} is a suitable stress amplitude, and *J* and *n* are fitted constants. Such frequency independence was reported by Lord Kelvin [1]. Rowett [2] studied torsion of steel tubes and observed frequency-independent power-law dissipation with *n*≈3. Kimball & Lovell [3] found *n*≈2 and frequency independence for 18 different materials, including metals, celluloid, glass, rubber and wood. Such frequency-independent power-law behaviour has been discussed by many others. Lazan [4] adopted a phenomenological approach and used power laws. Granato & Lücke [5] proposed an explanation based on dislocation pinning by impurity particles. Dawson [6] considered an unknown function of non-dimensionalized stress, formally expanded in a Taylor series using even powers only, leading by assumption to *n*=2 for small stresses.

However, as indicated above, we are interested in multiaxial stress states. No convincing mechanically based engineering model for the same is presently available. For example, the empirical laws in Lazan [4] do not identify the equivalent stress of equation (1.1) under an arbitrary triaxial load. Dislocation-based models as in Granato & Lücke [5] involve several parameters related to the crystal structure, yet to be translated into measurable external macroscopic model parameters. The approach of Dawson [6] also does not identify the role of multiaxial stresses in the dissipation model. As a final example, Hooker [7] proposed that the equivalent stress amplitude should be computed as
1.2where *I*_{1} and *I*_{2} are the first and second stress invariants, respectively, and *λ* is a fitted parameter. The above is motivated by the fact that it is a linear combination of distorsional and dilatational strain energies.

In contrast to the above, there is in fact a micromechanical modelling approach that seems promising for our purposes. This approach considers randomly distributed frictional microcracks within an elastic material, as opposed to the *ad hoc* prescription of equation (1.2).

The literature on frictional microcracks in elastic materials is rich. Kachanov [8] proposed an approximate analysis method based on a superposition principle for interactions of multiple cracks at moderate distances from each other. Aleshin & Abeele [9] presented a tensorial stress–strain hysteresis model due to friction in unconforming grain contacts. Their model pays close attention to the variation of actual area of contact under normal stress, but has many parameters (seven for uniaxial compression alone). Deshpande & Evans [10] studied frictional microcracks in the context of inelastic deformation and fracture of ceramics. Al-Rub & Palazotto [11] computed energy dissipation in ceramic coatings and found that frictional dissipation in microcracks contributes significantly to overall dissipation. More recently, Barber and co-workers have presented several papers on mechanics with frictional microcracks. Jang & Barber [12] discussed the dissipation in interacting microcracks using Kachanov's [8] approach. Barber *et al.* [13] studied a single frictional elastic contact subjected to periodic loading. The contact has an extended area, part of which sticks while the remainder can slip; some basic results about energy dissipation under periodic loading are obtained. Jang & Barber [14] examined the substantial effect of the relative phase of harmonically varying tangential and normal loads on the dissipation in an uncoupled frictional system. Barber [15] discussed discrete frictional systems under oscillating loads, and examined the conditions under which the steady-state solution retains a memory of the initial state. Individual cracks with surface roughness models for the contacting faces have been studied by Putignano *et al.* [16], who observed that, with microslip in variable regions but no gross slip at the contact surfaces, energy dissipation varies as the cube of the stress amplitude for small amplitudes.

None of the above studies have considered the net dissipation in a body, with spatially variable stresses, from a multitude of randomly oriented frictional microcracks. Here, we seek a macroscopic constitutive relation for such dissipation.

To this end, for simplicity, we assume that the cracks are small and far apart (non-interacting); that initial super-small microslips on crack surface asperities can be neglected, and the crack faces modelled as non-adhering yet geometrically flat; that all crack face frictions can be modelled using a single coulomb friction coefficient *μ*; and that the material remains linearly elastic. We assume that the normal and tangential loadings at each crack are in phase (as would occur for vibration of a body in a single mode). Although our initial formulation allows both non-zero mean stresses and non-uniform distributions of crack face orientations, our final formula assumes zero mean stress and uniformly distributed crack face orientations.

Under these assumptions, we develop an empirical formula with two fitted parameters for the dissipation per unit volume and per cycle of time-periodic triaxial stress, and use the solid object of figure 1*a* in a computational example.

The contribution of this paper may thus be viewed as the first assembly of the following tasks in one self-contained sequence: we use existing ideas about bodies with frictional microcracks, integrate their dissipation rate over all possible crack orientations, develop a constitutive model for the net specific dissipation per cycle, and demonstrate the use of the formula to compute the modal damping ratios of arbitrarily shaped objects.

## 2. Frictional dissipation in a single microcrack

The first step in our task is to compute the frictional dissipation in a single microcrack due to remote cyclic loading. For the two-dimensional case, an analytical treatment is given in Jang & Barber [12], including an analytical formula based on a single degree of freedom spring-block model for when the entire crack face sticks or slips as one. We have studied the same using FE calculations in both two and three dimensions. The formula based on the spring-block model, with one fitted constant, turns out to be highly accurate. These computations are outlined below (for details, see the electronic supplementary material).

Look at figure 2*a*, where the crack is circular, planar and parallel to the horizontal faces of a cube-shaped element. We have found in separate two-dimensional computations that removing tiny regions around the crack tips does not influence the net dissipation, and so material or geometric nonlinearities are neglected.

We can simplify the three-dimensional picture. To any given state of stress and deformation, an additional *σ*_{x} causes an additional uniform strain in the element but no added shearing at the crack face. Consequently, *σ*_{x} causes no slip and does not contribute to the energy dissipated. Similar arguments apply for *σ*_{y} and *τ*_{xy}. Thus *σ*_{x}, *σ*_{y} and *τ*_{xy} do not affect our results and are dropped. By contrast, *σ*_{z} will affect the frictional forces, while *τ*_{zx} and *τ*_{zy} will cause frictional sliding. Finally, by rotating the coordinate system about the *z*-axis, *τ*_{zx} can be made zero. This leaves just the normal stress *σ*_{z}=*σ* and the shear component *τ*_{zy}=*τ* contributing to frictional dissipation. A two-dimensional representation is shown in figure 2*b*.

Figure 2*c* depicts possible periodic normal and shear tractions, acting in phase, on the element of interest. Inertia is negligible (alternatively, frequencies are low), so the dissipation is frequency-independent. Our pseudostatic simulation results will actually apply to every periodic waveform where (i) the *changes* in far-field normal and shear stresses maintain a fixed proportion throughout the load cycle (i.e. the time-varying parts have similar waveforms) and (ii) there are only two points of stress reversal per cycle. See the electronic supplementary material for some possible waveforms. For convenience, we used a triangular loading pattern in our computations described below.

### (a) Finite-element simulations

Although we are interested in small flat cracks in three dimensions, we began with detailed two-dimensional (plane stress^{1}) simulations in ABAQUS. Subsequently, we carried out three-dimensional calculations with both a circular crack as well as a symmetrically loaded isosceles triangular crack. Details are presented in the electronic supplementary material. In particular, convergence was verified for both mesh size in space and load steps in time. Several load cases were run in each case with different values of stress amplitudes as well as the friction coefficient *μ*. Results for five cases in the two-dimensional simulation are given in table 1 (many more such numerical results, for both two- and three-dimensional cases, are listed in the electronic supplementary material).

We acknowledge that ABAQUS gave us *one* solution in each case, without proving uniqueness. Jang & Barber [14] have discussed uniqueness in a more complicated, but two-dimensional, setting. Our three-dimensional computations using ABAQUS had the specific goal of developing a physically defensible constitutive relation for damping. We neither sought nor noticed evidence of dynamic waves near the sliding crack face.^{2}

### (b) Dissipation formula using a spring-block system

We now present an analytical formula that captures every result of the kind exemplified in table 1. The formula is not new: it is given, in a different form, in Jang & Barber [12]. However, we include it below because it plays a key role in this paper. For details, see the electronic supplementary material.

We consider a spring and massless block system. The block slides on a frictional surface (coefficient *μ*), and is attached to a rigid wall through the spring. Periodic normal (*σ*) and tangential loads (*τ*) act on the block. When *σ*<0, there is no friction. To the extent that the entire crack face slips or sticks as one, this single degree of freedom model should be accurate: we will find below that it is.

Define
Note that *τ*_{a} and *σ*_{a}, being amplitudes, are positive by definition. The mean normal stress *σ*_{m} is taken positive when compressive (figure 2*c*). The mean shear stress *τ*_{m} affects the mean position but not the steady-state cyclic dissipation (see the first two rows of table 1, as well as the electronic supplementary material). The dissipation per cycle in this system can be shown to be
2.1which includes a single load- and friction-independent fitted constant *C*. The square brackets denote logical variables (equal to 1 if the inequality holds and 0 otherwise). All results for a given crack should fit this formula (as they do, below).

We observe that *D* in equation (2.1) depends only on two-dimensional quantities: the fitted constant *C* and the normal stress amplitude *σ*_{a}. Consequently, all other non-dimensional ratios held constant, *D* varies as the square of the stress.

Dissipation results from the two-dimensional FE simulations are plotted against the dissipation predicted by equation (2.1), with one fitted constant, in figure 3. The match is excellent. Similar matches, with a different *C* in each case, were obtained for the circular and triangular cracks in three dimensions (see the electronic supplementary material). Thus, equation (2.1) is acceptable for our purposes.

We now turn to the use of equation (2.1) and a Monte Carlo method to compute the *average* dissipation from a multitude of randomly dispersed and oriented microcracks within the material, assuming their interactions may be neglected.

## 3. Dissipation due to multiple cracks: Monte Carlo method

### (a) Random orientations

We now consider randomly oriented cracks, and for simplicity assume that all orientations are equally likely, i.e. the material is macroscopically isotropic. Geometrically, the normals () to the crack faces are uniformly distributed on the surface of the unit hemisphere (figure 4). These points were generated by first generating points uniformly distributed within the upper half of a cube, then discarding points that lay outside an appropriate sphere, and finally by projecting the remaining points radially outward onto the surface of the sphere.

### (b) Non-dimensionalization

Let the stress state of interest be . All orientations of the crack faces being equally likely, the coordinate system is irrelevant. Crack size and shape affect constant *C*, but we assume *C* is independent of and can be averaged separately. Here we take *C*=1; we can multiply by a fitted constant later.

Since the coordinate system used to describe **S** is irrelevant, it is simplest to think in terms of principal stresses (*σ*_{1}≥*σ*_{2}≥*σ*_{3}). For *σ*_{1}=*σ*_{3}, there is no shear stress, *ζ*<*μ*, and the dissipation is zero (see equation (2.1)). Accordingly, we assume *σ*_{1}>*σ*_{3} and first scale the stress so that *σ*_{1}−*σ*_{3}=1. Later, we will multiply back by the square of the scaling factor (recall the discussion following equation (2.1)).

We can visualize the time-harmonic part of the stress state (i.e. matrix **S**) and associated dissipation possibilities using Mohr's circles for three-dimensional stresses (figure 5*a*). For any given normal , the resultant shear (*τ*, assumed positive) and normal stress (*σ*), represented as a point (*σ*,*τ*) on the Mohr diagram, will lie in a region bounded by three circles [19]. The dissipation corresponding to any such point will be zero unless the point lies outside the friction wedge, corresponding to *ζ*>*μ* in equation (2.1), as indicated in figure 5*a*.

As indicated in figure 5*a*, stress state B can be reflected to stress state A, because **S** is multiplied by in any case. Accordingly, we can assume that the centre of the largest Mohr circle is on the non-negative real axis (*σ*_{1}+*σ*_{3}≥0).

Figure 5 also shows that, for any *μ*>0 and *σ*_{1}−*σ*_{3}=1, for *σ*_{1} sufficiently large, there is no dissipation. Conversely, for a state of simple shear, with *σ*_{1}=0.5, *σ*_{2}=0 and *σ*_{3}=−0.5, there is non-zero dissipation for any *μ*>0. In other words (figure 5*b*), the hydrostatic part of **S** affects the dissipation per cycle.

For clarity, we write down the sign change and scaling described above using a single formula. If the actual time-harmonic state of stress is with eigenvalues , then we use the scaled stress
3.1where the square brackets denote a logical variable as before. The principal stresses corresponding to **S** are denoted by *σ*_{1}≥*σ*_{2}≥*σ*_{3}. It is now assured that *σ*_{1}+*σ*_{3}≥0 and that *σ*_{1}−*σ*_{3}=1. The scaling factor is
3.2and the dissipation obtained using **S** will be divided by to obtain the dissipation due to . Later, for spatially varying stresses in a body vibrating in a given mode, we will use equation (3.1) repeatedly using a computer program.

At this point, we have non-dimensionalized the stress amplitude **S** by taking *C*=1 and using *k*_{f} as above. The non-dimensional **S** can now be specified by *σ*_{1}≥0.5 along with *σ*_{2}≥*σ*_{1}−1. These two non-dimensional principal stresses, along with *μ*>0, are inputs to the dissipation calculation.

### (c) Building block: a single average dissipation calculation

We have so far taken three steps towards our goal: (i) we have verified for three dimensions a formula for the dissipation at a single frictional microcrack due to arbitrary oscillating imposed far-field stresses, (ii) we have generated a large number of random normal directions , and (iii) we have non-dimensionalized and scaled the stress ‘amplitude’ matrix .

Now, for our fourth step, we find the average dissipation for a given scaled stress matrix **S**, mean stress **S**_{0} (non-time-varying and identically scaled), and friction *μ*, by sequentially considering the full random set of normal directions. For the results below, we have used 4.5 million such directions.

For each normal direction , we find the time-varying part of the traction vector . The normal stress amplitude *σ*_{a} then is
If a mean stress matrix **S**_{0} is given, the corresponding mean normal stress *σ*_{m} is similarly calculated using the same . Subtracting from **t** gives the shear component of the traction; and its magnitude is *τ*_{a}. The shear component is not computed from the mean stress **S**_{0} because it plays no role in cyclic dissipation, as mentioned above (§2*b*). Now, for the given , the dissipation is computed using the formula of equation (2.1) with the fitted constant *C* set equal to unity as explained earlier. The above dissipation is finally averaged over all the to account for the random orientation of the crack.

For simplicity, the non-zero mean stress **S**_{0} is not retained in what follows.

For initial demonstration, we consider uniaxial tension and simple shear. The principal stresses for tension are *σ*_{1}=1, *σ*_{2}=0 and *σ*_{3}=0. For shear, they are *σ*_{1}=0.5, *σ*_{2}=0 and *σ*_{3}=−0.5. These stresses are already normalized: *σ*_{1}−*σ*_{3}=1 in both cases. Dissipations for these cases are calculated for various *μ* values (0.01–1.2) and shown in figure 6. The ratio of the two varies significantly with *μ*. Note that we view *μ* as a fitted parameter in our two-parameter constitutive model (the other parameter is the overall multiplicative factor of *C*).

As an example of how *μ* might be fitted consider Robertson & Yorgiadis [20], who sought the same specific dissipation per cycle in two different loading conditions.^{3} For such equal dissipation to occur under both pure (simple) shear and pure extension, the ratio of the shear stress amplitude during torsional vibration to the normal stress amplitude during longitudinal vibration was found to be between 0.48 and 0.60 (for several materials). This statement of equivalence translates, for our model, into roughly 0.6<*μ*<1.1 by the following reasoning. Fix the longitudinal stress state at *σ*_{1}=1, *σ*_{2}=0 and *σ*_{3}=0, as above. Let the shearing stress state be *σ*_{1}=0.5, *σ*_{2}=0 and *σ*_{3}=−0.5 but only after normalization by some factor *k*; in other words, the applied shear stress would be of amplitude *k*/2. Since the shear stress amplitudes found by Robertson & Yorgiadis [20] are between 0.48 and 0.6, *k* lies between 0.96 and 1.2. If *k*=0.96, the dissipation in shear will be 0.96^{2}≈0.92 times the value in figure 6*a*. Alternatively, the ratio plotted in figure 6*b*, upon division by 0.92, should give 1, implying *μ*≈1.1. Similarly, if *k*=1.2, we find *μ*≈0.6.

## 4. Fitted formula

So far, the dissipation has been computed as a function of **S**, possibly a non-zero **S**_{0}, and *μ*, using a time-consuming Monte Carlo simulation.

However, we eventually want to compute the modal damping of a given object of arbitrary shape. For each mode, the stress state varies spatially. We cannot do Monte Carlo simulations for every point on the body. So our fifth step is to summarize the dissipation values obtained from Monte Carlo simulations, for the special case of **S**_{0}=**0**, using a quick multivariate polynomial fit.

### (a) A comment on prior efforts

To motivate our multivariate fitted formula, we first note some prior attempts at *ad hoc* modelling of material dissipation under multiaxial stress states. Recall equation (1.1), wherein a suitable equivalent stress amplitude needs to be defined. Damping under biaxial stresses has been studied by several authors, including Robertson & Yorgiadis [20], Whittier [21], Torvik *et al.* [22] and Mentel & Chi [23]. A review of these articles is given in the electronic supplementary material. All these authors considered at least one *ad hoc* definition equivalent to equation (1.2), possibly rearranged or differently normalized. But if equation (1.2) had general validity it would apply to our dissipation results as well, since these are derived from legitimate (though approximated) physics. To check the same, we can rewrite equation (1.2) as
4.1where the *λ*'s are fitted coefficients, and the assumed roles of the distortional and dilatational strain energies are clearly visible. In checking equation (4.1) against our dissipation results, we note that different *μ* represent different material behaviours, and so we should work with one *μ* at a time. Figure 7*a*,*b* shows least-squares fitted comparisons for two *μ* values. The poor match indicates both the inapplicability of equation (4.1) in general cases and motivates our multivariate polynomial fit below.

### (b) Inputs to the fitted multivariate polynomial formula

The scaled stress **S** and *μ* are inputs for our dissipation calculation. We will later divide the computed dissipation by (see equations (3.1) and (3.2)) to obtain the dissipation for the actual stress . We now introduce two new scaled variables.

First, define *χ*=*σ*_{1}−*σ*_{2}. The inequality *σ*_{1}−1≤*σ*_{2}≤*σ*_{1} becomes 0≤*χ*≤1. Figure 5*b* shows the non-zero-dissipation range of *σ*_{1} as 0.5≤*σ*_{1}≤*σ*_{1}_{cr}, where
4.2

For *σ*_{1}≥0.5 and 0≤*χ*≤1, with *μ* as a parameter, we now generate surface plots of the dissipation as computed from Monte Carlo simulations. We compute 11 such surface plots, for equally spaced *μ* values from 0.2 to 1.2. Two representative plots are shown in figure 8; another nine are given in the electronic supplementary material. The figure confirms that the dissipation does become zero for each *μ* when *σ*_{1} crosses *σ*_{1}_{cr} (equation (4.2)). We now introduce a final scaled variable
4.3such that there is non-zero dissipation only for *s*>0. We will now seek a single fitted formula for all these dissipation surfaces, given *s*, *χ* and *μ*.

### (c) Multivariate polynomial fit

Recall equation (2.1). Now setting *C*=1 (different *C* will be incorporated later), and for *β*=0 (zero mean stress), we write
4.4where the average is over all possible orientations of the crack face. We propose, for simplicity, a polynomial form for the fit:
4.5where the 's collectively denote 100 fitted coefficients. Our numerical fit will be much faster than Monte Carlo simulation, and given below in an easily portable form. In particular, *D*_{f} will be written as a product of three matrices, **ABM**, and the constitutive relation for damping will be
4.6Here *D* (as before) is the dissipation per unit volume and per stress cycle due to scaled stress **S**, *C* is a fitted scalar coefficient, [*s*>0] is a logical variable that ensures zero dissipation for *s*≤0 (see equation (4.3)), **M** is a column vector containing powers of *χ* as described below, **A** is a row vector containing products of powers of *s* and *μ* as described below and the matrix **B** contains fitted numerical coefficients.

We write **M** and **A** as follows:
4.7with **A**_{0}=[*s* *s*^{2} *s*^{3} *s*^{4} *s*^{5}] and **A**_{m}=*μ*^{m}**A**_{0}.

The fitted 20×5 matrix **B**, determined from a least-squares calculation, is given in appendix A. In the fit, the maximum absolute error as a percentage of the maximum for each corresponding *μ* is within 5.2 per cent, with typical errors being substantially smaller.

Figure 9*a* shows the quality of the fit using two surface plots for *μ*=0.4. One surface is from Monte Carlo simulation (accurate) and the other is from the fitted polynomial. The match is good. Another comparison is shown in figure 9*b*, where all the data points used in the fit are plotted, fitted value against original Monte Carlo value, along with a 45^{°} line. Within this plot are represented 11 equally spaced *μ* values from 0.2 to 1.2. The match is reasonably good, and can be improved if desired by using higher order polynomials or other fitting methods. But dissipation, similar to other non-ideal material behaviours involving friction, fracture and plasticity, is difficult to model accurately in any case; so we arbitrarily chose to limit the number of fitting coefficients to 100, arranged in a matrix that is easy to cut and paste.

An illustration of our dissipation calculation is now given for completeness. Let the state of stress be
an arbitrary choice. The principal stresses are . The scaling factor *k*_{f}=0.089 and the normalized eigenvalues work out to *σ*_{1}=0.964,*σ*_{2}=0.140,*σ*_{3}=−0.036, giving *χ*=0.824. We take *μ*=0.5. Then *σ*_{1}_{cr}=1.618, giving *s*=0.585. Now using equation (4.6) with *C*=1, we find *D*=0.014. Finally, the dissipation is *D*/(*k*_{f})^{2}=1.818 in appropriate units.

Two further aspects of the fit are mentioned here.

First, the polynomial fit occasionally predicts some small negative values, as suggested by the bottom left portion of figure 9*b*. We simply replace those negative predictions with zero, with negligible consequence because such values are both infrequently encountered and small.

Secondly, the dissipation computation using our fitted formula is very fast compared with the Monte Carlo Method. A total of 10 000 evaluations of the formula took about 0.2 s on an unremarkable desktop computer, where a single Monte Carlo evaluation with 4.5 million points (found separately to be large enough for the accuracy needed) took about 3 min.

## 5. Finite-element computation of modal damping

We have now completed all the steps needed to consider the effective modal damping of any given mode of an arbitrarily shaped object. We will illustrate our calculations using the solid body shown in figure 1*a*.

For completeness, we have included a brief introduction to relevant aspects of vibration theory in the electronic supplementary material. The key ideas are summarized as follows. The damping mechanism we have considered is nonlinear, but the damping is assumed to be small. For small damping, the damping plays no significant role except near resonance. Near each distinct resonant frequency (or natural frequency) of an arbitrary body, an effective damping ratio for the corresponding mode can be defined. This section is concerned with the computation of such effective *modal* damping values.

### (a) Effective damping ratio (*ζ*_{eff})

We begin with an elementary formula. A lightly damped harmonic oscillator of the form
5.1has damping ratio *ζ* which, to first order, is equivalent to
5.2where is the energy dissipated per oscillation, is the total energy of the system averaged over one cycle and the bar is to distinguish the energy from Young's modulus which will be discussed later. Equation (5.2) works only for lightly damped systems because is computed over a cycle by assuming a harmonic solution. However, it is general: it does not need a linear viscous model. The general unforced solution is then approximated as
For small *ζ*, we may often just write

For a simple analytical example, consider
5.3Assume first an approximate solution (neglecting the damping over one cycle) of
from the maximum potential or kinetic energy. The energy dissipation per cycle is, by this approximation,
5.4whence by equation (5.2) we have
A numerical solution of the nonlinear equation (5.3) with *c*=0.1 and initial conditions *x*(0)=1 and is shown in figure 10. A plot of e^{−ct/π} is also given for comparison and is seen to match the oscillation envelope very well.

The above example shows the utility of equation (5.2), which we will use below. The energy dissipation calculation below will use our fit of equation (4.6), suitably integrated over the entire vibrating object.

### (b) Finite-element prediction of effective damping ratio

Our FE computation of mode shapes and modal damping proceeds as follows. We have used the FE package ANSYS for modal analysis and related computations in this paper. Elementwise integrals below will use Gaussian quadrature [24].

The vibrating object of interest is first meshed using 10 noded tetrahedral elements (SOLID187) using automatic meshing within ANSYS. Modal analysis in ANSYS yields natural frequencies and mass-normalized mode shapes. Nodal displacements for each mode of interest are extracted using a small external program.

Using the element shape functions and the nodal displacements, the displacement field is computed and then differentiated to obtain strains, and thence stresses, at four Gauss points per element. At each Gauss point, the stress is used in conjunction with equation (4.6) to compute the dissipation per unit volume and per cycle. The dissipation in the element is then obtained by the usual weighted sum of its values at the Gauss points; and the same is added up for all elements to obtain the total energy dissipated per cycle in the vibrating object. This dissipation is .

is simply *ω*^{2}/2 because the mode shape is mass normalized.

Now the effective damping ratio is obtained using equation (5.2). Some further details are provided in the electronic supplementary material.

### (c) Results for an arbitrary solid object

We finally consider the solid object shown in figure 1*a*. It is not special, and merely represents an object that is difficult or impractical to treat analytically. The object is an unconstrained thick circular plate of radius 1 m and thickness 0.1 m, with a square hole. The edges of the hole are 0.4 m, and its centre is 0.5 m from the centre of the circle. Young's modulus (*E*), Poisson's ratio (*ν*) and density (*ρ*) of the material are arbitrarily taken as 100 GPa, 0.28 and 4000 kg m^{−3}, respectively.

A total of 32 418 elements were used for meshing the object. The effective damping was computed for the first three vibration modes. The first mode was shown in figure 1*b*. The second and third modes are shown in figure 11.

Noting that *C* has units of Pa^{−1}, we have arbitrarily chosen *C*=2*π*/*E*, where the dependence on *E* is motivated by the units and the 2*π* ensures that, for axial vibrations of a uniform rod, the effective damping is exactly equal to *D*_{f} in equation (4.5). In real applications, where *C* would be fitted from test data, the 2*π* would be replaced with a fitted constant.

We pause for a moment to take stock. Recall that our dissipation model has two fitted parameters, namely *C* and *μ*. Where *C* is an overall measure of dissipation and has been arbitrarily assumed to be *C*=2*π*/*E* above. The parameter *μ* governs the relative importance of different components of the stress in the material, as discussed earlier (recall figure 6). Here, we present simulation results for a range of *μ* values. Note that the absolute magnitude of the damping ratio for any mode depends on both *C* and *μ*, but the *relative* magnitudes (or ratios) of two modal damping ratios is a function of *μ* alone.

We now present our illustrative FE results. Values of *ζ*_{eff} against *μ* are shown for the first three vibration modes in figure 12. For the first two modes, the dependence on *μ* is near identical. However, for the third mode, it differs significantly (see the differences in deformation patterns in figures 1 and 11). Such differences might be used to estimate *μ* (recall the discussion following figure 6).

## 6. Discussion and conclusions

We began this paper with the goal of computing the modal damping ratios of a lightly damped object of arbitrary shape. We noted that, though natural frequencies and mode shapes of such objects are found routinely using modal analysis in widely available commercial FE packages, the damping ratio is not easy to compute. In fact, a study of the literature revealed that there were no satisfactory models for the per-cycle dissipation in engineering materials under arbitrary triaxial states of time-periodic stress.

In this paper, we have traced out a possible route for computing the modal damping ratios of arbitrarily shaped objects. We have first adopted a dispersed frictional microcrack-based model of dissipation in an otherwise elastic body. We have verified a single simple formula based on a sliding spring-block analogy that accurately describes the energy dissipation per cycle as obtained from both two- and three-dimensional computations. To account for the random orientations of the cracks, we have used Monte Carlo averaging for any given state of (time-periodic) stress. We have summarized the Monte Carlo simulation results for a range of the friction parameter *μ* and for arbitrary triaxial stresses using a multivariate polynomial fit. The polynomial, with three independent variables and 100 fitted coefficients, is fairly accurate and arranged for easy portability. Finally, we have chosen for demonstration a body of somewhat complex shape. Using commercial FE code, we have found its first few natural frequencies and mode shapes. For each mode, we have extracted the nodal displacements, computed stresses at Gauss points, used the fitted formula to estimate dissipation rates, integrated over the body to find the net dissipation and thereby computed the modal damping up to a fitted constant *C* and for a range of *μ*.

We anticipate that this line of work may have useful applications and extensions in future work. Clearly, such modelling provides a route to optimizing engineering component designs for damping. Additionally, such work may lead to new academic research towards incorporating residual stresses, other dissipation mechanisms, anisotropy in material properties or flaw distributions, interactions between flaws, etc. Finally, we hope that such constitutive relations might eventually be built into commercial FE codes, so that modal damping values could be computed and compared routinely along with natural frequencies and mode shapes.

## A. Appendix A. Matrix B for the constitutive model

For equation (4.6), the fitted matrix **B** was found to be

## Acknowledgements

We thank Arghya Deb, Vikranth Racherla and A. K. Mallik for discussions, and the Department of Science and Technology, Government of India, for financial support. Anonymous reviewers raised questions that helped improve the paper.

## Footnotes

↵1 Plane stress versus plane strain is equivalent if we are willing to redefine the elastic constants [17], and numerical values we use for these constants are notional in our case anyway.

↵2 An anonymous reviewer pointed us to Schallamach waves [18]. Such waves seem unlikely here.

↵3 A brief review of other multiaxial dissipation experiments is given in the electronic supplementary material.

- Received November 16, 2012.
- Accepted January 8, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.