This post has the following aims:

• Provide documentation and source code for a spherically symmetric wave propagation in a linear-elastic medium.
• Tell a story illustrating how this simple verification problem helped to validate a complicated rate-dependent and history-dependent geomechanics model.
• Warn against believing previously reported material parameters, since they might have been the result of constitutive parameter tweaking to compensate for unrelated errors in the host code.

Many years ago, a colleague of mine had developed a complicated history-dependent and rate-dependent constitutive model for limestone. The model had been extensively verified against analytical solutions for triaxial compression and extension at a variety of rates and at a variety of confinement levels. The first serious validation test for this model was to reproduce measured data for the rock’s response to a small explosive embedded deep within the rock. The data consisted of two measured velocity histories, one a small radial distance from the explosive and the other a large radial distance from the explosive. To avoid the modeling uncertainty of predicting blast pressures, the problem was set up so that the velocity measured at the small-radius accelerometer would serve as the boundary condition at the inner radius of a thick spherical shell. The model would be validated if it could then predict the material response measured by the accelerometer located at the larger (downstream) radius.

When my colleague set up and ran this problem, he was heartbroken that the predicted velocities did not come close to matching the measured velocities. He was mystified why his model would be that far off, given that it had been very carefully calibrated to match other data in similar loading modes. Whereas poor prediction results are usually caused by bad material models, it turned out that it was the host code’s treatment of artificial viscosity that was the culprit in this case. Once this was fixed, my colleague’s simulation matched the data very well (without any parameter tweaking).

So, how was the problem with artificial viscosity exposed? By setting aside the validation problem and solving an idealized VERIFICATION problem instead. Specifically, the material model was replaced with simple linear-elastic Hooke’s law, so that the code could be checked against the analytical solution to this problem documented in David Aldridge’s 2002 government report, aldridge_sphere_source. The semi-analytical solution in that report was evaluated using software provided by the author (AldridgeSphericalSourceVerificationSolver). It was found that the solid-mechanics FEM code was unable to match Aldridge’s semi-analytical result. Ultimately, it was found that the default artificial viscosity settings in the FEM code were corrupting the prediction. When the artificial viscosity settings were revised to match this idealized problem, the same settings were used in the validation problem to obtain much improved results.

Imagine how this exercise would have played out if it were not for the existence of Aldridge’s verification problem. Unscrupulous analysts might be tempted to tweak the viscosity parameters in the material model in order to improve comparisons with data. If this had occurred, it is conceivable that these incorrect viscosity values could become regarded as “correct” for that material without adequate justification. People coming along later would have been unjustly forced to explain why their values were different from these early results, but how could they know that the problem with the first endeavor was a numerical bug? As Prof. Dave Benson remarked in 2012, “The first girl to the ball isn’t necessarily the prettiest.”  In other words, the first reported result might not be the right result (so we shouldn’t let it bias current results).

Aldridge’s main contribution was that he provided source code (downloadable from this posting) which implemented the solution in a format that any engineer could easily figure out (just by reading the “readme” files). Engineers could quickly apply Aldridge’s code to their own problems as long as they had access to a legacy FORTRAN compiler (which comes standard on most Linux systems).  Non-experts (even undergraduates) used Aldridge’s software to obtain the response functions from a broad range of boundary conditions, and they did it in typically less than one day. This is a major advance that sets a standard: we should not satisfy ourselves with merely obtaining a solution. We need to deliver our work to the scientific community in a package that doesn’t require weeks to learn.  It would be nice if someone were to transform his Fortran code into other more commonly known languages (like matlab, python, and/or C). Any takers?

Additional background information (and another report, Denny1991_ExplosionSeismicSource) contributed by Jim Kamm:

The problem [here referred to as]  “Aldridge’s verification problem” has been used at LANL (and presumably elsewhere) well before the the publication of Aldridge’s report…This problem is known within the LANL/LLNL communities as the “Blake problem”, in reference to Blake’s 1952 paper.  As cited by Aldridge and other authors, the genesis of this problem goes back further, notably to Jeffreys and to Sharpe, and, so, might be more appropriately named differently; for whatever reason, however, the “Blake problem” name has stuck at LANL/LLNL and, indeed, has for decades.  The general solution is also to be found in the well-known text of Achenbach, “Wave Propagation in Elastic Solids”, 1976.  Additionally, this problem is known within that branch of the seismology community dealing with seismic source functions;  see \S II of the attached paper by Denny & Johnson (1991). Moreover, there have been various lab reports through the years detailing mesh convergence studies of linear elastic models in hydrocodes on this problem.