## Abstract

This paper reports on a generalization of Lamb's problem to a linear elastic, infinite half-space with random fields (RFs) of mass density, subject to a normal line load. Both, uncorrelated and correlated (with fractal and Hurst characteristics) RFs without any weak noise restrictions, are proposed. Cellular automata (CA) is used to simulate the wave propagation. CA is a local computational method which, for rectangular discretization of spatial domain, is equivalent to applying the finite difference method to the governing equations of classical elasticity. We first evaluate the response of CA to an uncorrelated mass density field, more commonly known as white-noise, of varying coarseness as compared to CA's node density. We then evaluate the response of CA to multiscale mass density RFs of Cauchy and Dagum type; these fields are unique in that they are able to model and decouple the field's fractal dimension and Hurst parameter. We focus on stochastic imperfection sensitivity; we determine to what extent the fractal or the Hurst parameter is a significant factor in altering the solution to the planar stochastic Lamb's problem by evaluating the coefficient of variation of the response when compared with the coefficient of variation of the RF.

## 1. Introduction

Wave propagation through heterogeneous media, owing to its relevance in all of the natural sciences and most engineering fields, is an important area of fundamental and applied research. While a range of various stochastic models has been developed [1,2], the approximations do not hold when random fields (RFs) with fractal and, in addition, Hurst characteristics are involved. Indeed, such media are typically encountered in nature [3] and a quite novel mathematical tool to represent and rapidly simulate them has recently been developed in the form of Dagum and Cauchy RFs [4,5]. In fact, the fractal and Hurst characteristics can not only be captured with them, but also decoupled. In this paper, we use these RFs to model mass density in an elastic half-space subjected to a transient normal line load on its surface.

Classical Lamb's problem is defined as an elastic, homogeneous, isotropic half-space subjected to a transient, normal line load at the surface [6–9]. We generalize that problem by keeping the same geometry of the system, half-space subjected to a normal line load, but adding disorder to the mass density field. Other studies of heterogeneous Lamb's problems exist, but the material properties are a function of space or depth [10,11] or periodic medium [12]. In this study, as the RFs have fractal and Hurst characteristics, they better represent the fields encountered in nature [3,13,14].

Effectively, this is a stochastic Lamb problem where we want to determine whether the fractal dimension, the Hurst parameter, or some combination thereof, is the most significant in altering the wave motion, when compared with the constant mass density case of the classical Lamb's problem.

While a deterministic field equation can be written in the form
**f** is a source or forcing function, and **u** is a solution field or displacement vector. A stochastic field equation, by introducing randomness to the differential operator, is expressed as
*ω* is a single outcome of sample space, *Ω*, and indicates randomness, defined over a probability space *σ*-algebra, *P* denotes a Gaussian measure. In other words, *B*(*ω*) parametrized by sample events *ω* of the *Ω* space. The basic objective is determination of the ensemble average response 〈**u**〉, where 〈⋅〉 denotes stochastic expectation, which is defined as
**u**〉 satisfies
*ρ*, defined over a probability space *σ*-algebra, *P* denotes a truncated Gaussian measure, and realizations are on a domain *ρ*(**x**)>0 ∀ **x**.

In this paper, we introduce spatial randomness to the mass density field, thus leading to a stochastic Navier equation of the form
*t* is time, **x** is the undeformed coordinate, **u** is the displacement vector, *λ* and *μ* are Lamé's parameters, *ρ* is the mass density which is an RF in the material domain, *Ω* is the sample space of elementary events [15,16]. This problem can be studied analytically—usually in the frequency domain in terms of the stochastic Helmholtz equation—through perturbation, Rytov's or other methods and approximations [17,18]. Unfortunately, the correlation functions of Cauchy and Dagum RFs lack explicit Fourier transforms (i.e. spectral densities) making analytical solutions practically impossible, except for one-dimensional, elastostatic problems [19,20]. Therefore, the only way to proceed is by running Monte Carlo simulations of the stochastic Lamb's problem. That is to say, we shall draw our conclusions from the results of a relatively large number of RF realizations (or RFs generated).

This paper is organized as follows: §2 introduces the necessary concept of RFs, §3 introduces cellular automata (CA), §4 describes the geometry of the problem and properties, §5 describes the results and §6 gives the conclusion.

## 2. Random fields

In this paper, we study Lamb's problem on RFs of mass density having a trivial (or white noise) and non-trivial covariance structure. The white noise case is considered a reference case to help draw conclusions from the non-trivial covariances. Considering that a wide range of natural phenomena including geological formations (e.g. [3]) possess fractal and Hurst characteristics, two very interesting RF models have been developed over the past 20 years: Cauchy & Dagum [5,21]. In general, the fractal dimension can be described as a roughness measure that ranges from 2 to 3 in value for two-dimensional systems. Typically, the larger is the fractal dimension the rougher are the field's realizations. The Hurst parameter, or Hurst exponent, defines the long-range persistence of a system. In general, if *H*∈(0,0.5) the system is said to be anti-persistent (e.g. an increase in value is typically followed by a decrease or vice-versa). If *H*∈(0.5,1), the system is said to be persistent (e.g. an increase is followed by an increase, or a decrease is followed by a decrease). For *H*=0.5, the system is said to reflect a true random walk without any long-range persistence.

Perhaps not surprisingly, given the richness, or complexity, of patterns that can thus be reproduced, Cauchy and Dagum models lack explicit Fourier transforms, so that spectral densities are not available in algebraic closed forms. Recent work has examined rods and beams with such RF properties [19,20]. Explicit elastostatic (i.e. not even elastodynamic) analyses, involving very lengthy formulae, have been possible in such one-dimensional systems, but in two- and three-dimensional one needs to resort to numerical solutions.

RFs are characterized by their covariance function, *Z*(**x**_{1}) and *Z*(**x**_{2}) are random variables [15]. The covariance function, *C*(**x**_{1},**x**_{2}), measures the strength of linear correlation between the two random variables, *Z*(**x**_{1}) and *Z*(**x**_{2}). If the random variables are uncorrelated, then *C*(**x**_{1},**x**_{2})=0. If *Z* is Gaussian, then *Z*(**x**_{1}) is independent of *Z*(**x**_{2}) if and only if *C*(**x**_{1},**x**_{2})=0.

In the following sections, we work with RFs for which *Z* can be taken as isotropic and wide-sense-stationary (WSS). Here, the WSS assumption restricts the mean to be constant and the auto-correlation function, and consequently the covariance function, to be independent of shifts in *r*=∥**x**_{1}−**x**_{2}∥. We consider two types of RFs, correlated and uncorrelated. Uncorrelated or white noise RFs have the following covariance function:
*δ* is the Kronecker delta function. White noise RFs are known as uncorrelated RFs because *C*_{WN}(*r*)=0 ∀ *r*≠0. White noise RFs are considered in this study to serve as a control and to provide a benchmark to compare with the response Cauchy and Dagum RFs.

The Cauchy RF covariance function is

The Cauchy and Dagum RFs are unique in that they capture and decouple spatial features with a fractal dimension (*D*) and Hurst parameter (*H*). The relationships linking *D*, *H* with *n*, *α*, and *β* for both Cauchy & Dagum RFs [22] are given by
*n*, for two dimensions, is equal to 2. Unfortunately, owing to the restrictions on *α* and *β*, Dagum RFs cannot capture as many pairs of *D* and *H* as compared to Cauchy RFs. The RFs are generated in this study using the RFs package of R [23].

## 3. Cellular automata

### (a) Basics

John von Neumann first introduced CA [24] and Leamy *et al.* adapted CA to elastodynamics [25,26]. As we see below, for a homogeneous system, discretizing with a rectangular grid, CA's governing interior equations are mathematically equivalent to subjecting classical elasticity's governing equation (in displacement formulation) to the central difference finite difference method. However, the governing partial differential equations of classical elasticity are absent in CA. CA also produces straightforward treatment of common boundary conditions based on physical considerations.

The state of each rectangular cell is dependent on the state of cells that share an edge or vertex with the cell. A cell's state is defined as the cell's deformation and velocity, in both the *x*- and *y*-directions. In CA, traction boundary conditions require an additional layer of cells assigned to the defined stress. Displacement boundary conditions are assigned directly.

For convenience, we introduce the following abbreviation for finite difference operators:
*i*,*j*) denote the position of CA cells, *m* is either *x* or *y*, Δ*x* and Δ*y* is grid spacing and, *u*_{x} and *u*_{y} are the *x* and *y* components of vector **u**, respectively.

To further emphasize the similarities of the finite difference method and CA, an example of the governing equation is given below. The balance of linear momentum in the *y*-direction for the special case of homogeneous Lamé's parameters, *λ* and *μ*, reproduced from [26] is given by
*z* is the thickness, and *ρ* is the mass density. The velocity in the *y*-direction is *v*_{y} and a derivative with respect to time is denoted by an overdot. Tensile forces on the top and bottom faces are denoted by

### (b) Time stepping algorithm

The Euler method [26] is used to advance in time and it is given in table 1. The mass density at point (*i*,*j*) is *ρ*(*i*,*j*), and at timestep *n*, the internal force is *m* is either *x* or *y*. The form of the internal force, in the *y*-direction, is given in equation (3.3). The external force is

## 4. Problem parameters

Lamb's problem, a half-space subjected to a normal impulse-type load is considered here (figure 1). This problem has been studied analytically and experimentally [6–9]. Since this problem generates the Rayleigh (R), or surface wave, in addition to the pressure (P) and shear (S) waves, it is an ideal problem to consider. In this paper, we consider a two-dimensional plane-stress approximation to the problem.

The ideal homogeneous, material modelled is CR-39 as it is used by Dally in an experimental study [6]. A study comparing the experimental results with the results from CA and a non-local numerical method for homogeneous mass density fields can be found in [27]. Nishawala *et al*. [27] also includes a convergence study. CR-39 has an Young's modulus of 3.85 GPa (559 ksi), Poisson ratio of 1/3 and homogeneous mass density 1300 kg m^{−3}. This value of mass density is the mean value used in the RFs below. Pressure, shear and Rayleigh wave speeds are 1826, 1054 and 969 m s^{−1}, respectively. The input load for the numerical methods is a triangular pulse as shown in figure 2. The load is applied over a region six nodes wide. The applied impulse load has a pulse width of 20 μs and an amplitude of about 20.7×10^{3} N.

The domain for our numerical models is 400/1024 m wide (about 0.39 m), 200/1024 m high (about 0.195 m) with a thickness of 6.655 mm (consistent with [6]). The impulse is applied at the midpoint of the top edge, as shown in figure 1. Over 92 μs, the end time of the simulation, for a homogeneous system we expect the pressure wave to travel about 0.17 m, which is well within the domain. In other words, reflections are not of concern. The bottom corners are fixed to prevent rigid body motion. With grid spacing of about 0.977 mm, the domain is discretized with 400 nodes in the *x*-direction and 200 nodes in the *y*-direction. To simulate the two-dimensional geometry, there is only one node in the *z*-direction. A timestep of 0.125 μs is used. For this study, we only look at the right-travelling waves.

Alongside the results, the analytical solution, using the equations of classical continuum mechanics and elasticity, referred to as ‘theoretical’ are plotted, along with the CA response to the homogeneous system. Also plotted are vertical black lines to mark the expected locations of the pressure (P), and surface or Rayleigh (R) waves for the homogeneous system.

For each *α*, *β* pair, for either Cauchy or Dagum RF used below, 128 realizations are used to evaluate the response. For each realization, we apply the same external force and capture the surface response at 92 μs. We then calculate the mean and standard deviation (s.d.) of the response as a function of *x* along the top surface as shown in figure 1. Using this information, we then calculate the coefficient of variation of the response (CV_{R}) using the following relationship:
_{R}, loses its meaning. Therefore, we also calculate the signal-to-noise ratio of the response (SNR_{R}) which is given by
_{R} is greater than the coefficient of variation of the RF (CV_{RF}) through the Monte Carlo method with 128 realizations. In other words, can the disorder of the response be greater than the disorder of the system and if so, does it depend on the fractal or Hurst coefficients?

The relationship for CV_{RF} is given by
_{RF} and mean_{RF} are the standard deviation and mean of the RF, respectively. Both the s.d._{RF} and mean_{RF} are values chosen when the RFs are generated using the R programming language.

## 5. Results and discussion

### (a) White noise response

CA is subjected to white noise (WN) RFs. The mass density is taken as a WSS RF, as defined above, with the mean 1300 kg m^{−3} with various levels of fluctuation. This mean value, corresponding to a photoelastic medium of [6], is chosen so as to make a clear comparison with the previous results [27], for a perfectly homogeneous medium, possible. Figure 3 shows examples of heterogeneous mass-density fields used here, with a CV_{RF}=0.124. Each figure is a different realization and different coarseness. *Coarseness* is defined as the ratio between the mass-density field and the CA node density. For example, for a ‘16 by 16’ field, a 16 by 16 square of CA nodes shares the same mass density. Three different CV_{RF} are used for each coarseness. The RFs have a CV_{RF} of 0.124, 0.0877 or 0.062. For each coarseness and each CV_{RF}, 128 realizations are generated.

Figure 4 shows the results for CV_{RF}=0.124 and coarseness of 4 by 4, 16 by 16 and 64 by 64. Vertical black lines mark the expected locations of the pressure (P) and surface or Rayleigh (R) waves for the homogeneous system. Figure 4*a* plots the mean and s.d. of the response along with the homogeneous and theoretical results. Figure 4*b* shows the CV_{R} plot. Noting that since CV_{R} loses its meaning when the mean goes to zero, we plot the SNR_{R} in figure 4*c*. We see that for fine fields, the mean response of 128 realizations is very close to the homogeneous response and we also see the s.d. increase with coarseness. Hence, the CV_{R} increases and the SNR_{R} decreases as the fields become coarser. For the coarsest field, we see that CV_{R}>CV_{RF}.

The effect of CV_{RF} as well as coarseness can be seen in figure 5. Figure 5*a* depicts the minimum positive value and the maximum negative value of CV_{R} as a function of coarseness. The minimum positive value is associated with the pressure wave and the maximum negative value is associated with the Rayleigh wave. Figure 5*b* plots the maximum and minimum SNR_{R} values as a function of coarseness. As shown in figure 4, as the RFs become coarser, CV_{R} tends to increase and SNR_{R} tends to decrease. Furthermore, as CV_{RF} decreases, CV_{R} decreases as well. Also, for relatively fine RFs we see that CV_{R} for the P wave is greater than CV_{R} for the R wave. The relationship reverses for coarser RFs. Furthermore, the relationship between CV_{RF} and CV_{R} seems to be independent of CV_{RF}. The one notable exception is the CV_{R} for the 64 by 64 P wave for CV_{RF}=0.062.

### (b) Cauchy random field response

Figure 6 shows nine different realizations of Cauchy RF with CV_{RF}=0.124, for different settings of the parameters *α* and *β*. In figure 6, (*a*,*d*,*g*) is *β*=0.2 (*H*=0.9), (*b*,*e*,*h*) *β*=1.0 (*H*=0.5) and (*c*,*f*,*i*) *β*=1.8 (*H*=0.1). (*a*–*c*) is *α*=1.8 (*D*=2.1), (*d*–*f*) *α*=1.0 (*D*=2.5) and (*g*–*i*) *α*=0.2 (*D*=2.9). From the figures, we can clearly see the effect of changing fractal dimension and Hurst parameter. The figures towards the top-left have clear regions of high and low values, whereas the figures towards the bottom-right have large and small values better distributed.

Figure 7*a* plots the mean and s.d. of the response of 128 realizations for *α*=1.8 and for several values of *β*. Figure 7*b*,*c* plots the CV_{R} and SNR_{R}, respectively. The electronic supplementary material contains more plots. From the figures, we see that as *β* increases, CV_{R} decreases. Comparing figure 6*a*–*c* with white noise RFs in figure 3, figure 6*a* shows that there are significant areas of high and low density clustered together, similar to a coarse white noise RF. As *β* is increased, as given in figure 6*b*,*c*, the areas of high and low values become smaller, similar to a fine white noise RF. The effect of increasing *β*, or decreasing *H*, is similar to that of creating a finer RF.

In the electronic supplementary material, we provide plots for a constant *β* and changing *α*. Unfortunately, unlike the previous case, we do not see a significant change in the results as we change *α*, or *D*. The electronic supplementary material also contains a plot where *α* and *β* are held constant and we change CV_{RF}. We see that as CV_{RF} increases, CV_{R} increases, as expected. We also note that if |CV_{R}|<CV_{RF} is satisfied for one CV_{RF}, then it is also satisfied for other CV_{RF}.

Other combinations of *α*,*β* have been evaluated as well to help determine the boundary between CV_{R} being ‘less than’ or ‘greater than’ CV_{RF}. Figure 8 plots the results for over 100 values of *α* and *β*. Figure 8*a* plots the result for the pressure wave. An ‘o’ signals, for that point, the maximum SNR_{R} exceeded SNR_{RF} and a ‘x’ denotes that the maximum SNR_{R} did not exceed SNR_{RF}. A dotted line is provided to show the approximate border between the two results. The same notation is used for the Rayleigh wave figure (figure 8*b*). To see the entire *α*, *β* space, see the electronic supplementary material. The contour plot in figure 8*a* plots the maximum SNR_{R} and the contour plot in figure 8*b* plots the absolute value of the minimum SNR_{R}.

We can see that for *α* greater than, about, 0.6, there is little change in the plot as *α* increases. For *α* less than 0.6, we see that for relatively small changes in *α* and *β*, the values of SNR_{R} can change dramatically. We also note that there are more ‘x’s in figure 8*b* than in figure 8*a*. This primarily may be due to two reasons: (i) interference of the P and R waves due to back-reflections of the P wave and (ii) the elliptic motion of the particles in the R wave.

### (c) Dagum random field response

For Dagum RFs, we use the same mean and CV_{RF} as in the Cauchy case. Example RF realizations, CV_{R} and SNR_{R} plots can be found in the electronic supplementary material. The trends found in Cauchy RFs from changing the fractal and Hurst parameters also occur for Dagum RFs.

Figure 9 plots the results for over 30 *α*,*β* pairs to help determine the boundary between CV_{R} being ‘less than’ or ‘greater than’ CV_{RF}. Figure 9*a* plots the result for the pressure wave. An ‘o’ denotes, for that point, the maximum SNR_{R} exceeded SNR_{RF} and a ‘x’ shows that the maximum SNR_{R} did not exceed SNR_{RF}. A dotted line is provided to show the approximate boarder between the two results. The same notation is used for figure 9*b* for Rayleigh waves. Also included are contour plots of the maximum SNR_{R} and figure 8*b* plots the absolute value of the minimum SNR_{R}.

We can see that for changes in *α*, there is some change in the SNR_{R}. It is possible to have SNR_{R} to be less than SNR_{RF} for larger values of *α* and *β*. However, smaller values of *β* are required for a smaller *α*. Similar to Cauchy RFs, there are more ‘x’s in figure 9*b* than in figure 9*a*. Again, as stated in the case of the Cauchy RF above, this could be due to P wave reflections and/or the elliptic motion of the R wave. See the electronic supplementary material for more information.

## 6. Conclusion

This paper generalizes Lamb's problem on an elastic half-space to random mass density fields with fractal and Hurst characteristics. The motivation for this study is the fact that a multitude of media in nature fall in that category and it is important to know the sensitivity of response of this classical elastodynamic problem to spatial fluctuations relative to the homogeneous medium assumption. Solutions of this stochastic multiscale problem are carried out using the CA method. For a reference, we also evaluate the response of half-space with uncorrelated random mass density fields while varying their coarseness. For finer mass density fields, we see little to no change in the response. The mean response, of 128 realizations, is very close to the homogeneous result. As the RF becomes coarser, we begin to see changes in the response which is reflected by the increasing CV_{R}. Also note that the CV_{R} of the R wave tends to be smaller than the P wave for finer RFs and for coarser fields, the CV_{R} for the R wave is larger than the P wave.

An extensive study of the sensitivity of Lamb's problem response to material random fluctuations is carried out for both, Cauchy and Dagum RFs. The sensitivity is measured in terms of the coefficient of variation of response (CV_{R}) versus the coefficient of variation of RF (CV_{RF}), the cause of scatter and departure from the homogeneous medium. In general, the relationship between CV_{R} and CV_{RF} is independent of CV_{RF} except for RFs that are very coarse.

Overall, for smaller values of *β* the CV_{R} tends to be greater than the CV_{RF}. Since a small *β* corresponds to a high *H*, i.e. to pronounced rare events, this is expected. A higher value of *H* is similar in structure to coarse WN RFs. We also see that larger values of *α* tends to support a larger range of *β*'s, where the CV_{R} is greater than CV_{RF}. However, for *α* greater than about 0.8, the relationship looks to be independent of *α*.

Next, for larger values of *β*, or smaller *H*, we see that the mean response of 128 realizations is very close to the homogeneous result. Which implies that a wave response close to the homogeneous result does not necessarily suggest a homogeneous medium. Also, CV_{R} tends to be smaller as *β* increases. As *β* grows, similar to finer WN fields, we tend to see CV_{R} decrease faster for the R wave than for the P wave. This could be due to P wave reflections and/or the elliptic motion of the R wave.

Generally, for finer RFs, or for *β* values greater than 0.55, CV_{R} would be less than CV_{RF}. For coarser RFs, *β* values less than 0.55, and *α* greater than 0.8, CV_{R} would be larger than CV_{RF}. For *α* less than 0.8, even smaller values of *β* are needed so that CV_{R} would be larger than CV_{RF}. Effectively, in the stochastic Lamb's problem, we find that the fractal dimension has a weaker effect than the Hurst parameter on altering the wave motion, when compared with the constant mass density case of classical Lamb's problem.

## Data accessibility

Additional plots in the electronic supplementary material. Scripts used to create data can be found at https://dx.doi.org/10.6084/m9.figshare.3979056.

## Authors' contributions

V.N. and M.O.-S. conceived the problem, interpreted results and wrote this manuscript. V.N. performed all CA simulations and calculations in consultation with M.O.-S. M.L. implemented the CA numerical model. M.O.-S. and E.P. consulted on Random Field Models.

## Competing interests

We have no competing interests.

## Funding

This work is partially supported by the NSF under grants CMMI-1030940 and IIP-1362146 (I/UCRC on Novel High Voltage/Temperature Materials and Structures).

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3579824.

- Received August 15, 2016.
- Accepted November 4, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.