ABSTRACT
Title of thesis: INVERSE SPECTRAL METHODS IN ACOUSTIC NORMAL MODE
VELOCIMETRY OF HIGH REYNOLDS NUMBER SPHERICAL
COUETTE FLOWS.
Anthony Mautino, Master of Science, 2016
Thesis directed by: Professor Vedran Lekic
Department of Geology
Experimental geophysical fluid dynamics often examines regimes of fluid flow infeasible for
computer simulations. Velocimetry of zonal flows present in these regimes brings many challenges
when the fluid is opaque and vigorously rotating; spherical Couette flows with molten metals
are one such example. The fine structure of the acoustic spectrum can be related to the fluid’s
velocity field, and inverse spectral methods can be used to predict and, with sufficient acoustic
data, mathematically reconstruct the velocity field. The methods are to some extent inherited from
helioseismology. This work develops a Finite Element Method suitable to matching the geometries
of experimental setups, as well as modelling the acoustics based on that geometry and zonal flows
therein. As an application, this work uses the 60-cm setup Dynamo 3.5 at the University of
Maryland Nonlinear Dynamics Laboratory. Additionally, results obtained using a small acoustic
data set from recent experiments in air are provided.
INVERSE SPECTRAL METHODS IN ACOUSTIC NORMAL MODE
VELOCIMETRY OF HIGH REYNOLDS NUMBER SPHERICAL COUETTE
FLOWS
by
Anthony Robert Mautino
Thesis submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Master of Science
2016
Advisory Committee:
Professor Vedran Lekic, Chair/Advisor
Professor Daniel Lathrop, Co-Advisor
Professor William McDonough
c© Copyright by
Anthony Robert Mautino
2016
Acknowledgments
First and foremost, with regard to this thesis, I would thank my advisor Ved Lekic. Ved
was always available, usually on a daily basis, to work with me. He was of considerable guidance
and help in both the research that went into this work as well as the drafting of this thesis. I thank
Ved greatly for both my research assistantships as well as arranging my teaching assistantships.
I would also like to thank the rest of my committee. Dan Lathrop has served as a co-advisor
throughout this thesis and offered much input. Bill McDonough has always made himself available
for discussing science both related and unrelated to this endeavor, as well as encouraging me to
have a more positive outlook in science.
I thank everyone else on the faculty who took time to discuss science and teach me things
while I was here, especially Phil Candela. I thank all the fellow students who gave me engaging
discussions, especially my fellow lab member Chao Gao.
I thank Todd Karwoski for all the IT support, and Michelle Montero for all her assistance
with the formalities of being a graduate student.
And of course, I thank my Mother for all her support.
ii
Table of Contents
List of Tables v
List of Figures vii
List of Abbreviations x
1 Background 1
1.1 Laboratory Analogs of Flows in the Earth’s Core . . . . . . . . . . . . . . . . . . . 1
1.1.1 Laboratory Simulations vs. Computer Simulations . . . . . . . . . . . . . . 1
1.1.2 Spherical Couette Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Parameters for Spherical Couette Flows . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1 Reynolds Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Ekman Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2.3 Rossby Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.4 ω Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.3 Velocimetry of Fluid Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.1 Challenges of Traditional Techniques of Velocimetry . . . . . . . . . . . . . 11
1.3.2 Acoustic Normal Mode Inverse Spectral Methods . . . . . . . . . . . . . . . 12
1.4 60-cm Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.1 Laboratory Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.2 Modeling the Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2 Previous Work in the Field 24
2.1 Analysis and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2 Contrast with this Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3 Mathematical Formulation 43
3.1 Acoustic Pressure Perturbation in Stationary Homogeneous Media . . . . . . . . . 43
3.2 Exploiting Axisymmetric Geometry for Computational Efficiency . . . . . . . . . . 44
3.3 Variational Form of the Acoustic Helmholtz PDE in Pressure . . . . . . . . . . . . 45
3.4 Discussion of Variational Scalar Formulation . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Modal Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4 Numerical Formulation 51
4.1 Finite Element Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2 Solving for the Stationary Modes and Stationary Frequencies . . . . . . . . . . . . 53
4.3 Discussion of FEM Predictions for Stationary Modes . . . . . . . . . . . . . . . . . 57
4.4 Solution for a Spherical Annulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4.1 Modes of the Spherical Annulus . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4.2 Computing the Solution to the Radial Wave Equation . . . . . . . . . . . . 58
4.5 Comparing Theory and FEM for Spherical Annulus Solution . . . . . . . . . . . . 60
iii
5 Linear Theory of Acoustic Splitting for Rotating Flows with Shear 63
5.1 Wave Equation inside a Vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2 Uniform Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.3 Soundspeed Uniformity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6 Relating Fluid Rotation to Acoustics 71
6.1 Flow Characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.2 Parametrization Style for the Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.3 The Forward Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.4 Boundary Conditions of the Fluid Flow . . . . . . . . . . . . . . . . . . . . . . . . 74
7 Techniques of Inversion 75
7.1 Constrained Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2 Minimizing Spread Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.2.1 Basic Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.2.2 Dirichlet Spread Function Minimization . . . . . . . . . . . . . . . . . . . . 81
7.2.3 Backus and Gilbert Spread Function Minimization . . . . . . . . . . . . . . 81
7.3 Adaptive Parametrizations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.3.1 Adaptive Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.3.2 Non-coefficient Model Parameters . . . . . . . . . . . . . . . . . . . . . . . 85
8 Applications to Ambient Noise Data 87
8.1 Data Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
8.1.1 Array Configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
8.2 Parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.2.1 Approach to Parametrizing Fluid Rotation . . . . . . . . . . . . . . . . . . 101
8.2.2 Cylindrically Radial Cubic B-Splines and Neumann Modes . . . . . . . . . 103
8.2.3 Choice of Reference Frame . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.3 Processed Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.4 Inversions for Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.4.1 Constrained Optimization Inversions . . . . . . . . . . . . . . . . . . . . . . 109
8.4.2 Consistency with Alternate Basis Functions . . . . . . . . . . . . . . . . . . 125
8.5 Alternate Methods of Inversions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.5.1 Minimization of Spread Functions . . . . . . . . . . . . . . . . . . . . . . . 129
8.5.2 Interpolating on a Fixed Grid . . . . . . . . . . . . . . . . . . . . . . . . . . 132
8.6 Mean Flow Approximation Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
8.7 Alternate Basis Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
9 Concluding Discussion 136
9.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
9.1.1 Application of the Finite Element Methods . . . . . . . . . . . . . . . . . . 136
9.1.2 Bulk Fluid Rotation Recovery . . . . . . . . . . . . . . . . . . . . . . . . . . 136
9.2 Looking Forward . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
9.2.1 Liquid Sodium Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
9.2.2 Poloidal Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
9.3 Recommendations for Future Improvements . . . . . . . . . . . . . . . . . . . . . . 138
9.3.1 Inversion Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
9.3.2 Poloidal Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
9.3.3 Modeling Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
9.3.4 Meshing Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
9.3.5 Test Functions and Shape Functions . . . . . . . . . . . . . . . . . . . . . . 140
9.3.6 Inclusion of Magnetic Data in Constraints . . . . . . . . . . . . . . . . . . . 141
9.3.7 Adaptive Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
Bibliography 142
iv
List of Tables
4.1 FEM calculated frequencies at c=344 m/s for the 60-cm setup. As can be seen, the
multiplets are split due to the shaft’s presence. Values have been rounded to the
nearest natural number. The data the methods are applied to is as high as 2.1kHz;
the reason for the range presented. . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2 Normal Mode Frequencies computed by FEM in column 3 agrees quite well with
the theoretical semi-numerical result in column 2, as seen in the relative errors in
column 4. The differences between frequencies of the elements in the multiplets are
quite small as well; column 5. The NA values are reflecting that only ms = 0 is
possible when l = 0. Note: the estimated relative errors are based on the unrounded
values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.1 A comparison of my FEM predictions with observed frequencies in the stationary
medium by targeted chirp, and errors are a small fraction of 1%. The error is likely
even slightly lower, since the soundspeed was calibrated using the lowest frequency,
rather than the highest, or some least squares fit to all. The measured acoustic data
in the column for observed frequency were performed by Matthew Adams. . . . . . 67
5.2 Comparison of the same mode set from the table 5.1 , at various rotation rates
where the inner and outer boundaries are spinning together, spun up to achieve
uniform rotation of the fluid. The ∆ values are the absolute shift of the modes from
the center frequency in Hz. The denominators on the observation ∆ values are the
rate of medium’s rotation in Hz. The leftmost column represents the expected shift
from center frequency due to 1 Hz uniform fluid rotation.The percent errors do not
take any account of the observational uncertainties or modeling uncertainties. The
missing data in the first row is because measurement could not be reliably made due
to that portion of the spectrum at the high rotation rate of the experiment being
very noisy. Under the presumption that there should be no splitting at 0 Hz then,
the linearity implies that the ratio of the splitting to rotation rate should remain
invariant during experiments and equal to the finite element prediction. Note: This
table does not give signs to the shifts, since signs were not measured for the data,
only the magnitude of the splitting was measured. Acoustic data measurements
made by Matthew Adams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.1 Mode wavenumbers and stationary frequency prediction at c=344 m/s. These are
the modes used to perform the inversions. Additionally, the frequency shift from
center frequency due to uniform rotation with the outer shell at 6Hz has been
provided, rounded to 1 decimal place, since it does not appear in the presented
splitting data, it has been subtracted off. The negative sign indicated that −|ms|
has a higher frequency than +|ms|. Values rounded to nearest tenth of a Hz. The
splitting is thus approximately twice these listed values. The first column provides
the index used in plots of data fits in this chapter, for reference. They are always
presented in this order herein. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
v
8.2 Frequency shift from center frequency in Hz. Top row is Fi in Hz, and each row
is a mode pair given by the first column. Measurements are in Hz. Outer Sphere
is rotating at 6 Hz and the inner sphere rotation in Hz is given at the top of each
column. Each row is for one particular mode. The data has had the correction
due to being in a rotating frame removed from it. For instance, the first row shows
(0,1,1) has relatively small shifts, but its shift due to the Coriolis effect of being in
a rotating frame is almost 4 Hz. The value stated is the residual after that amount
was subtracted off. Refer to table 8.1 for the approximate Coriolis correction of the
rotating frame that has been removed. . . . . . . . . . . . . . . . . . . . . . . . . . 106
8.3 Errors (Hz) relative to the mean that are used for performing inversions. Table
layout presented in the same form as table 8.2. Values are rounded off to nearest
tenth in the table. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
8.4 Fraction of measurement that the error radius represents, as a percent. Often the
uncertainty is small relative to the absolute shift that includes the shift induced by
the rotating frame, but much larger relative to the residual shift. Table layout same
as in table 8.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
vi
List of Figures
1.1 Schematic of the 60-cm setup and configuration in the lab. Diameter of the outer
sphere is approximately 60 cm and of the inner sphere is a approximately 20 cm.
Diagram courtesy of Matthew Adams. . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.2 Interior of the 60cm Setup taken apart showing electronics mounted inside, before
assembly. The setup tries to minimize the invasiveness of the electronics. The axis of
rotation, embodied by the physical axle, is vertical during operation. Photographs
courtesy of Matthew Adams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.3 Meridional cross section of the geometry, showing the proportions used to model
the cross section. The z-axis is the axis of rotation where the axle and shaft are. . 23
2.1 Reproduction of a figure from Triana et al. [2014], showing anemometer profile and
level of consistency with inversion profiles, as well as posterior uncertainty estimates.
s is cylindrical radius, ro is the outer sphere radius. . . . . . . . . . . . . . . . . . . 30
2.2 Reproduction of a figure from Triana et al. [2014], averaging kernels for several
locations used in the inversion. The white cross hairs are showing the location in the
domain, which has been broken into thousands of small tiles as basis functions, and
the colors are showing the averaging kernel for it. Note the high non-local diffusion,
especially for the bottom left sub-figure. The reconstructed rotation rate at the
point marked, will be a weighted average of the rotation rates at other locations in
the domain, and the weights are given by the coloring of the domain. . . . . . . . 33
3.1 Meridional cross section of pressure modes (1,0,0) and (2,0,0). The color represents
the relative amplitude. The nodal curves in white. . . . . . . . . . . . . . . . . . . 47
3.2 Meridional cross section of pressure modes (0,1,0) and (0,2,0). The color represents
the relative amplitude. The nodal curves in white. . . . . . . . . . . . . . . . . . . 49
3.3 Meridional cross section of pressure mode (1,3,2). The color represents the relative
amplitude. The nodal curves in white. Since m=2 is known from how the PDE is
solved, the nodal curves can be counted. See that n=1, and l − |ms| = 1, and so
l = 3. No inconsistencies or ambiguities result from adopting this naming convection
in the altered geometry within the frequency range and small wavenumbers used in
this study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
8.1 Comparison of acoustic spectrum with different frequency resolutions. The black
curve has been created with a frequency resolution of 0.3 Hz whereas the red curve
has a frequency resolution of 0.05 Hz. Note: magnitudes are linear, not in decibels.
The clusters of peaks could be due to poloidal flow, noise envelopes, mechanical
noise, changing mean flow state during time integration, or some combination thereof. 91
vii
8.2 Power spectra plot of acoustic frequency(horizontal axis) vs. inner core rotation
rate ( vertical axis ) near the modes (0,1,1) and (0,1,-1). Color is relative amplitude
in log scale. The spectra is computed by taking the difference at two microphones
at equal latitude spaced apart by 180 degrees in terms of azimuth. This shows drift
as the temperature increases inducing center frequency shift, and shows that the
dominant effect on splitting is the outer spinning at 6 Hz rather then the greater
splitting that occurs due to additional rotation of the inner core. There could also
be some effect of poloidal flow. In this regime centrifugal force is not expected to
have any significant effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
8.3 Similar figure and axes as in figure 8.2. Power spectra plot of power spectral density
vs. inner core rotation rate near the mode(0,1,0). This shows drift as the temper-
ature increases inducing center frequency shift. Color is relative amplitude in log
scale. The spectra is computed by taking the sum of two microphones at equal
latitude spaced apart by 180 degrees in terms of azimuth. . . . . . . . . . . . . . . 96
8.4 Power spectra plot of frequency vs. inner core rotation rate near the mode (2,3,1).
Intensity of the plot is power spectral density, relative amplitude in log-scale. Note
splitting changes with rotation rate, and also the drift of the center frequency. . . . 97
8.5 Power spectra plot (difference between microphones) of power spectral density vs.
inner core rotation rate near the modes (0,2,1) and (0,2,-1). This shows drift as the
temperature increases inducing center frequency shift, and shows that the dominant
effect on splitting is the outer spinning at 6 Hz rather then the greater splitting
that occurs due to additional rotation of the inner core. The lower plot is a power
spectral density estimate placed on a linear scale, for inner core rotation rate of -40
Hz. It is a relative amplitude, not calibrated to any particular loudness level. . . . 99
8.6 RMS misfit of the ratio of the discrepancy between predicted frequencies from inver-
sions and the measured values of frequencies to the error weights. The rotation rate
of the outer sphere was 6 Hz, and the inner sphere’s rotation rates were -4,-6,...,-40
Hz and shown in the horizontal axis. For the purposes of assessing overfitting vs.
underfitting the data, the vertical axis values should be scaled by a factor of between
2 and 4 for converting to uncertainty on the mean. Any choice of a particular value
for scaling would be somewhat arbitrary. . . . . . . . . . . . . . . . . . . . . . . . . 109
8.7 Histogram with Gaussian fit for an ensemble of 5000 inversions for Ro=-4.67, for
the angular momentum quantity described in figure 8.8. The distribution was suffi-
ciently Gaussian and similar in standard deviation and mean for smaller ensembles
that it allowed for sufficient accuracy in the weighting of the data when performing
fits. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
8.8 The azimuthal angular momentum as a percent of the uniform rotation case at
6Hz, depicted by triangles. The variation with respect to sampling the data within
the uncertainty interval around presented as error bars, and the curve depicting a
purely qudratic relationship with Rossby number after offsetting by 100% at uniform
rotation, i.e. Ro=0. The coefficient in the quadratic term is ≈ −0.62 . . . . . . . . 112
8.9 The fit of %V6 is a second order polynomial in Ro with an intercept of 100% at Ro=0.114
8.10 The fit of %Ω6 is linear but does not go through the 100% at Ro=0. . . . . . . . . 115
8.11 Inferred ωeffect, described in chapter 1, as a function of Ro. . . . . . . . . . . . . . 116
8.12 Estimating the location of high shear. The horizontal axis represents the location
of the vertical line separating the two piecewise constant basis functions, as inferred
by minimizing the misfit to the acoustic data; described in text. Though a line is
also possible to fit the data nearly as well, the square root function looked like a
better fit and also goes through zero, which is consistent with solid body rotation of
the fluid. The vertical axis has been non-dimensionalized by dividing by the radius
of the outer sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.13 Data misfit for the shear zone parametrization by a sigmoid, using similar layout as
figure 8.6. This corresponds to the data plotted in figure 8.13. . . . . . . . . . . . . 118
viii
8.14 Reconstructed inversions at each rotation rate. The rotation rate of the inner sphere
has been placed on top of each profile, but all other labelling has been omitted to
reduce clutter. The outer sphere is spinning at 6 Hz. . . . . . . . . . . . . . . . . . 120
8.15 Reconstructed flow profile for Ro=-4.67, i.e. Fi = −22 Hz and Fo = 6 Hz. There
appears to be a jet forming on the equator of the inner sphere. . . . . . . . . . . . 121
8.16 Data prediction from reconstructed flow profile for Ro = -4.67, i.e. Fi = −22 Hz
and Fo = 6 Hz. The data is plotted in red with the error weights used in the inversion.122
8.17 Relative variability profile for Ro=-4.67, i.e. Fi = −22 Hz and Fo = 6 Hz. The units
of Hz is consistent with the mean flow error radius. This is only variation across an
ensemble though, not necessarily reflective of the uncertainties related to how well
the basis can match the flow in very localized regions. . . . . . . . . . . . . . . . . 123
8.18 Reconstructed flow profile for Ro=-4.33, i.e. Fi = −20 Hz and Fo = 6 Hz. Any jet
is less pronounced in this inversion than in the one for Fi = −22 Hz. . . . . . . . . 124
8.19 Reconstructed Profiles using only 4 cubic splines, as discussed in text. Recovered
profiles look similar to how the profiles in figure 8.14 would look if they were averaged
vertically and homogenized in the z variable based on that average. Figure similar
labelling convention as figure 8.14. . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
8.20 Data misfit analysis for basis consisting of 4 cubic splines, discussed in text. Figure
similar to figure 8.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
8.21 Result of minimizing Dirichlet spread function for 2D cubic spline basis. It has many
realistic features, but also has some overshoot in the sense that it super-rotates the
outer shell in a way that does not seem plausible. . . . . . . . . . . . . . . . . . . . 130
ix
List of Abbreviations
R spherical radius variable
Ro spherical radius of outer shell
Ri spherical radius of inner core
r cylindrical radius from axis of rotation
z vertical signed distance from the equatorial plane
φ azimuthal angle
θ polar angle, usually colatitude
ρo density
c soundspeed
P fluid pressure
p acoustic pressure perturbation
ψ pressure perturbation at zero azimuth
v fluid velocity field
u acoustic material displacement perturbation
ω acoustic angular frequency
f acoustic frequency in Hz
δω perturbation in angular frequency
Ω fluid angular velocity
Ω fluid rotation rate in Hz
Ωi inner sphere angular velocity
Ωo outer sphere angular velocity
Fi inner sphere rotation rate in Hz
Fo outer sphere rotation rate in Hz
k radial wavenumber
n overtone number
l polar degree
ms azimuthal order
Re Reynolds’ Number
Rm Magnetic Reynolds’ Number
Ro Rossby Number
Ek Ekman Number
ω omega effect
C forward operator
d data vector
m model vector
Rm model resolution matrix
F Fourier transform operator
BFM Barrier Function Method
FEM Finite Element Method
FDM Finite Difference Method
FWHM Full Width at Half Max
MHD Magnetohydrodynamics
NSE Navier-Stokes Equations
OLA Optimally Localized Averages
PDE Partial Differential Equation
RMS Root Mean Square
SCF Spherical Coutte Flow
x
Chapter 1: Background
1.1 Laboratory Analogs of Flows in the Earth’s Core
1.1.1 Laboratory Simulations vs. Computer Simulations
Experimental Fluid Dynamics provides simulations inaccessible or unmanageable on the
computers currently available; synthesizing flows with high Reynolds numbers and high magnetic
Reynolds numbers appropriate to test theories of the Earths core is not feasible [Lathrop and Forest,
2011]. For instance, numerical simulations of flows in the Earth’s core use unrealistically high
viscosity for the actual regimes they are modeling in order to maintain numerical stability [Dormy
et al., 2000]. Early numerical models such as the canonical Glatzmaier and Roberts [1995] model
were able to recover the Earth’s dominant magnetic axial dipole, similar power spectra we observe
on Earth, and other observed phenomena [Glatzmaier et al., 2004] while simulating flows lacking the
turbulence thought to be present in the Earth’s outer core. As the numerical models improved with
computing power and more of parameter space was explored computationally, parameter values
closer to that of Earth were used. The magnetic fields they produced in the simulations were
not more ”Earth-like”, but were actually less “Earth-like”[Christensen et al., 2010]; lowering the
viscosity relative to the Earth’s modeled rotation rate and size more closely to the Earth’s theorized
ratios, made the magnetic fields overly dominated by the axial dipole not dipole dominated at all.
The ratio of kinematic viscosity to Coriolis forces is termed the Ekman number Ek, and is
discussed herein. It is thought to be possibly as low as Ek = 10
−15 for the Earth’s outer core,
yet computer simulations typically employ values at least 8 orders of magnitude larger, and the
“Earth-like” behavior occurring at 10 orders of magnitude higher near Ek = 10
−4 [Christensen
1
et al., 2010]. Though the kinematic viscosity in the outer core is one parameter that is not known
very accurately, or at least not widely agreed upon [Secco, 1995], most estimates yield an Ekman
number much smaller by many orders of magnitude than is used in numerical simulations that
produce Earth-like magnetic fields.
Numerical modeling on the computer is considered to be many orders of magnitude away
from the Earth in terms of several other key parameters: Reynolds number Re and Magnetic
Reynolds number Rm. The Reynolds number is a measure of the ratio of inertial forces to viscous
forces, and the magnetic Reynolds is a measure of the ratio of induction to magnetic diffusivity. In
the numerical simulations, added viscosity and magnetic diffusivity in the form of a lower Reynolds
and magnetic Reynolds number can suppress characteristic features that would be present in real
flows, and also obscure whether features observed in numerical simulations would persist in the
presence of turbulence and lower viscosity; properties associated with authentic representatives of
the regimes being studied. Similar problems exist for the magnetic simplifications made as well,
but the issues with numerical simulations are not purely due to computing constraints.
Christensen et al. [2010] point out that even though computing power available places
bounds and constraints on the parameters and parameter combinations employed in Magneto-
hydrodynamic (MHD) simulations of the core, these parameters can be varied within the com-
putationally feasible region of parameter space. Yet, many published numerical simulations and
models do not justify the exact values used in their simulations. It is not clear from the mathe-
matics that the phenomena produced by such models are expected to persist within a reasonable
neighborhood in parameter space of the parameter values chosen. Certainly, parameters that can
be explored are limited by computational limitations, but their precise selection might be made to
reproduce Earth-like behavior.
Another approach to simulating hypothetical models of core dynamics is using laboratory
analogs. Laboratory models that perform actual experiments with real fluids. Experimental fluid
dynamic models of planetary dynamos, though imperfect analogs as well, approximate core dy-
namics more closely than computers can by exploring a much vaster region of parameter space.
2
Laboratory experiments such as those carried at University of Maryland’s liquid sodium experi-
ments, far surpass computer simulations in matching many of the characteristic parameters needed
to test theories of the Earth’s core [Lathrop and Forest, 2011].
It remains an unsolved problem how the Earth’s magnetic field is sustained. Laboratory
MHD experiments seek to replicate key features of a planetary dynamo, a self generating and
sustaining magnetic field induced by the flowing electrically conductive fluid as in the Earth’s outer
core; an output of magnetic energy sustained by the kinetic energy of fluid flow. For instance,
Parker [1955] described a model whereby a feedback process in the Earth’s core could produce
a persistent magnetic field in the presence of differential rotation. It is thought that chemical
convection in the outer core causes buoyant material to rise radially from the Earth’s center, and
is deflected azimuthally by the Coriolis effect induced by the Earth’s rotation, then in turn swirled
in helical patterns as the process evolves. It is also thought that to a lesser extent the convection
is driven by torque-induced precession as well as tidal forces and thermal convection. At high Rm
the magnetic field lines are frozen in to the fluid, and are carried with the rotating flow.
For magnetohydrodynamic experiments related to core studies, it is especially important
then to be able to infer the azimuthal flow so as to quantify the toroidal magnetic field induced
from the twisting of poloidal magnetic field lines crossing through the core, the so called ω effect.
The entrainment of frozen-in poloidal magnetic field lines by a flow can induce deformation of
magnetic field lines by deflecting them in the toroidal direction as the field lines are sheared and
carried by the vortex’s differential rotation. The α effect—the other component of the feedback
loop—being the reconversion of toroidal field lines back into poloidal field lines; from the helical
winding as it is carried with the convecting flows. Inferring the toroidal component of the internal
magnetic field from sensors on the boundary is difficult magnetically, but can be inferred from the
azimuthal fluid velocity field if it can be measured.
Large-scale features of the flow, such as average rotation rate in the mean zonal flow, ex-
tent of shear zones, super-rotation, and large-scale profiles of specific angular momentum in the
azimuthal component of the mean flow, correlate with magnetic phenomena. The ability to char-
3
acterize the fluid dynamics within the system would offer a means to analyze what aspects of the
fluid dynamics enhance magnetic field gain. Such an ability thereby allows a more directed search
of parameter space when studying planetary dynamos. MHD rotational flows are also relevant to
the study of the stars and our sun. Producing these types of flows in the lab by convection with
the same means the Earth’s core and stars do would be difficult, and more importantly, would
evolve quite slowly.
Spherical Couette Flows (SCF) offer another means by which similar flows can be established
and evolved in the lab, on realistic time scales. SCF are driven by two rotating concentric spherical
boundaries encapsulating the fluid medium, and offer a means through which relevant flows can
be realized in conductive fluids such as molten metals. Sodium is a commonly used metal since
it is highly conductive electrically, and sodium melts at approximately 98 degrees Celsius; it is
also quite cheap. Plasma-based experiments are also possible and offer a larger value of magnetic
Reynolds number; liquid sodium offers a better analogue to planetary dynamos whereas plasmas
may better approximate stellar dynamos. Even laboratory analogs do not match stellar dynamos
very well though, with either fluid. See Lathrop and Forest [2011] for a full discussion thereon.
1.1.2 Spherical Couette Flows
The differential rotation in SCF is produced quite different from how a convection-driven
systems. Therefore not all inherent phenomena necessarily have an analog in convection-driven
differential rotation. Imaging fluid flow for SCF experiments is important to geophysical theories
of the core for an additional reason: it is highly relevant to know when the observed phenomena
stem from effects particular to SCF and not necessarily convection-driven dynamos or other MHD
rotating flows. Being able to assess the velocity field can help determine if the flows being generated
are dissimilar to those that can be generated by convection-driven dynamos, and add to the
discussion on the importance of the concomitant magnetic observations during those flows. Here
I focus on the hydrodynamic aspects of SCF.
The mathematical formulation and physical expression for SCF is not scale invariant; there
4
are effective constraints placed on the fluid medium by its composition and state variables. The
geometry and velocity alone are insufficient to fully describe the dynamics; the amount of viscosity
relative to the velocity and length scales also has considerable affect on the flow, as previously
described.
1.2 Parameters for Spherical Couette Flows
1.2.1 Reynolds Number
One way of quantifying the ratio of viscous forces to inertial forces induced by the relative
fluid rotation is the Reynolds number. High Reynolds number flows have weak viscous forces
relative to inertial forces. More functionally, the Reynolds number Re serves as a proxy for the
relative dominance of two terms in the incompressible Navier-Stokes Equations (NSE) description
of flow in Eulerian form:
2∂tv + (v · ∇)v = − 1
ρ0
∇P + ν∆v
∇ · v = 0
(1.1)
where v is the velocity field, P is the pressure, ρo is the density, and ν is the dynamic viscosity. For
simplicity of exposition, the force of gravity has been neglected here, as well as other body forces
due to being on the Earth, and also the forcing at the boundary induced by the moving surface.
Also, this work does not use the fully coupled magnetic equations described, such as described by
Taylor [1963], nor provide a rigorous examination of what creates the particular functional forms
of the flows. Instead, the investigation focuses on reconstructing flows and associated properties
mathematically from acoustics observations, drawing on the physics of the flow only when necessary
to reduce non-uniqueness of results. The focus on High Reynolds number flows is not a restriction
of the acoustic method employed in this work. The methods contained herein work equally well
or better on laminar flows, but they are perhaps more useful for geophysical laboratory studies
related to the core that are difficult to simulate on computers.
The Reynolds number is defined so as to quantify the ratio of (v · ∇)v to the quantity
5
ν∆v in a single parameter, as the ratio of the square of the characteristic velocity associated
with the flow over the characteristic length scale to the dynamic viscosity ν. When the Reynolds
number is very high, and the system is exhibiting a relatively steady mean flow, the viscous term
and the time differential of velocity term become less significant, and the term (v · ∇)v becomes
predominantly balanced by the pressure gradient term - 1ρ0∇P . In a rotating reference frame, the
balance is divided between the advection relative to the reference frame and the Coriolis force;
when the dominant term is the Coriolis force, the flow is said to be geostrophic. The equation is
non-dimensionalized by choosing units so as to scale the characteristic length scale L and speed V
of the flow using x¯ = xL for the spatial variables and v¯ =
v
V for the velocity scaling, and similarly,
setting t¯ = VL t and P¯ =
1
ρoV 2
P . Upon substitution yields
∂tv¯ + (v¯ · ∇x¯)v¯ = −∇x¯P + ν
V L
∆x¯v¯ (1.2)
which can be rewritten as
∂tv¯ + (v¯ · ∇x¯)v¯ = −∇x¯P + 1
Re
∆x¯v¯ (1.3)
Then, presuming the turbulence is not sufficiently unstable that ∆x¯v¯ balances the down-
weighting due to Re−1, the last term 1Re∆x¯v¯ is relatively weak if the Reynolds number is sufficiently
large. Computationally, Re−1 reduces the diffusion of small length scale variation in the velocity
field, and so on a computer 1Re being very small, leads to numerical instabilities.
In numerical schemes for NSE at high Reynolds number, small-scale features over short spa-
tial wavelengths are created that become unstable computationally. As short wavelength features
begin to emerge in the velocity field, high curvatures are created in those locations. Additional
viscosity is added numerically so as to complete a feedback loop whereby as those emerging cur-
vatures are created, the spatial variation in those locations is damped out in the velocity field
evolution. Mathematically, small-scale variations increase the Laplacian where those curvatures
exist, and those regions are diffused spatially by the added viscosity. Thus additional viscosity pre-
6
vents unstable growth of the small-scale features. However, the added viscosity has many negative
consequences as well.
Added numerical viscosity removes the effects on the simulated flow that arise from the
dynamics occurring at shorter length scales that would be present in the actual physical system
being emulated. It also increases the amount of diffusion at larger length scales in the simulation,
smoothing the velocity field even for the mean flow; though certain modern methods try to some-
what alleviate that problem by using spectral methods and controlling the diffusion at different
length scales more directly. Most importantly, because NSE are highly nonlinear, the suppression
of these small-scale features by added viscosity causes poor approximation even of the mean flow
and larger spatial wavelength averages in solutions. Added viscosity to time evolutions of spatially
dependent variables has the effect of low pass filtering the evolving variable spatially at each time
step. For most linear partial differential equations such as the elastic wave equation, constantly
low pass filtering the spatial variations of the evolving variable only reduces accuracy at length
scales being filtered appreciably. That is because the linear differential operators involved admit
a spectral decomposition whereby each spatial wavelength can be evolved separately: a normal
mode decomposition. However, for a nonlinear system such as NSE, evolution of different length
scales cannot be so easily decoupled. The constant low pass filtering of the velocity field dur-
ing the simulation means that the simulation will not be accurate even at the length scales not
directly filtered appreciably from time step to time step. Similar issues also arise from adding
unrealistic magnetic diffusivity, and also the interaction between the fully coupled magnetic and
hydrodynamic equations. These are highly relevant flaws in numerical simulations of geophysical
and astrophysical dynamos. It is critical in such simulations to accurately capture the manner in
which energy is exchanged between components of the system at different length scales.
Laboratory experiments seek to sidestep the unrealistic effects of the artificial damping out
of turbulence and damping out of small-scale features of the velocity field. The length scales over
which turbulence are active is far smaller than the domain discretization or spatial resolutions
employed by most computer models; a problem less pronounced by using analogous physical fluid
7
dynamics that can be generated in the laboratory. Even if sufficient computer memory to store
a model of sufficient resolution to evolve at the length scales of turbulence is available, and suffi-
cient computational resources to overcome the corresponding small time step required are present,
there is still the matter of efficiency. Laboratory data can be collected from very long duration
simulations in real time, as opposed to computer simulations that may take weeks or months or
even years to produce seconds or minutes of simulation for the same regime as the laboratory ex-
periment. Computer simulations have one distinct advantage over laboratory experiments though:
direct access to the physical variables throughout the entire domain at each time step. Taking
advantage of the benefits of laboratory experiments then, requires developing methods to infer
those variables in the manner needed.
Creating flows in the laboratory more analogous to those theorized for core-like flows and
magnetic behavior, requires the use of electrically conducting liquids or plasmas and devices less
amenable to traditional means of inferring the advecting fluid’s velocity field. The kinematics,
temperatures, and vessel designs involved with doing such experiments on practical time scales
make velocimetry more difficult, as does using liquid metal or using any other opaque fluid. This
work presents alternative means of inferring velocity remotely from the boundary of the vessel,
focusing on zonal flows in a core like geometry. For a discussion of modeling the Earth’s core with
liquid sodium flows, refer to Adams et al. [2015]. This work does not address in any detail the
direct relationship to the Earth, its geometry or location in parameter space in regard to various
quantities, or the strength of the analogy between SCF and the Earth’s core.
1.2.2 Ekman Number
For SCF, when observing in a rotating frame, the characteristic velocity scale used in the
Reynolds number is relative to that frame, and another term helps characterize the flow. Noting
the analogy between Vlab and the quantity ΩoL, the Ekman number for SCF is defined by
Ek =
ν
2ΩoL2
(1.4)
8
where Ωo is the rotation rate of the outer sphere, and l is the outer radius Ro. The Ekman number
can be interpreted as the ratio of viscous forces to Coriolis forces. The derivation for it parallels
the derivation for Reynolds number, after substituting the Coriolis force into the left hand side of
the momentum equation in NSE. For sodium experiments, the Ekman number can be less than
10−7; for air during the experiments this work produces inversions from is approximately ·10−6.
For the Earth’s outer core, the Ekman Number is thought to be even as low as 10−14 or less, or
several orders of magnitude higher depending on the value of kinematic viscosity used, one of the
least understood parameters for the outer core; see [Secco, 1995, Olson, 2011, see].
1.2.3 Rossby Number
Another useful parameter for SCF is the Rossby number, which is defined slightly differently
in SCF than in Atmospheric sciences. The Rossby number in SCF is given by
Ro =
Ωi − Ωo
Ωo
(1.5)
where Ωi is the angular velocity of the inner sphere in the SCF. It can be interpreted as the ratio
of the inertial forces to Coriolis forces, the length scale having been obscured due to algebraic
simplification in the expression. The net angular velocity appearing in the numerator is very
important for numericists. It should be noted though, that only a very small portion of the fluid
cavity will follow the inner sphere rotation rate. Even outside that boundary layer, the fluid tapers
off quite quickly in most regimes, and the Rossby number might not adequately characterize the
bulk fluid flow. To within an order of magnitude accuracy, though, it may be said that
Re =
Ro
Ek
(1.6)
The relation follows by noting the numerator of the Rossby number expression is the angular
velocity of the fluid relative to the rotating frame, and can be used to express the velocity scale
in the Reynolds number, and then after cancellation of the length scale factors the relationship
9
follows trivially.
Often it is convenient to work with inverse Rossby number Ro−1. This is because when the
outer sphere is stationary and only the inner is rotating, Ro would be infinite. In the limit that
Ro−1 goes to positive or negative infinity the flow becomes geostrophic in the canonical idealized
sense and the Taylor-Proudman theorem applies throughout the cavity. As inverse Rossby tends
to 0, the Coriolis effects on the flow become negligible and the rotation induced by the inner sphere
dominates as if the outer sphere were having no effect, except perhaps in a thin viscous boundary
layer. The regimes analyzed in this study do not quite approximate either endmember though.
1.2.4 ω Effect
Shear in SCF of MHD experiments is thought to enhance magnetic field gain by converting
poloidal magnetic field components into toroidal magnetic field components. Consider an azimuthal
velocity field given by Ω(r, z). For simple axial dipole dominated magnetic fields, the forcing term
involving the poloidal field in the induction equation for the toroidal component is described
approximately by
∂tT = −(r∂rΩ)(∂zP ) + 1
R
(∆− 1
R2
)T = 0 (1.7)
where ∆ is the Laplace operator, and T is the toroidal component and P is the poloidal component
in the Helmholtz decomposition of the field; the simplicity coming from the fact that the magnetic
field has zero divergence. In this form it can be seen that shear in the rotating fluid by way of
r∂rΩ allows the poloidal component to serve as a forcing term to the toroidal component of the
magnetic field. This is why the zonal flow velocimetry is important to testing many hypotheses
related to magnetic field gain in MHD experiments. To that end,
ωeffect(r) := r∂rΩ(r) (1.8)
10
1.3 Velocimetry of Fluid Flows
1.3.1 Challenges of Traditional Techniques of Velocimetry
In situ observations of experimental and natural fluid flows present many difficulties; they
are disruptive to the flow, acquisition can be difficult or economically prohibitive, and the coverage
of the domain is often either limited or greatly asynchronous.
With such direct measurements, flow diversion and turbulent wakes around the instrument
of observation can disrupt the flow, as well as initiate dynamics in the flow that persist. This
is especially true at scales of experimental fluid dynamics. The invasive nature of simultaneous
direct sampling of the fluid domain with a dense sensor array is discouraging. The alternative
of sampling fewer portions of the domain introduces then issues of aliasing and sampling bias.
More remote, simultaneous, and less disruptive means of interrogation of the fluid are preferable.
Remote observation of fluid flow can be accomplished by several means.
Traditional methods of fluid velocimetry that are less obstructive mainly include optical
techniques such as Particle Image Velocimetry (PIV), Laser Doppler Velocimetry (LDV), and also
ultrasonic techniques such as Ultrasonic Doppler Velocimetry (UDV). PIV measures the displace-
ment between visible particles in the fluid over time to infer the velocity field. LVD and UDV
utilize the Doppler shift the moving particles induce in waves reflected off those particles, optical
waves in the case of LDV and acoustic waves in the case of UDV. Each of the traditional methods
has advantages, but also many disadvantages; in many cases they do not even work.
PIV and LDV methods require that the material not be opaque in the part of the light
spectrum the instrument sees in. Flowing liquid metal is an opaque material, and is especially
difficult to probe using optical techniques, making both PIV and LDV impractical. PIV is arguably
only negligibly disruptive to the flow, but it is certainly invasive. Reliable performance with PIV
requires particles inside the domain matched in density, resistant to corrosion, in the case of liquid
metal resistant to melting, and distributed in a manner conducive to the needed coverage. It is
extremely difficult to match the particles density to that of the fluid, especially when the fluid
11
does not have static state variables, and ultimately over relatively large time scales the particles
either sink, rise, or obey transits in the moving fluid that do not sample the medium uniformly.
Constantly replacing the particles or resetting them is cumbersome and economically discouraging
in terms of labor and time consumption, and even potentially dangerous.
UDV has been shown to work in liquid metal [Takeda, 1991], and does not require addi-
tional particles placed in the medium; [Takeda, 1991] suggests air bubbles as sufficient to explain
the reflections of the sound waves. In liquid sodium, the fluid used in many magnetic related
experiments, oxides floating in the medium are thought to serve as the objects that reflect the
sound waves from which velocity is inferred [Nataf et al., 2008]. In liquid sodium filled SCF, UDV
has been shown to be difficult [Sisan et al., 2004], partially motivating this study. UDV relies
on very high frequency acoustic waves. That means that UDV requires a wetting of the surface
between the instrument and the fluid. Ultrasonic waves can be disrupted by layers of gas forming
between the fluid and the container boundary where the acoustic source is, a layer thicker than
the wavelengths of the ultrasonic sound; without wetting, the waves will tend to reverberate inside
the gas layer more so than transmit through the sodium and back. Also, it is thought that the
tracer oxides may not distribute evenly due to centrifugal and gravitational effects; presumably
the oxides have a slightly different density that would potentially bias them more toward one part
of the cavity than another.
In order to sample larger portions of the domain, UDV LDV and PIV require either large
amounts of instrumentation or instrumentation that can have its target direction modified con-
tinuously. Such modification is challenging for large vigorously rotating setups with mounted
instrumentation. Methods that sample substantial portions of the domain and whose sensors do
not favor a particular directionality or sampling line are preferable.
1.3.2 Acoustic Normal Mode Inverse Spectral Methods
Acoustic normal mode inverse spectral methods offer a means of mathematically recon-
structing fluid flow remotely, requiring only observation at a few locations along the boundary of
12
the medium. They are capable of imaging the entire medium on relatively short time scales. By
virtue of their indirect sampling of the whole medium, they thus also avoid the problems asso-
ciated with spot measurements or sampling bias towards a narrow beam as in UDV. The longer
wavelength nature of low frequency normal modes also allows them to couple to the fluid medium
even in the presence of a gas layer between the container wall and the fluid. This work focuses on
one such method, the adaptation of normal mode seismology, an inverse spectral method, to fluid
velocimetry.
Acoustic inverse spectral methods use the transformation to the spectrum of the acoustic
free oscillations to invert for the fluid velocity, analogously to how normal modes of the Earth
are used in seismology to infer properties of the bulk Earth. Normal mode seismology uses sensor
arrays on the Earth’s surface or in the shallow subsurface to infer bulk properties of the deep Earth
such as density and soundpseed, the excitation sources primarily being large earthquakes. This
form of remote sensing known as helioseismology was adapted by the heliophysics community to
study the interior of the sun [Thompson et al., 2003]. Helioseismology uses the analog of sound
waves in the sun’s interior, acoustic waves whose primary restoring force is pressure—to image the
rotation of the sun’s interior. SCF in the regime where azimuthal flow dominates has an analogous
geometry and flow to the interior rotation of the sun. Both seismic tomography of the Earth and
oceanic tomography often employ acoustic normal mode based inverse theory, but also supplement
their knowledge with other data such as travel times of body waves propagating from a source to
a receiver, to reveal smaller scale features. Like helioseismology, this work only uses normal mode
spectra.
Pressure perturbations from a speaker or source device induce propagation of small ampli-
tude compression waves in specific frequency bands. If the source is sinusoidal and delivers enough
energy at or near the frequencies of modes of free oscillation for the cavity, then resonance occurs
and the amplitude of sound swells at that frequency. That is, when the small acoustic source
device such as a speaker drives a piece of the medium at or near one of the natural frequencies
of vibration of the sounding cavity, large amplitudes of vibration at that natural frequency occur
13
inside. The distribution of power across the acoustic spectrum has localized peaks at the resonance
frequencies. The location of those peaks can be modeled as function of the velocity field, and both
a predictive model for the acoustic spectra based on a flow as well as inversion schemes can be
created therewith. Many aspects of the flow and state variables can shift the frequencies, but
in particular it is the splitting between modes with nonzero azimuthal wavenumbers +|ms| and
−|ms| that is the most informative and least ambiguous effect of the azimuthal flow. This will be
addressed in detail throughout this work.
The acoustic modes of vibration associated with these spectral peaks can be excited in several
manners. Combing through a broad frequency range is possible with a chirp source function such
as the function sin(λt ·ωit), which for sufficiently small lambda delivers an instantaneous sinusoidal
source with frequency ω(t) = λtωi, with ωi being the initial frequency the chirp starts at, and λ
being the chirp rate. In cases when the flow generates its own sound or when other sources of
ambient noise are present, active sources may not be necessary; though I would argue preferable.
Such ambient noise produces less accurate and relatively noisy spectra, since amplitude envelopes
are created across different frequency bands possibly due to time varying excitation or time varying
attenuation. Also, the signal to noise ratio is not as high. For instance, for a simple sinusoidal
amplitude envelope in the time domain, the effect on spectra amounts to convolving with two delta
spikes about the origin in the frequency domain.
Suppose a simple example such as
γ(t) = sin(pit)sin(2pift) (1.9)
then
Fγ(t) = Fsin(pit) ∗ Fsin(2pift) (1.10)
where F denotes the Fourier Transform. Then once looking at |Fγ |2, the convolution has taken the
approximately Dirac impulse that would occur due to a resonance at frequency f and places two
Dirac impulses 1 Hz apart centered on f. Since not every band has the same envelope nor are the
14
bands easily disentangled. This effect cannot be removed so easily. Shortening the time window to
see more short term behavior accrues defocusing as a result of reciprocal spreading in the frequency
domain; i.e. there is an implicit time window function and reduced frequency resolution on top of
any defocusing for purely physical reasons.
In this work, both chirp based data and ambient noise data are used. Chirp-based data,
targeted to excite modes by chirping through a band around the resonance, are used to assess
relevant aspects of the theory and match between the numerical modeling and the laboratory
setup data. Ambient noise based data sets with larger numbers of computed splittings are used to
perform inversions on data.
Small acoustic pressure perturbations to the medium have negligible affect on fluid flow.
Sound waves manifest as compression perturbations in the medium and are approximated as having
no affect on the flows in this work for that reason.
The central physical phenomenon of interest in this work is the velocity distribution of
the fluid inside the acoustically resonant cavity; specifically, the mean flow, the time averaged
azimuthal velocity profiles, presumed to be predominantly azimuthal and axisymmetric in SCF.
As such, often it is more convenient to work with angular velocity, and rotation rates will be spoken
off in place of velocities.
Meridional flows and other inherent features are presumed sufficiently less vigorous in these
flows of interest, and are not included in the modeling of the flow explicitly. Though meridional
flows and other types of poloidal velocity contributions are not explicitly incorporated into the
physical model of the flow for purposes of inversion, discussion of their effects in certain regimes and
relation to how they modify the accuracy of the inversions is considered. In cylindrical coordinates
with z aligned with the axis of rotation, for a flow component in the r,z directions, to first order
the splittings of the modes are not affected. Higher order effects are possible if the flow is vigorous
enough to appreciably distort the mode shapes, and are discussed in later sections of this work.
Here, I will briefly discuss the dominant effect on acoustic splittings, the fine structure of the
acoustic spectrum caused by azimuthal flows.
15
A standing acoustic wave inside the rotating flow within the cavity is transformed by Coriolis
forces and advected due to the transit of the fluid on which it is superimposed. It also is affected
to a lesser degree by centrifugal forces and body forces, as well as slight anomalies in wave speed in
the medium; these are shown to be sufficiently small that they can be neglected. Modes that have
an azimuthal wave number that is non-zero, have a spinning pressure pattern. When the fluid is at
rest, the frequency observed at the boundary comes from an invariant pressure pattern spinning,
the relative amplitudes crossing the sensor with the time evolution. It is important to stress that
it is not the material that spins, but, rather the relative amplitude pattern of the modes that does.
Cylindrically radial motions and lateral motions of the fluid caused by the acoustic wave
are transformed by the Coriolis force, causing the mode to be offset slightly; the eigenvector is
perturbed and also gets a small imaginary component, but the frequencies remain real. Advection
affects splitting due to mode’s spin being relative to the fluid, to first order. Whether the spin is
retrograde or prograde affects whether the spinning mode will spin faster or slower once the fluid
starts moving. If the sensor/microphone is at rest, and the fluid motion advects the spinning mode
by spinning the fluid in the same orientation as the mode, the mode spins faster in terms of the
frequency observed at the sensor/microphone since it is spinning relative to the fluid medium. If the
spin of the mode is opposite the fluid motion, the fluid motion slows down the spin, and lowers the
frequency. As such, the frequencies observed in the lab frame are not the same frequencies observed
in the rotating frame, except nominally for the spherically symmetric modes whose pressure pattern
does not spin.
Since there can exist acoustic standing waves with markedly different modal shapes in the
cavity, their frequency shifts and spectral alterations can differ from one another depending on the
distribution of the fluid velocity; in other words, they have different sensitivities to flow. With
appropriate observations of specific acoustic mode frequencies, and constraints from the theoretical
physics governing the flow, the fluid velocity field can be reconstructed by inversion. Resolution and
accuracy of that reconstruction depends largely on the number, complimentary spectral sensitivity
of the employed modes, and quality and accuracy of the measured spectra. The reconstruction is
16
of the mean flow, where the time integration is on a time scale related to the time window of the
data acquisition and spectral estimations. Within that time window, flow is modeled as having
an invariant mean profile with superimposed transient and possibly turbulent components that do
not fluctuate so vigorously as to significantly warp or displace the spectra from one that would be
produced by the mean flow in their absence.
It is the intent of this study to develop methods and software applicable to future laboratory
experiments, capable of handling irregular geometries associated with actual setups, and to produce
inversions demonstrating the capabilities and limitations with current data. SCF setups in labs
require shafts to spin the inner core, provide stability, hold the inner core in place. This work
abandons the approximation of a spherical annulus/shell, and performs a Finite Element Method
to more accurately incorporate those structures, so as not to distort the flow nor incur undue error
in the forward problem. The formulation would work equally as well for a Taylor Couette flow,
a flow induced by two concentric rotating cylindrical surfaces, or any other azimuthal flow with
symmetry about the axis of rotation. Coupling magnetic equations to NSE formulation of the fluid
flow is not necessary since air is not electrically conducting to any relevant degree.
Posterior error analyses for inversions often state their uncertainties in terms of the model
parameter uncertainty or covariance. In these types of reconstructions though, what is potentially
of the greatest significance is not the amount that the model parameter values would change as a
result of uncertainties in the data. The greatest significance is how accurately the true mean flow
can be reconstructed by the form of parametrization used. In the absence of the shaft and related
asymmetric missing pieces and sensor apparatus, the sensitivity kernels (discussed in chapter 5) are
even about z = 0 for an axisymmetric azimuthal flow. This means that an odd component to a flow
would produce no splitting contribution. Since any axisymmetric azimuthal flow can be written as
the sum of a component even about z = 0 and a component odd about z = 0, and the component
that is odd can be scaled arbitrarily without affecting the acoustics to first order, it is effectively
invisible to the spectra. For this reason, and the fact that the flow is not expected have it’s mean
azimuthal flow deviate from axisymmetry or evenness about z = 0 substantially, this work confines
17
itself to studying the even component of the flow. Slight asymmetries exist in the experimental
setup that do possibly lift that ambiguity, but are small enough that the effects would be conflated
with the noise and uncertainties inherent in the data. A feature that occurs superimposed on
top of the flow in only one of the top or bottom hemispheres, will be averaged equally into both
the top and bottom hemispheres. That must be taken into account when interpreting features in
inversions.
Splitting between the modes with azimuthal wavenumbers +|ms| and −|ms| are linear in
the azimuthal rotation of the fluid, for the regimes examined herein. The splittings can thus be
expressed as a composite effect of the splittings induced by a collection of basis functions whereby
the weighting each basis flow’s induced splitting is weighted in the overall splitting by the same
weighting of the basis function in the synthesis expansion of the flow. A set of splittings can be
predicted from any flow, and an inverse problem can be formulated for the coefficients of the basis
flows based on measured acoustic spectra.
This work primarily concerns itself with large-scale properties of the mean flow. It focuses
on developing methods that can applied when only a limited number of spectral measurements
are avilable, constrained to using only a few sensors that are perhaps not ideally located for more
sophisticated array processing techniques. Testing of methods on data is confined to air, since
sodium data is not currently available.
1.4 60-cm Setup
1.4.1 Laboratory Configuration
This work uses as a case study data collected from a 60 cm diameter titanium vessel filled
with air, the inner core radius approximately 1/3 the outer, the outer spinning at 6 Hz and the
inner core spinning in the range of -4 to -40 Hz. The vessel is described in detail in [Adams, 2016].
The 60-cm setup, Dynamo 3.5 at the University of Maryland Nonlinear Dynamics Labora-
tory, is a SCF setup. It has a diameter to the inner boundary of the outer shell of 60 centimeters.
The dimensions used in the modeling are Ro = 0.3048 and Ri = 0.1000. The outer enclosure is
18
made of titanium, whereas the inner core is copper. The inner sphere is spun by an axle that
rotates with it during the experiments, and the axle is enclosed inside a shaft for support, though
the shaft is welded to the outer shell and spins with the outer sphere. Small pieces of the shaft
covering the axle were removed where it meets the inner sphere. In strictest terms, these additional
components violate the idealization inherent in SCF. The inner sphere must not only be spun, but
must also be held in place. One difference between the 60-cm setup and other experimental SCF
setups, is the shaft is much larger and it spins with the outer sphere rather than with the inner
sphere. The axle spins with the inner sphere, though barely exposed. As can be seen in figure 1.2,
the electronics placed on the boundary for microphones are quite small, especially when the fluid
near the outer shell in those regions is not deviating from the rotation rate of the outer.
19
Figure 1.1: Schematic of the 60-cm setup and configuration in the lab. Diameter of the outer sphere
is approximately 60 cm and of the inner sphere is a approximately 20 cm. Diagram courtesy of
Matthew Adams.
20
Figure 1.2: Interior of the 60cm Setup taken apart showing electronics mounted inside, before
assembly. The setup tries to minimize the invasiveness of the electronics. The axis of rotation,
embodied by the physical axle, is vertical during operation. Photographs courtesy of Matthew
Adams.
21
The housing encasing the structure, as seen in figure 1.1, is related to safety measures for
the sodium experiments. The red “T” indicates the placement of the thermometer in the cavity.
The device was originally used for liquid sodium experiments. The acoustic recordings collected for
this work were performed by Matthew Adams, who also built many of the devices used to retrofit
the device with the acoustic equipment for sound production and recording. For a full discussion
of the setup, refer to Adams [2016].
1.4.2 Modeling the Geometry
The geometry is modeled as being an axisymmetric volume of revolution of the cross section
depicted in figure 1.3. The shaft is approximately 44% the radius of the inner sphere, and has
approximately 75% as large a surface area within the cavity’s interior. The axle is approximately
25% the radius of the inner sphere. The sphericity of the spheres is exceptional, and the shells on
the boundary of the fluid only deviates from being spherical by one part in one thousand, insofar
as the curvatures at welded joints and other aspects are ignored along the shaft [Lathrop, personal
communication].
22
Figure 1.3: Meridional cross section of the geometry, showing the proportions used to model the
cross section. The z-axis is the axis of rotation where the axle and shaft are.
23
Chapter 2: Previous Work in the Field
2.1 Analysis and Discussion
Triana et al. [2014] presented a study of spherical Couette flow using air in a small vessel
whose outer shell’s radius was approximately 15cm, the inner core radius approximately one third
of that, and the inner rotating and speeds up to 36 Hz; the outer sphere stationary. Since the
outer shell was stationary, anemometer readings could be taken easily. Flow profiles taken by
the anemometer showed that the flows were quite similar at the two elevations measured, which
can be interpreted to suggest significant vertical invariance of the flow. Triana et al. [2014] also
attempted an inversion of acoustic mode frequencies for the flows. The inversions utilized damping
to regularize the inverse problem, since for the flow parametrization used, it was under-determined.
While many aspects of the inversion procedure were straightforward, the nuance of applying it to
laboratory SCF merits some discussion.
In particular, methods employed in the self described “proof of concept” study carried out
in Triana et al. [2014] could not yield an unambiguous choice of regularization parameters, which
were ultimately chosen based on the anemometer data [Triana, Personal Communication]. In this
section I will discuss that work, which is to the best of this Authors knowledge, the only other
published work in experimental SCF using helioseismic methods. To do so, I will first recast the
inversion method formulation in a different form, so as to foster a theoretical analysis rather than
a computational approach.
Consider the traditional weighted least squares problem, cast in an optimization form with
objective function Φ. Let C be the forward operator, d be the observed data vector with uncer-
24
tainties modeled by σ, and m the model vector of fitting parameters over which to minimize
Φ(m) =
1
2
‖Wd(Cm− d)‖2, (2.1)
where Φ(m) is called the data misfit in this context and Wd is a weighting of the data by some form
of inverse data standard deviation matrix, or confidence measure based the data uncertainties. In
the case of uncorrelated errors, a diagonal matrix with reciprocals of the uncertainties of each
entry in the data vector d is most typical, which was what was used in Triana et al. [2014]. Wd is
analogous to an inverse data covariance matrix for data uncertainties that are standard deviations
of a sampling of a random process. To prevent unphysical solutions and lift the ambiguity of the
answer owing to the under-determined nature of the problem, regularization was used.
Regularization by damping is a method that seeks to minimize the Φ in the presence of a
penalty function ‖D(m)‖2 as follows:
Φ(m) =
1
2
‖Wd(Cm− d)‖2 + ‖D(m)‖2, (2.2)
where  controls the trade-off between data misfit and regularization, and D is typically some
linear function of the flow such as a differential operator like the Laplacian or second derivative
with respect to some spatial coordinate(s). The trade-off arises from the relative minimization of
the two terms in the expression for the objective function Φ; the first term is often called the data
misfit. As can be seen from Eq. 2.2 the choice of  is only scale invariant with respect to the flow
magnitude when Wd is not a function of the scale of the flow; i.e. when the uncertainties do not
change as the magnitude of the flow or flow pattern changes, possibly due to changing the rate of
the inner sphere. Though it cannot be judged from the composite uncertainties presented in Triana
et al. [2014], examination of the actual data collected in the study shows that the uncertainties
are not invariant as a function of the flow. It is implicitly understood that if the number of data
observations changes, that  must be scaled accordingly as well, else the trade-off does not have
the same proportional weighting. More importantly, how  is chosen is critical. To motivate that
25
discussion, first I will address other features of the formulation.
It can be seen that the damping in Eq. 2.2 is not done point-wise or in any localized sense;
it is expressed as a norm, an integral over the volume of the cavity in the case of Triana et al.
[2014]. As a result, different choices of  will not only cause trade-off with data misfit, but also
cause the degree of regularity in the domain to trade-off from one region to another. That behavior
is benefitial for a large data set, since localizing in a region at the expense of increasing penalty
function term trade-offs with the localized increased curvature’s ability to decrease the data misfit.
With small data sets that do not have good coverage across the whole domain, this behavior results
is ambiguity due to many possible ways of localizing high curvatures that fit the data. Ambiguity
is created when solutions produce similar values of the objective function even with the included
penalty term. Many methods attempt to lift the subjectivity of regularization parameters that lift
the nonuniqueness, such as optimally localized averages and minimization of functionals defined on
the resolution matrix of the inversion. Triana et al. [2014] did not use such functional minimization
explicitly, but suggested optimally localized averages as a form of improvement. I will address both
the method used in Triana et al. [2014], and also why there are similar difficulties encountered with
that suggested improvement as well.
In the case of the Tikhonov regularization, Triana et al. [2014] chose to damp curvatures
in the spherically radial direction and the polar angle direction. Various metrics were used as
motivation for the choice. First, it is important to be aware of certain aspects of inverse spectral
methods, to put those metrics in context.
Inverse spectral methods are quite dissimilar from fitting surfaces or curves to control points
or direct measurements of the profile being reconstructed. In the case of these types of acoustic
methods, how well the inversion is able to predict the acoustics data is quite different from how
well it is able to predict the flow. If measurements of the flow at points inside the cavity are used,
then enforcing smoothness on the solution amounts to preventing oscillation in the reconstructed
profile; the inversion is thus interpolating between data points measured inside the flow while
allowing tolerance to the fit at measured data points when it reduces the curvature of the solution
26
overall, based on some metric to judge the trade-off between the two ideals of fitting the data
enforcing smoothness. When the data measured is indirect, as in acoustic spectral methods, and
not measurements of the variable that is being reconstructed, problems can occur.
In theory, with sufficient data, the quantities directly measured and the properties recon-
structed by inversion are in a strong correspondence and this subtle issue is of little concern.
Helioseismology typically uses several hundred measured modes to obtain results considered accu-
rate in the outer convective zone, and arguably in the upper most portion of the radiative interior,
losing resolution closer to the poles. In the experimental 30cm setup, those analogous regions of the
cavity the method is best suited for were nearly at rest, and all interesting features were primarily
localized within the tangent cylinder. Furthermore, fewer than 30 modes were used. Therefore,
the problem tackled by Triana et al. [2014] is not an easy inversion to perform, regardless of the
method or approach.
Many disciplines use a prior model of the quantity being inverted for, using damping of
curvature via covariant differentiation relative the background model; extreme damping then just
reproduces the background prior model. In this case, solutions trade-off between endmembers.
One endmember, corresponding to no damping, perfectly reproduces the data at the expense of
being unphysical or fitting the noise. Another endmember, corresponding to maximum damping,
perfectly reproduces the background prior model by damping out all deviation therefrom. Damping
in several variables can be formulated as a single variable of damping, with a relative weighting
between the components of the damping operator being another tuning parameter.
Triana et al. [2014] used a Tikhonov regularization method that damped the second deriva-
tives of the spatial coordinates 1) the spherical radius R, and 2) the polar angle θ. Damping in
the spherically radial R and polar angle second derivatives, means that latter endmember is a flow
linear in R and linear along each polar latitudinal curve. The first linearity in R is very unrealistic
at the high Reynolds number in the flows inside the 30 cm setup during the experiments. The
second linearity in polar angle is somewhat bizarre in that context, though the constraint of the
evenness of the solution about z=0 imposed in the study modifies that constraint to tend toward
27
being invariant along polar angles θ as the damping tends to infinity; again, that means the so-
lution is spherically symmetric though, which is quite unrealistic unless the viscosity was many
order of magnitude higher. Such a regularization scheme is not a problem when there is a plethora
of data, in which circumstance the damping need only enforce continuity or some small amount
of regularity in the reconstruction. However, it does pose problems when using so few observed
modes.
High shear zones, which are present in MHD studies, present another problem in that
context. Regions of high shear in the flow will be defocused when damping is based on curvatures
because high shear zones have sharp curvatures; it is unavoidable in the absence of choosing a
biased parametrization or biased damping strategy. With a small data set, there is not enough
resolving power that a smoothed and reduced shear in the solution would significantly degrade
data misfit. Furthermore, the extent to which the shear zone is smoothed out also trades off with
smoothness in the rest of the domain. If a particular mode is much more sensitive to a region
than the rest, the solution will tend to produce a flow that poorly fits that mode’s splitting, and
then correct the issue by placing larger flow in the region(s) it is relatively uniquely sensitive. The
expense to the misfit of the other mode spectra is negligible in such a case. Regularization can
only prevent such behavior when it smooths the solution sufficiently that it cannot oscillate over
length scales as large as the smallest region of reduced complimentary sensitivity in the modes
used.
So, with small amounts of observed modes and non-ideal resolving capabilities, there is no
way to preserve shear unless additional biasing assumptions are made. Since the penalty function
in Eq. 2.2 is an integral in this context, it will then trade-off with curvature in other regions of the
domain. With limited data, that is difficult to overcome and reconstruct the shear correctly. The
variables chosen, the spherical radius and polar angle variables, also penalize a geostrophic flow or
one that is somewhat vertically invariant more so than a spherically symmetric flow linear in the
radial and/or polar variable. For spherical Couette flows well beyond the laminar regime, a strong
component that is vertically invariant is expected, especially outside the tangent cylinder and away
28
from any jet instabilities after the velocity field is time averaged. The anemometer data in Triana
et al. [2014], even though it only sampled limited locations, also supports that contention.
It would be perhaps better to use the cylindrical coordinates r and z for curvature damping
when the flow has a strong vertically invariant component, applying light r curvature damping,
and heavy z curvature damping. Either way, it raises the issue of whether a feature is being
created or reconstructed. For example, z curvature damping may create cylindrical symmetry in
the reconstruction when it is less pronounced in the actual flow, whereas in the case of spherical
parameters used in R, θ curvature damping may erode the cylindrical symmetry in the inversion
despite cylindrical symmetry being dominant in the actual flow. With only several dozen splittings
for an inversion, one should not expect to obtain similar inversions using both methods without
additional information or severe loss of resolution, because with scant data and gaps in coverage,
the damping must do more than enforce continuity. The damping must constrain the solution by
supplementing the data to a considerable degree. The ultimate difficulty with constraining the
solution by damping in the context of very limited data, is the choice of damping parameters is
largely subjective even when more advanced means are used to select them.
Ultimately, damping parameters for the Tikhonov regularization in Triana et al. [2014] were
chosen so as to minimize the misfit of the inversions with the anemometer profiles of the flow
[Trianna, Personal communication], and the motivating analysis did not yield an unambiguous
justification for the damping parameters. As a result, the inversion fit the profiles well even inside
the tangent cylinder where the mode set chosen had relatively little sensitivity. Figure 2.1 shows
the result for two inversions, the Tikhonov and semi-spectral Bayesian approaches; it also shows
what the anemometer data inferred as the rotation profile.
29
Figure 2.1: Reproduction of a figure from Triana et al. [2014], showing anemometer profile and level
of consistency with inversion profiles, as well as posterior uncertainty estimates. s is cylindrical
radius, ro is the outer sphere radius.
30
If the inversion had been done as if the profile was unknown, and still recovered these profiles,
that would suggest robustness of the choice of inversion parameters and perhaps an over estimation
of the levels of uncertainty inherent in the shaded region in figure 2.1 expressing the posterior
error bounds. However, since damping parameters were chosen minimizing RMS error along these
profiles, it is actually somewhat expected that this result is obtained. In fact, it is actually
somewhat surprising it doesn’t fit the anemometer data better. I will address why it doesn’t fit
better shortly. The motivation for modal based acoustic velocimetry in MHD experiments though,
is to be able to infer velocity without a tool such as anemometers or UDV; as such means are often
not available. Even if an experiment can supplement its acoustic observations with traditional
UDV in a manner analogous to the anemometer recordings, other issues arise, and are discussed
below.
As can be seen from figure 2.1, even from that choice of damping parameters the posterior
error distribution shows a large distribution of velocity profiles consistent with the data under the
stated uncertainties, though the covariance cannot be gauged from the plot. The bounds of the
distribution is so large that it does not fit inside the plot they presented. One possible reason for
this, is that minimizing the root mean square error along the anemometer profile with respect to
cylindrical radius r placed equal weighting at each cylindrical radius r despite different r values
being associated with different volume of fluid, and the modes used being much more sensitive to
flow outside and away from the tangent cylinder. An inversion performed this way should fit the
data better where the modes have more resolving power, and the data will be better predicted
by the reconstructed flow. Personal investigation with this mode set found that reconstruction
outside the tangent cylinder was somewhat robust to changes in damping parameters, but inside
the tangent cylinder could produce very large changes in the profiles relative to very small changes
in the damping parameter. If the flow is not approximately vertically invariant then, even if UDV
measurements were used in an analogous manner as the anemometer data there is no guarantee
it will accurately reconstruct the flow anywhere besides along the UDV profile. On the other
hand, if the flow is presumed to be relatively vertically invariant, and one has a UDV profile
31
along a cylindrical radius, then the importance of an acoustic inversion is greatly diminished. The
reconstructed profiles in figure 2.1 are thus at best a test of consistency.
The work also showed other criteria like averaging kernels and error magnification. Averag-
ing kernels are rows of the resolution matrix, though as a function of space rather than arranged
as a vector stack’s dual. They quantify the diffusion from basis function to basis function as a
result of the inversion operator mapping from the model space to itself. More specifically they
quantify the diffusion into a basis element from the others by giving the relative weighting of flow
in the region of space contributing the reconstructed flow. The reconstructed model parameter is
a weighted average of the flow, which if expressible in the parametrization associated to the model
parameters, amounts to a weighted average of the model parameters describing the flow. When
the model parameters are the values of flow at grid points or a fine tiling of the domain, then
m¯rozo ≈
∫
rz
Ko(r, z)mrz (2.3)
where m¯rozo is the recovered estimate for the model parameter mrozo , and the integral is carried
out over the entire cross section after the volume integral has been reduced by axisymmetry of the
parametrized flow to an area integral.
Averaging kernels are most useful when to some extent the averaging kernel is a slightly
defocused form of a Dirac delta distribution, so that the inversion only blurs flows slightly and
in a localized manner. That is, the value of flow produced by an inversion at a point or element
can be interpreted as the average in the neighborhood around it, usually some type of distribution
for a weighting function localized around the point or element. In such a case, the inversion can
be understood a the flow that would result from blurring the true flow by a low pass filtering
operation. The averaging kernels were quite delocalized and irregular in Triana et al. [2014], as
can be seen by a reproduction of a figure from that study in figure 2.2. And in motivating the
choice of damping parameters, chose one particular small region that was less irregular looking, ; it
did not take into account the rest of the domain in the analysis. Averaging kernels are more useful
when the coverage of the domain is more uniform, yielding a measure of defocus that damping
32
induces. Using one location that is less defocused is perhaps, at best, useful for studying that
region in the cavity; though it is not actually sufficient even for studying one region in the cavity
either, since that region will not be affected by the inversion’s diffusion into that grid point or tile
the region is nominally associated with, but also by diffusion to and from the points contained in
that region, possibly from very nonlocal sources. See figure 2.2 and note the features.
Figure 2.2: Reproduction of a figure from Triana et al. [2014], averaging kernels for several locations
used in the inversion. The white cross hairs are showing the location in the domain, which has been
broken into thousands of small tiles as basis functions, and the colors are showing the averaging
kernel for it. Note the high non-local diffusion, especially for the bottom left sub-figure. The
reconstructed rotation rate at the point marked, will be a weighted average of the rotation rates
at other locations in the domain, and the weights are given by the coloring of the domain.
33
The bottom left averaging kernel in figure 2.2 is quite problematic in terms of resolutionand
exaggerates the difficulty with even defining a width for the distribution associated with the point
at the white cross hairs. These observations suggest that using a set of basis functions that are a
fine tiling of the domain, such as in Triana et al. [2014] that used many thousands of basis functions
tiling the (R,θ) plane in the quadrant it parametrized, is overly ambitious in terms of resolving
power. Indeed, the damping has the effect of removing any small-scale features. Assuming the
flow is even about z=0, and having n data points, one should not expect to be able to have a
solution that has more resolution than n controls points in the upper half of a meridional cross
section could produce an interpolant for.
The situation is somewhat worse actually. The coverage of the spectrum that a mode set
used to perform inversions with has, is altogether different from the coverage of the domain that it
can resolve well. Many modes have similar sensitivities. It must be stressed that having n modes
spread out in terms of the acoustic spectrum, whose splitting has been measured, is not analogous
to having a uniform grid of n control points that in situ measurements of the domain yielded.
It is perhaps best understood in terms of the analogy of fitting a curve to a set of data points.
Similar sensitivities to flow are analogous to data points being very close to each other, such that
they do not so much add spatial information, but, rather, simply reduce the uncertainty at the
spatial location that they mutually inform about. Even using a basis that has fewer functions
in it than the number of data point, traditional weighted least squares inversion cannot always
be performed by an explicit inverse, or at least inverting matrices so ill-conditioned that direct
inversion is highly inadvisable, and still requires some form of regularization or truncated singular
value decomposition to invert for the model parameters. Many methods try to foster forms of
objectivity for choice of damping parameters.
Methods such as minimizing the Dirichlet spread function or the traditional Backus and
Gilbert approach of optimally localized averages are popular choices for regularization. They
actually try to choose the damping parameters so as to minimize the diffusion/blurring caused by
inversion in an ill-posed inverse problem. Though these were not used in Triana et al. [2014], since
34
they were recommended as an improvement, and are commonly used in both deep earth seismic
tomography as well as helioseismology, they warrant comment. Averaging Kernels as depicted in
figure 2.2 would normally dissuade one from exploring that option, however the depicted averaging
kernels were for parameters chosen to fit the anemometer rather than what had the best optimally
localized averages. If a resolution and damping parameter combination is chosen to allow a fit
along the anemometers, and that flow is not well resolvable with the data and mode choice, it is
to be expected that the averaging kernels would look bizarre. For that reason and other matters
to be discussed, the approach seemed like a possibly viable option, especially if the expectation of
resolution were decreased and a less fine basis were used.
The Backus and Gilbert approach essentially tries to localize diffusion to each region of
the cavity during inversion, using a basis that tiles the domain with small localized functions
such as a Cartesian or spherical grid. The Dirichlet approach tries a similar method, but only
penalizes the diffusion itself, and does not weight differently diffusion between distant regions and
diffusion between neighboring regions. As such, the Dirichlet approach tends to preserve smaller
scale features better, but as a result also produces artifacts, such as overshoot. Both methods
have noble goals. Unfortunately, though they proceed by optimizing each respective property just
described, with small amounts of data and domain coverage there is a wide chasm between doing
something optimally and doing something well. Also, they are not without a certain amount of
subjectivity.
There is little merit in preventing defocus in many regions in the cavity. For instance, the
region that the averaging kernel selected in Triana et al. [2014] as motivation for its inversion, is
not one where you would expect localized variations in flow. Even if that averaging kernel had
ended up lifting the ambiguity of the choice in damping parameters, it would have done so at the
expense of the regions in the cavity where there are possibly localized variations in flow. Such a
region was examined in that study, and not used due to how nonlocalized the diffusion was.
If the inversion is carried out on a 2D slice and assumes axisymmetry, then to weight each
volume element of the cavity equally requires weighting the damping functions by the cylindrical
35
radius; the integrals are done on a 2D slice and the true domain is the volume of revolution obtained
from it: dV = 2pirdA. As a result, the trade-off in defocus optimization by volume fraction is
biased toward regions of the cavity where there will be less variation; there is no reason to prevent
blurring of somethign that is already relatively homogenized. Choosing a weighting function then
that is optimal is difficult if features present in the inversions are to be deemed recovered features
rather than enforced features. Some approaches to Backus and Gilbert also impose additional
constraint equations that try to enforce conservation of the reconstructed variable. That is, if a
basis function is diffused by the resolution matrix, that diffusion should be conservative in some
manner. The difficulty with SCF though, is that the rotation is not a conserved quantity. It
is tempting perhaps to try and enforce some type of conservation of angular momentum, but as
the number of choices increases in the inversion, the solution the Backus and Gilbert approach
provides becomes less objective. Many other disciplines have the benefit of heuristically adapting
the tunable parameters, and such methods flourish.
Triana et al. [2014] also tried to show consistency with the error magnification operator. An
error magnification operator essentially gauges a similar idea for the data as the averaging kernels
do for the model parameters. It gauges whether the inversion process inflates errors. An inversion
operator that overfits the data given the data uncertainties, will have an error magnification less
than 1, and an inversion that does not predict the flow within data uncertainties will have an error
magnification greater than 1; on average. Error magnification operators though, are more complex
when the errors are not invariant from data measurement to data measurement.
Triana et al. [2014] simplified the analysis by assumed uniform errors and defended the
choice by asserting that the errors were all roughly equal, which was not the case nor was that
assumption that actually needed to be approximately true in order to make the mathematical
manipulation made. To factor out the error magnification operator from a weighted sum, the work
claimed that the errors in each splitting were all roughly equal; in which case an operator that
varied with each component of the sum coudl be moved outside the sum. They actually ranged
between 2 and 35 and were on the order of 10% often. The factorization actually requires not
36
the errors to be all roughly equal, but for the square of the errors to be roughly equal, i.e. 4 to
be roughly 1225, which is clearly not the case in that context. The assumption may not be so
severe in its effect though, since Triana et al. [2014] used composite data from different rotation
rates normalizing by the inner sphere. As such, its error weights that stemmed from standard
deviations are actually conflating the standard deviation with some form of standard error of the
mean possibly. If the inversion is for the mean state of flow composited, the uncertainty is much
smaller than what was used; treating them as standard deviations implicitly asserts that the same
thing was being sampled unbiased and the difference was due to random errors. It is not clear to
this author what the benefit is of using composite data from different rotation rates. It should be
noted though, that average error magnification is perhaps not a good metric, especially when there
is modeling uncertainty due to neglecting the shaft/axle. Many modes are often quite easy to fit
very well, in terms of their splittings measured, by a plethora of qualitatively different solutions
whose differences only contribute negligibly to the overall misfit, whereas some modes are quite
hard to fit simultaneously, and produce larger errors. Error magnification being 1 on average is
not only somewhat arbitrary, but does not assess the distribution of misfits.
Similar issues relate to the ”semi-spectral” method employed in Triana et al. [2014] as
well. In the ”semi-spectral” method employed in Triana et al. [2014], the damping parameters
chosen were justified as noting the parameters chosen ”obtain a model that provides a good fit
to the measured mode splittings”; but in that type of inversion, there should be an increase the
goodness of fit as the tunable regularization parameters are decreased. Fitting the data adequately
is insufficient justification. Changing the damping parameters even slightly in a method such as
this with few modes, can produce very large differences in the recovered profiles, especially inside or
near the tangent cylinder. There are no unambiguous criteria for how to choose these parameters
presented. Subjective inversion parameter choices are not ideal, especially when the test case is
one for which the profile is reasonably well known as in this case, or anemometer profiles exist.
The most obvious question, is if the parameters in the inversion were chosen subjectively, the
anemometer profile was used to select the damping parameters, then why doesn’t it fit much
37
better along the anemometer profiles? There seems to be some possible causes.
The splittings presented graphically in Triana et al. [2014] as a function of rotation rate of
the inner sphere,upon my examination were not representative of the rest of the splitting data.
There are significant deviations from linearity that cannot be readily explained by a less precise or
changing motor rotation rate variability, which is visible in the acoustic spectrum. Even excluding
the lower rotation rate data, the higher rotation rate splittings do not yield a (Ωi, δf) curve that
approximately goes through the origin, even with creative approaches to data weighting based on
the listed uncertainties. That is to be expected if the flow is changing qualitatively. Splitting
from the m=0 mode is to be expected due to the shaft temperature drift or poloidal flow. An
extrapolated inferred splitting of the +|ms| from the −|ms|mode at zero rotation rate, the quantity
contained in the data set examined, for m nonzero, is not as easily interpretable, though it would
have been be easy to verify directly. Triana et al. [2014] explains this by noting the shell is
not perfectly spherical, but slightly prolate. With fluid at rest, though it is quite obvious that
+|ms| should not split from −|ms| due to asphericities that are axisymmetric about the shaft
since they satisfy the same exact PDE, it is perhaps not as obvious that asphericities that are not
axisymmetric cannot do it either. Just as in the perfectly spherical case, it is expected that the
+|ms| and −|ms| pair should be time reversals of each other. Though it cannot be said that they
spin in opposite directions exactly due to the non-ideal geometry, they should still evolve as mirror
images with respect to each other with respect to the time axis, and hence have the same observed
frequency. The only thing that seems a plausible explanation for them being split at rest would
be a body force placed on the system, such as the Earth’s Coriolis force. Gravity and centrifugal
effects are two weak in this regime to produce the splitting. Perhaps some other factor involving
the mechanical device that spun the shaft is not included in the model?
The simplest explanation, assuming splitting was not measured between mismatched modes,
is that the flow is changing, and despite being linear about some state, is not linear about the fluid
being at rest. The only other possibility is that the splitting was measured between modes with
different values of |ms|, which would be offset at zero by the distance between rest frequencies;
38
that seems unlikely given the description of the method used. Though perturbation theory can
be made to work by linearizing about another state than the fluid at rest, the kernels used for
inversion were those about the fluid at rest. It is not likely that poloidal flow components caused
the shift, but even if they had, the linearity about the higher rotation rate would still not resolve
the issue and poloidal related kernels would need used. Asphericities may have diverted flow into
the (r,z) directions, but that would only affect center frequencies appreciably, not splitting unless
the poloidal flow was so vigorous as to distort the mode shape; that seems unlikely.
The forward problem of perturbation theory does not require linearity in the inner sphere
rotation rate. It only requires the splittings are linear in the mathematical scaling of the flow from
some reference state. The method used in Triana et al. [2014] required more than the first order
perturbation theory regime, it required the flow be linear with respect to the inner sphere’s rotation
rate, else the data at different rates could not be composite by normalizing the splittings by the
inner sphere rotation rate. The kernels used implied that linearity is the form that mandates that
zero splitting occurs at rest, with no affine shift corresponding to an offset at zero rotation rate
of the inner sphere. The inversion procedures ignored the shifts and only used the slopes. To be
clear, it is not just ignoring the offset due to splitting from the ms = 0 modes as a result of the
axle, but also ignoring the offset of the trend in splitting between +|ms| and −|ms| at 0 rotation
rate of the inner sphere, as discussed below.
Regardless of whether the flow actually scaled with inner rotation rate, it is not clear what
is to be gained from using multiple rotation rates to perform an inversion, and the variation in
the data suggests different rotation rates would yield a substantially different answer when used
individually.
In summary, it is not clear to this author that the data and anemometer do actually support
the forward problem being established in Triana et al. [2014]; at least not established in a manner
that would allow splittings measured at different rotation rates to be composite as they were.
Furthermore, since data was not acquired for the low rotation rates on many modes, which means
that it is not possible to track their splittings back to the zero rotation rate to ascertain if the
39
linear relationship between inner rotation rate and splitting was robust. The intercept on many
was several Hz. (n,l,m)=(0,4,1) had an intercept of 3.55 Hz with an error estimate of 0.18 Hz,
but at most climbed to a splitting of 9.46 Hz. Refer to section 3.5 for a discussion of (n,l,m).
Many had intercepts 1, 2, and 3 Hz as intercepts. (0,13,5) had an intercept of over 8 Hz with an
uncertainty of less than 1 Hz. Using composite splitting data then, means the flow goes to rest at
0 Hz rotation, yet there is a body force on the fluid splitting the modes possibly. It is not clear
such offsetting affects could be large enough to cause the observed discrepancies.
2.2 Contrast with this Work
The Experimental setup for which most of of the acoustic data used in this thesis uses
has a much larger shaft in terms of cylindrical radius, though has much smaller deviations from
spherical symmetry otherwise. If the offsets in the linear trends in splitting were somehow due
to lack of axisymmetry, they should not be expected for the 60cm setup. Such geometry related
concerns support finite element methods over approximations of the geometry as being spherical,
as acknowledged and recommended in Triana et al. [2014].
The benefits of FEM are more than matching predicted frequencies of the modes of the
stationary medium and preventing distortion to the inversions. They also improve the forward
model’s linearity. Even very small perturbations to linear systems can require higher order correc-
tions when too many idealizations are used in the modeling. For instance, no matter how accurate
the first order correction to the spectrum can be predicted for a rotating sphere of fluid, if the
actual geometry is elliptical, then not only would one have to add a correction for ellipticity, but
if severe enough, also add a second order correction. The second order correction to the rotational
splittings taking into account the coupling between the rotational affects and the affects of ellip-
ticity; the first order correction due to rotation does not make a good linear forward model, since
it was not computed in a geometry that well approximates the experimental setup.
The microphone configuration used in Triana et al. [2014] is far superior to the setup in
the 60cm setup though. It had 4 microphones equally spaced in the azimuthal direction, which
40
allowed a much more robust array method to assess the magnitude of the azimuthal wave number
|ms|. The measured spectra were also made using a loud chirp, and the audio files are relatively
clean, unlike the 60-cm data.
When only the inner sphere is rotating, one can assume a fluid rotation profile that is sign
invariant, since rotating the inner sphere in one direction presumably does not produce rotation
in the opposite direction of the fluid, at least on larger length scales when time averaged. The
Kernels have a single sign across the domain, and hence the splittings should all have uniform sign.
That is not necessarily true in a rotating frame where the acoustics are recorded. For flows near
uniform rotation for instance, the splittings are more so dominated by Coriolis forces rather than
advection, when measured in the rotating frame where the fluid is nearly at rest throughout the
bulk of the cavity.
The largest difficulty with the 60cm setup, compared to the 30cm setup in Triana et al.
[2014], is that the 60cm data has no empirical validation for the velocity fields inside. The data
has no traditional velocimetry measurements to assess the accuracy of inversions are available.
That is not altogether different from the situation in deep Earth seismology or Helioseismology,
but given the relatively accessible nature of the 60cm setup, it is unfortunate. The difficulties with
getting direct measurements of velocity is actually a motivation for this study. The developed
methods are in response to that difficulties. That difficulty coupled with the small size of the data
sets, has motivated a considerable amount of conservativeness to the approach, forsaking many
nuances that could possibly be teased out of the data, in favor of producing the most reliable and
robust results. In chapter 8, a modest data set that is available is used, and these types of issues
are addressed.
The methods employed in Triana et al. attempted to utilize the tools of helioseismology and
deep earth seismology for laboratory velocimetery, but the methods employed in those field utilize
a great number more modes and observations and use prior models of the expected structure to
linearize from. These reservations expressed regarding Triana et al. [2014] are meant to reflect the
difficulties with acoustic inversions of this sort, and motivate other possibilities that are explored
41
in this thesis. The most obvious solution might appear to be to fit the cavity with an array capable
of automating the identification of hundreds of modes. While that sounds quite feasible for studies
in air such as this work and Triana et al. [2014], since the actual goal is to be able to adapt the
techniques to pre-existing liquid sodium experiments and similar devices, this work follows the
mindset of extracting as much information robustly as possible from minimal data in chapter 8.
The development of the tools for calculating the forward operator is amenable to any number of
modes and many other methods of inversion.
42
Chapter 3: Mathematical Formulation
The mathematical formulation reduces the general 3D problem to a 2D problem by exploiting
the axisymmetric nature of the domain. This allows greater accuracy at a small fraction of the
computational cost of performing this operation in 3D. Accurately computing the modes of the
stationary media is importanct since they can be used to compute the first order effects of rotating
flows on the spectra.
3.1 Acoustic Pressure Perturbation in Stationary Homogeneous Media
In terms of material displacement, the acoustic Helmholtz equation for sound waves in a
stationary homogeneous fluid has the form
ω2u = c2∆u, x ∈ V
u · n = 0, x ∈ ∂V
(3.1)
where u is the material displacement vector, ω is the angular frequency of acoustic oscil-
lation, c is the speed of sound in the cavity, the domain V with boundary ∂V , and n is the unit
normal to the boundary. The boundary condition idealizes the boundary as being perfectly rigid
and an infinite mechanical impedance contrast. This work focuses on the cases where that is a
realistic approximation, and does not include an impedance term or coupling to the elasticity of
the bounding material.
By taking the divergence of the above equations, one finds the formulation in terms of
pressure perturbation
43
−ω2p = c2∆p, x ∈ V
∂np = 0, x ∈ ∂V
(3.2)
where p is the perturbation in pressure relative to the background pressure, and ∂n is the normal
derivative at the boundary wall. A pressure formulation of the Helmholtz equation allows the
problem to be expressed as a scalar partial differential equation, thereby reducing computational
requirements. Inclusion of a boundary term does not seem warranted for experiments in air, but
for other applications the inclusion of a boundary term ∂p∂n =
iρω
Z may be warranted, where Z is the
reciprocal of the field admittance at the boundary. Treating the housing of the fluid as a coupled
elastic structure could be warranted in some situations involving liquid sodium and certain types
of containers, but is not addressed herein. With the vanishing Neumann boundary conditions the
Laplacian has a discrete spectrum on this type of domain.
3.2 Exploiting Axisymmetric Geometry for Computational Efficiency
The geometry of the cavity is exceptionally axisymmetric, and is so modeled. The axisym-
metry of the domain enables a partial separability of the Helmholtz equation—in the azimuthal
variable. The 3D Helmholtz equation in pressure for each azimuthal number m is mapped to an
equivalent elliptic partial differential equation in 2D axisymmetric geometry. This is done by using
cylindrical coordinates, and an auxiliary wave function ψ(r, z) satisfying
p(r, φ, z) = ψ(r, z)eimφ (3.3)
Thus for each azimuthal number ms,
∂
∂φ = ims, using units so that c is uniform and unity,
the Helmholtz equation reads
−ω2ψ = ∆rzψ − m
2
r2
ψ
∂nψ = 0
(3.4)
This PDE can interpreted as the pressure profile at φ = 0. Now the problem has been
44
transformed from being 3D vector valued on a 3D domain, into being scalar valued on a 2D
domain. In this study, only values of ms = 0, 1, ...5 were used. The identification of the azimuthal
number ms for a mode is thus performed apriori, and is part of the formulation.
3.3 Variational Form of the Acoustic Helmholtz PDE in Pressure
The computational benefit acquired by the reduction to a scalar PDE on a 2D domain also
results in a form that is not a standard wave equation. I next place it in a more standard elliptic
form:
∇rz · (r∇rzψ)− m
2
r
ψ = −ω2rψ
∂nψ = 0,
(3.5)
For an admissible test function ζ defined on the domain V , multiplying both sides by ζ and
integrating by parts yields
−
∫
A
r∇rzψ · ∇rzζ dA−
∫
A
m2
r
ψζ dA = −ω2
∫
rψζ dA∫
∂A
r∂nψ · ζ dA = 0,
(3.6)
where the first integrals are carried on the area cross section A on which this 2D PDE is defined,
the last being carried out on ∂A, the boundary of A . The absent boundary term is due to the
vanishing Neumann condition ∂nψ = 0 along the boundary of the domain.
3.4 Discussion of Variational Scalar Formulation
The reduction of the problem to 2D, as done in the previous section, has enormous benefits
that are difficult to overstate. Having solved the problem in 3D in Comsol Multiphysics, I found
obvious drawbacks worth mentioning.
In 3D the number of elements needed, even when using higher order basis functions, is
much larger. Though it is still manageable, many of the tasks involving changing the boundary
geometry over a range of values to assess the kind of uncertainties associated with their tolerances
45
would have taken much longer. Also, it became a chore finding the correct solver to calculate the
eigenvalues accurately since not having the separation of the modes by value of |ms| often produced
a spectrum that was nearly degenerate in the sense that modes with different |ms| had the same
frequency, and in the frequency range of interest. Not only does the 2D formulation prevent that
in the frequency ranges used herein, but it also makes the identification of the ms wavenumber
automatic, since it is part of the PDE as formulated. It allows rapid visual identification of the
mode from looking at the meridional cross section while knowing |ms|. Wokring in 2D comes with
drawbacks though.
First, the formulation in the previous section requires axisymmetry. In the case of 30-cm
setup used in Triana et al. [2014], I was never able to fit the stationary mode frequencies as well
as for the 60-cm setup even after adding ellipticity to the shells after physical examination, by
assuming axisymmetry. Even without ellipticity correction in the formulation, the errors were
never larger than 1% and usually much less. If a mode frequency is at 4kHz though, 1% amounts
to an absolute error of 40 Hz. In the midst of a cluttered spectrum of mode frequencies that can
make identification difficult. Such absolute errors would be even more problematic for a setup
like the 60cm setup where there is no array configuration to determine |ms| by cross correlation
of microphones at different azimuths; they are 90 degrees apart in terms of azimuth in the 60-cm
setup.
3.5 Modal Nomenclature
The altered geometry of the cavity is not so severe as to require abandoning the traditional
wave numbers n, l, and m to characterize the modes, for the spatial wavelengths of the lobes used
in this study at least. Each mode will be identified by a triple (n,l,m). For very high frequencies
there could arise ambiguity with using this notation and the standard interpretation.
Briefly, the overtone number “n will represent the number of spherical nodal curves, a
contour of nodal points for each spherical ray. Examples are depicted in figure 3.1.
46
Figure 3.1: Meridional cross section of pressure modes (1,0,0) and (2,0,0). The color represents
the relative amplitude. The nodal curves in white.
The degree l and order ms are analogous to the wavenumbers appearing in spherical har-
monics. |ms| is the number of the surface nodes that are longitudinal. l− |m| then is the number
of nodal curves emanating from the inner core to the outer core in the 2D meridional cross section
the reduced wave equation is solved on. Throughout, I use the nodal curve counting convention.
This is also the motivation for letting n range from 0, since it signifies the number of nodal curves
being counted; rather than n starting at 1 as many conventions prefer. Once a lobe of the spherical
annulus modes is of similar size as the deviations in the cavity from being a spherical annulus,
ambiguities with naming can arise, but since those modes are not used in this study, the matter
will be addressed here no further. See figures 3.2 and 3.3 for simple examples of the counting
47
principle. When identifying modes visually on the meridional cross section, |ms| is known based
on how the PDE is solved, and so the nodal curves can be counted, and n and l − |ms| can be
inferred. Since |ms| is known, l is then inferred. It should be noted that the naming convention
in no way affects the science, only how the modes are named. Only really does the value of m
matter, which comes from the PDE and is unambiguous.
48
Figure 3.2: Meridional cross section of pressure modes (0,1,0) and (0,2,0). The color represents
the relative amplitude. The nodal curves in white.
49
Figure 3.3: Meridional cross section of pressure mode (1,3,2). The color represents the relative
amplitude. The nodal curves in white. Since m=2 is known from how the PDE is solved, the
nodal curves can be counted. See that n=1, and l − |ms| = 1, and so l = 3. No inconsistencies
or ambiguities result from adopting this naming convection in the altered geometry within the
frequency range and small wavenumbers used in this study.
50
Chapter 4: Numerical Formulation
4.1 Finite Element Formulation
The flexibility of the Finite Element Method (FEM) to accommodate irregular geometries,
and allow for future refinements to the geometric model of the cavity, make it ideal for simulating
acoustic waves in experimental setups. The more precise the modeling of the geometry, the more
accurate the center frequencies will be predicted, and the easier it will be to identify modes. This
is crucial when the sensor array is limited. More importantly, the presence of altered geometries
within or near the tangent cylinder—where arguably the more interesting variations in flow might
be observed, the more accurately inversions will infer the flow in those regions by using FEM. The
variational formulation does not require uniform soundspeed. However, soundpseed variation is not
utilized in this study since not only are the spectra observed consistent with a uniform soundspeed,
but also significant soundspeed variations are not expected at the velocities and thermal behaviors
present in the experiments.
I use a piecewise triangular mesh sufficiently fine to allow an accurate boundary and conver-
gence in the solutions, consisting of over one hundred thousand triangles. The particular meshing
algorithm I employ is based on the methodology set forth in Persson and Strang [2004], with minor
modifications. I set a minimal bound on the regularity of the mesh and mandate a triangle quality
greater than 0.75 for every triangle, with the triangle quality metric given by
µ(T ) =
4AT
√
3
h21 + h
2
2 + h
2
3
(4.1)
where h1 h2 and h3 are the side lengths of the triangle T, and AT is the triangle area. More than
51
80% of the triangles have a µ value greater than 0.99 though; the metric equals one when the edge
lengths are all equal. I first create a coarser mesh, then refine it twice by way of breaking every
triangle into 4 subtriangles via placing new vertices between edge midpoints of the original triangle
and connecting them. The large initial mesh size prior to refinement is chosen so as to ensure when
the nodes along the boundary of the domain are moved to the boundary specified by the descriptive
geometry after subdivision, the regularity of the triangles only changes negligibly. That is, the
refinement process also improves the boundary accuracy while not eroding the regularity of the
surrounding triangles.
I use isoparametric elements, where the basis functions and test functions are piecewise
linear splines; “tent functions”. Each basis function is unity at a single node, and linearly decays
to zero at the neighboring nodes, and is zero outside that neighborhood. Given such a discrete
basis {ψj : j = 1, 2, ..., N} for the solution space, any solution ψ can be approximated in the form
ψ =
N∑
j=1
ajψj (4.2)
where the coefficients aj : j = 1, 2, .., N are to be determined that best approximate the solution
with respect to the basis {ψj : j = 1, 2, ..., N}.
Deriving from Eq. 3.6, the stiffness matrix elements are given by
Kjk = −
∫
A
r∇ψj · ∇ψk dA−
∫
A
m2s
r
ψjψk dA (4.3)
Similarly, the mass matrix is given by
Mjk =
∫
A
rψjψk dA (4.4)
The integrals for the matrix elements for the stiffness matrix are assembled in two parts. The
first integral is the same regardless of ms. For that contribution, a centroid rule for estimating the
integrand is used, since the gradients on the piecewise linear elements make the integrands piecewise
52
constant inside each element using an exact Gauss quadrature rule. For the other integral and
the mass matrix, the following scheme has been used: the coefficient functions are interpolated
to their centroids on each triangle and moved out of the integrals, and the remaining quadratic
integrands ψjψk are then integrated by an exact Gauss quadrature rule. The presence of the axle
in the model alleviates the singularity in the
m2s
r factor.
4.2 Solving for the Stationary Modes and Stationary Frequencies
The modes and frequencies of the stationary medium are solved for using sparse Arnoldi
iteration [Saad, 1980], carried out with ARPACK [Lehoucq et al., 1998] via interface with Matlab.
The assembled finite element matrices are symmetric, so the algorithm used is Implicitly Restarted
Lanczos Method (IRLM).
Since a different system is solved for each value of ms, there is greater spacing between
the adjacent modes throughout the low frequency band of the spectra used, even greater since
the eigenvalues solved for are actually −ω2; alleviating any numeric issues associated with near
degeneracies. The regularity of the mesh, and the low frequency range being used, means that
convergence of the eigenvectors and eigenvalues with respect to triangle density also implies con-
vergence of the solutions; there are no significant effects of round-off for the precision needed; the
matrix elements having very similar order of magnitude. The larger source of error, though still
quite small, comes from idealizations in the geometric modeling.
The actual experimental setups such as the 60-cm setup contain small deviations from the
descriptive geometry used to construct the mesh. The joint of the shaft and outer sphere is not
a sharp corner, but a fillet. Furthermore, the shaft has missing pieces, which were estimated in
the descriptive geometry from mechanical drawings and photographs, but the actual components
contain all additional pieces not modeled. The microphones and speakers also take up space, though
not sufficiently large to merit incorporating them into the model. These unmodelled aspects will
introduce quite small errors. All else held equal, they are quite insignificant to the inversion
methods, but even a 0.25% error at a mode frequency of 2kHz amounts to a 5Hz shift in terms of
53
absolute frequency, which can make it difficult to identify a mode based on its frequency and the
limited knowledge of ms, especially when the sensor locations are suboptimal.
Near degeneracy or small spectral spacings between modes for the same value of ms may
cause higher order effects to become significant. In the data this work has applied its methods
to though so far, clusters of such modes are often difficult to measure the splittings of since the
available array techniques on the 60-cm setup do not adequately distinguish them. Thus, all the
modes used of the same value of m are distantly spaced in the frequency spectrum, and the issue
is circumvented by default.
The (n, l,ms) identification was the performed by viewing a meridional cross section pro-
duced by the FEM, and then based on counting nodal curves, as explained in the previous section.
All modes in this work were checked visually by examining the pressure pattern on the meridional
cross section to ensure there were no mistakes in identifying n or l of the associated multiplet from
the FEM solutions. Several automation methods were developed, such as using operators on the
FEM solution and inferring the (n, l,ms) from the eigenvalues, as well as cross correlation with
the spherical annulus solution. The methods produced different results for high frequencies, since
cross correlation with the spherical annulus modes is inconsistent with the nodal curve naming
convention in some instances for high l values for example. Cross correlating in iteratively with
modes as the shaft is grown from the spherical annulus was infallible at coinciding with the nodal
curve interpretation of (n, l,ms) identification; visual identification was simple enough to not risk
incurring an error associated with not checking though.
Table 4.1 shows the results of the FEM for the modes up to 2.5 kHz. It uses a soundspeed of
c = 344 ms . Stationary frequency accuracy was high enough using only a single mode for calibration
of soundspeed that least squares methods were unnecessary. The actual stationary frequencies are
not used in inversions for the regimes applied to in this study, since the stationary frequencies are
only needed when the speed of rotation is much higher, and the Helmholtz equation’s linearity in
−ω2 must be used rather than linearizing the perturbations δω2 to perturbations in δω; discussed
in chapter 5. Even in the case that pertubrations in the square of the angular frequency are prudent
54
to use, the error incurred by approximating the center frequency is quite small for plausible rotation
rates for these types of studies.
55
n l ms f (Hz) n l m f (Hz) n l ms f (Hz)
0 1 1 346 0 7 5 1605 1 6 2 2167
0 1 0 365 0 7 0 1695 2 3 0 2172
0 2 1 570 1 4 1 1702 0 10 3 2173
0 2 2 593 2 0 0 1728 1 6 1 2179
0 2 0 612 1 4 2 1734 0 10 4 2181
0 3 1 774 1 4 3 1749 0 10 10 2181
0 3 2 807 1 4 4 1750 0 10 9 2181
0 3 3 809 0 8 1 1760 0 10 8 2181
0 3 0 837 0 8 2 1768 0 10 7 2181
1 0 0 932 2 1 0 1790 0 10 6 2181
0 4 1 970 2 1 1 1791 0 10 5 2181
0 4 2 1009 0 8 3 1795 1 6 3 2198
0 4 3 1014 0 8 4 1798 1 6 4 2203
0 4 4 1014 0 8 8 1798 0 6 6 2203
1 1 1 1045 0 8 7 1798 1 6 5 2203
0 4 0 1056 0 8 6 1798 2 4 1 2281
1 1 0 1071 0 8 5 1798 0 10 0 2308
0 5 1 1165 1 4 0 1847 2 4 2 2313
0 5 2 1204 2 2 1 1898 0 11 2 2321
0 5 3 1213 0 8 0 1908 2 4 3 2329
0 5 4 1214 2 2 2 1922 2 4 4 2330
0 5 5 1214 1 5 1 1938 1 6 0 2348
1 2 1 1233 2 2 0 1940 0 11 3 2359
1 2 2 1263 0 9 2 1949 0 11 1 2361
0 5 0 1270 1 5 2 1957 0 11 4 2371
1 2 0 1313 0 9 1 1964 0 11 11 2372
0 6 1 1361 1 5 3 1979 0 11 10 2372
0 6 2 1395 1 5 4 1981 0 11 9 2372
0 6 3 1409 1 5 5 1981 0 11 8 2372
0 6 4 1410 0 9 3 1985 0 11 7 2372
0 6 6 1410 0 9 4 1990 0 11 5 2372
0 6 5 1410 0 9 9 1990 0 11 6 2372
1 3 1 1461 0 9 8 1990 1 7 2 2375
0 6 0 1483 0 9 7 1990 1 7 1 2404
1 3 2 1500 0 9 6 1990 1 7 3 2411
1 3 3 1507 0 9 5 1990 1 7 4 2419
0 7 1 1560 2 3 1 2060 1 7 7 2419
0 7 2 1583 1 5 0 2083 1 7 6 2419
1 3 0 1587 2 3 2 2099 1 7 5 2419
0 7 3 1603 2 3 3 2105 2 4 0 2452
0 7 4 1605 0 9 0 2124 0 12 2 2508
0 7 7 1605 0 10 2 2135 0 11 0 2520
0 7 6 1605 0 10 1 2155 0 12 3 2543
Table 4.1: FEM calculated frequencies at c=344 m/s for the 60-cm setup. As can be seen, the
multiplets are split due to the shaft’s presence. Values have been rounded to the nearest natural
number. The data the methods are applied to is as high as 2.1kHz; the reason for the range
presented.
56
4.3 Discussion of FEM Predictions for Stationary Modes
The shaft lifts the degeneracy of the multiplets, as can be seen from the predictions in table
4.2 . Better inter multiplet agreement occurs for very high frequency modes though since they are
less sensitive to the shaft often. However, modes are more affected by the presence of the shaft
are the ones that carry the most information on the tangent cylinder’s interior zonal flow, so they
are of high utility for reconstructing zonal flow inside the tangent cylinder. Helioseismology for
instance, has evolved to use G-modes—modes whose restoring force is gravity—for imaging deeper
structure. Pressure modes are quite useful for inferring the azimuthal angular velocity in the sun’s
convective zone and to some extent the upper radiative zone, yet have weaker resolution for the
deeper interior of the sun and near the poles.
4.4 Solution for a Spherical Annulus
In the absence of a shaft, axle, or other components that present themselves in experimental
work, the acoustic Helmholtz equation in a rigid spherical annulus is possible to solve via separation
of variables in spherical coordinates. Solving this form of the problem provides a comparison of
the frequencies without the shaft or axle to those of the actual setup. Additionally, comparison to
my FEM algorithm with altered parameters provides another quality assurance for the numerics.
4.4.1 Modes of the Spherical Annulus
To avoid confusion, in this section the variable r will stand for the spherical radial variable,
and U will be used for the radial wavefunction variable. Consider a solution of the form
P(r, θ, φ) = U(r) Θ(θ) Φ(φ) (4.5)
where P is pressure, and has been notated in this way so as to avoid confusion with the Legendre
polynomial symbol. The functions Θ and Φ are invariant with respect to the size of the inner
57
sphere, and they are given by
Θ(θ) = Pmsl (cos θ)
Φ(φ) = eimsφ
(4.6)
where P are the associated Legendre polynomials, and the product of these azimuthal and polar
functions are spherical harmonics Y ml (θ, φ) = P
m
l (cos θ)e
imφ. The radial wave functions cannot
be found purely analytically though. They can be expressed in terms of spherical Bessel functions,
but the boundary conditions require the wavenumbers to be computed numerically, in all but a
few specific ratio of radii of the outer boundary to the inner. The solutions are of the form
Unl(r) = J
sp
nl (r) + bN
sp
nl (r) (4.7)
where Jsp and Nsp are spherical Bessel functions of the first and second kind respectively. These
solutions are analytic if k and b are known, but b and k must be calculated numerically for most
aspect ratios RoRi .
The radial ordinary differential equation in spherical radial variable r is given by
d2U
dr2
+
2
r
dU
dr
− l(l + 1)U
r2
= −k2U (4.8)
where l is the degree and k is the wavenumber, the ratio of the acoustic angular frequency to the
soundspeed.
4.4.2 Computing the Solution to the Radial Wave Equation
Since the inner core is roughly 1/3 the diameter of the outer, no singularities or near
singularities are present in spherical coordinates. For each value of l I solve for the wavenumber
k using the boundary condition dUdr = 0 at r = Ri and r = Ro. The frequency of the mode then,
is the product of the wavenumber and the soundspeed. Each frequency is implicitly a function
of l but has no ms dependence in this geometry. Several methods are employed to check the
accuracy, so as to have a very accurate ground truth solution. Both a finite element formulation
58
and finite difference formulation in 1D are used, being highly accurate since the equations are
defined over a 1D interval only. Both the FEM and FDM reduce to eigenvalue problems and are
solved using sparse Arnoldi iteration. Visual inspection of the wavefunctions produced, as well
as performing a grid search over a neighborhood of k and the inferred b value was performed to
ensure high accuracy, by way of the analytic solutions; a simple means of ensuring the very high
number of nodes used did not produce any numerical round-off concerns. Both the 1D FEM and
1D FDM were sufficiently accurate that the slight changes in k observed were likely due to numerics
associated with the fitting of k and b simultaneously. Error in frequencies on the modes presented
here, are estimated to be less than 0.01 Hz in error, much more accurate for most, especially low
l values.
The finite element formulation uses the following equivalent form
d
dr
(r2
dU
dr
)− l(l + 1)U = −k2r2U (4.9)
then proceeds with using piecewise linear splines and exact Gauss quadrature to formulate the
stiffness matrices Kl : l = 0, 1, ... and the mass matrix M. FEM incorporates Neumann boundaries
must more naturally, as the variational formulations works by performing integration by parts on
the PDE and the normal derivative at the boundary being zero nullifies the contribution of the
boundary term. With FDM ghost nodes must be used at least implicitly in how the boundary is
handled, which means that to use a few nodes implies handling the length of the last node spacing
in a way so as to not incur error. FDM presume the solution goes completely flat at the boundary
over the last interval. The two methods produce the same answers at sufficient grid resolution
though, and so FEM was adopted.
Finite difference operators were constructed using centered differences and ghost nodes at
the boundaries to enforce the boundary conditions. This alleviates fitting the mode to the geometry
with both k and b simultaneously. At these low frequencies the accuracy and speed of solving the
ODE via iterative eigenvalues problems is exceptional, in either method.
59
4.5 Comparing Theory and FEM for Spherical Annulus Solution
In table 4.2 I compare the obtained frequencies with those from my finite element code
with no axle or shaft; the shaft radius was set to 10−6 (m) in order to prevent singularity when
evaluating −m
2
r ψ. Presented in table 4.2 are the 1D FEM coded values, the 2D FEM values, the
relative error of the 2D from the 1D, along with the maximum absolute deviation in Hz that the
ms 6= 0 modes deviated from the ms = 0 mode in the 2D finite element solutions. Each multiplet
is labeled by the nSl convention.
60
Multiplet f1D (Hz) f2D (Hz)
|f2D−f1D|
f1D
‖f2D − fm‖∞ (Hz)
0S1 356.93 356.93 1 · 10−7 0.02
0S2 593.70 593.71 5 · 10−6 0.05
0S3 808.86 808.87 1 · 10−5 0.07
1S0 931.01 931.02 1 · 10−5 NA
0S4 1013.76 1013.77 1 · 10−5 0.12
1S1 1055.18 1055.20 2 · 10−5 0.05
0S5 1213.50 1213.51 1 · 10−5 0.15
1S2 1264.20 1264.23 3 · 10−5 0.12
0S6 1410.21 1410.23 2 · 10−5 0.21
1S3 1507.12 1507.16 3 · 10−5 0.23
0S7 1604.90 1604.94 2 · 10−5 0.20
2S0 1733.86 1733.97 6 · 10−5 NA
1S4 1749.74 1749.80 3 · 10−5 0.32
2S1 1797.58 1797.70 7 · 10−5 0.79
0S8 1798.10 1798.15 3 · 10−5 0.06
2S2 1922.90 1923.02 6 · 10−5 0.10
1S5 1981.15 1981.12 3 · 10−5 0.09
0S9 1990.13 1990.21 4 · 10−5 0.08
2S3 2104.57 2104.70 6 · 10−5 0.08
0S10 2181.21 2181.32 5 · 10−5 0.10
1S6 2203.02 2203.13 5 · 10−5 0.11
2S4 2398.23 2329.96 6 · 10−5 0.10
0S11 2371.50 2371.66 6 · 10−5 0.12
1S7 2419.03 2419.17 5 · 10−5 0.13
Table 4.2: Normal Mode Frequencies computed by FEM in column 3 agrees quite well with the
theoretical semi-numerical result in column 2, as seen in the relative errors in column 4. The
differences between frequencies of the elements in the multiplets are quite small as well; column 5.
The NA values are reflecting that only ms = 0 is possible when l = 0. Note: the estimated relative
errors are based on the unrounded values.
61
The errors shown in 4.2 are quite small. The fit with theory is quite good. The near
divergent behavior of 1r in the absence of the large shaft is likely to blame for the discrepancy as
much as discretization error and other concerns. A finer mesh was not used to compensate for that
behavior. The mesh was kept as similar to the one used in the computations with the shaft as
possible, in order to make as accurate of a comparison as possible. The modeling uncertainty likely
far outweighs error due to discretization. The mesh is denser than needs be for accurate frequency
identification, but the same mesh is used for flow parametrization and Kernel calculation, which
requires smoother parametrization and numerical differentiation when computing kernels and can
introduce a small amount of error as well. If a mesh was used at the minimum resolution level
needed to accurately predict the needed frequencies for the 60-cm setup, it would have been far
coarser and more disagreement would have occurred within the multiplets. The mesh was also used
for the flow parametrization though, and for numerical integration associated with computing the
sensitivity kernels and forward operator discussed in the next chapter. A much finer mesh was
used than the acoustics mandated, and so there was very high accuracy. The resolution level would
perhaps be inadvisable in terms of efficiency for an analogous problem in 3D that could not be
reduced via axisymmetry or to a scalar potential such as pressure perturbation.
62
Chapter 5: Linear Theory of Acoustic Splitting for Rotating Flows with Shear
5.1 Wave Equation inside a Vortex
The homogeneous medium wave equation in the presence of a material velocity field within
the medium is given by
u¨ = ∆u + 2(v0 · ∇)u˙, (5.1)
where again, units have been chosen so that the soundspeed is set to unity. This leads to the
perturbed Helmholtz equation:
− ω2u = ∆u + 2iω(v0 · ∇)u, (5.2)
and the resulting first order correction in ω:
δω =
−i ∫
V
u†(v0 · ∇u)dV∫
V
|u|2dV . (5.3)
These expressions do not have a density distribution or soundspeed term weighting the
integrands, which would be required if the variations induced by the flow were significant relative
to the resolution of the inversion being sought.
If v0 is azimuthal, i.e. purely directed in the cylindrical coordinate φˆ, then
v0(r, z) = rΩ(r, z) φˆ (5.4)
63
and thus for a mode with spin ms in an azimuthal flow
(v0 · ∇)u = imsΩ u + Ω z× u (5.5)
The first term is associated with the affect of advection, and the second term is associated
with the affect of Coriolis forces. In the case of uniform rotation, for instance, advection makes
the fluid medium rotate such as to increase the frequency of the prograde mode spinning about
the same axis, and decrease the frequency of the mode spinning retrograde to the fluid motion.
In the this work I will use the convention of the −|ms| mode spinning prograde and the +|ms|
spinning retrograde. This is because the frequencies of the modes are considered to be positive,
the distinction between ±ω not being made. And so, if the azimuthal envelope eimsφ is modulated
by eiωt, then the combined effect is a spinning envolope given by ei(ωt+msφ). The pattern of the
mode can be seen to spin by tracking the phase value α by ωt + msφ = α, i.e., φ(t) = α − ωt.
The pattern is spinning backwards at an angular rate of ω radians per second. Inserting Eq. 5.5
into Eq. 5.2 also produces a sum of two self adjoint operators like the Laplacian; proving that
they do not destroy the discrete eigenstate structure necessitating a singular perturbation theory
interpretation is beyond the scope of this work, though insofar as one assumes the problem can
be discretized, i.e. written as a matrix formulation, it is a triviality and follows from the spectral
theorem.
Next, I modify these canonical expressions to be expressed in the derived pressure mode
cross sections that I have calculated via FEM.
Consider the auxillary function η(r, z) given in terms of the auxillary wave function ψ
representing the acoustic pressure pertubration at φ = 0.
η(r, z) := r(
∂ψ
∂r
)2 +
m2
r
ψ2 + r(
∂ψ
∂z
)2 (5.6)
η simplifies the expressions for the corrections. The calculation for the first order correction to
angular velocity reduces to
64
δω = −m
∫
A
Ω(r, z)η(r, z)dA∫
A
η(r, z)dA
+m
∫
A
Ω(r, z)∂|ψ|
2
∂r dA∫
A
η(r, z)dA
(5.7)
where A is the cross section at φ = 0. I derived this from applying first order perturbation
formula for the eigenvalue −ω2 for a normalized material displacement mode solving the Helmholtz
equation, and then using the approximation
− ω2u = − 1
ρo
∇p (5.8)
where the variable p is the acoustic pressure perturbation such that P = Po + p. The expression
is then fully linearized the by the Taylor approximation of
ω =
√
ω2o + δω
2 (5.9)
which implies that
δω =
√
ω2o + δω
2 − ωo (5.10)
which yields
δω = ωo(
√
1 +
δω2
ω2o
− 1) (5.11)
expanding
δω ≈ ωo(1 + 1
2
δω2
ω2o
− 1) (5.12)
hence
δω ≈ δω
2
2ωo
(5.13)
What is measured is δω, or δf . δω2 = (ω + δ)2 − ω2. Using the approximation of Eq.
5.13 introduces an error of approximately (δω
2)2
8ω3o
. For the data used in this study with the highest
rotation rates and highest splittings, the error is less than 1%, on average 0.25%, and for the
modest rotation rates used in inversions in chapter 8 even lower. No significant improvement can
be made by trying to incorporate the finer correction by adjusting based on center frequency,
65
especially with the uncertainty associated with measuring the center frequency or the temperature
acoustically.
As can be seen from the frequency correction expressions, there is no soundspeed variable
or material parameter such as density or bulk modulus. This implies that to first order the
temperature drifts of the system that change the soundspeed inside the cavity, will only shift the
center frequencies and not alter the magnitude of the splitting about the center frequency. In
the first order regime, measuring the splitting between +|ms| and −|ms| will not require drift
correction. For the regimes of flow being studied, and the magnitudes of splitting observed, the
small fraction of a percent the perturbation formulation in terms of −ω2 would offer would be offset
by the error incurred due to ωo center frequencies requiring drift correction, and possibly higher
order corrections. It can also be seen that ms = 0 modes do not split due to azimuthal flow, to
first order. The spin of the mode, is purely a spin of the pressure relative amplitudes, not a spin of
the material medium, and so the axisymmetric ms = 0 modes do not split by virtue of remaining
invariant along any azimuth. This assumes a relatively uniform temperature and pressure inside
the cavity.
The integrals in Eq. 5.7 can be viewed as linear operators on the rotational flow, with
kernels defined by
κ(r, z) = −ms η(r, z)∫
A
η(r, z)dA
+ms
∂|ψ|2
∂r∫
A
η(r, z)dA
(5.14)
The shift in angular frequency of a mode ψ with spin ms ins the presence of an azimuthal flow
having rotational profile Ω(r, z) is given by
δω =
∫
A
κ(r, z)Ω(r, z)dA (5.15)
The kernels κnlms are computed numerically by evaluating the expressions at triangle centers,
components of the gradient derived from the gradient of the tangent plane of the piecewise linear
solution over each triangle in the mesh. Testing of theory was done in several stages using chirp
based data. Uniform rotation, where the fluid is spun up to rotating like a rigid body was used.
66
Analysis for differential rotation and the presence of shear is carried out for the actual data used
in inversions by way of the ms = 0 modes in chapter 8.
5.2 Uniform Rotation
Uniform rotation of the fluid inside the cavity is often called ”solid body rotation”, since
the material moves analogous to a rigid object under rotation. If the container is spun up slowly in
the appropriate way, both inner and outer boundaries rotating in unison, approximately uniform
fluid rotation results. Using data collected in Adams [2016] I test my software against the results.
The data was collected in the rotating frame, meaning that it is nominally at rest in that frame.
Relative to the rotating frame, the equations simplify to
− ω2u = ∆u− 2iωΩzˆ× u (5.16)
Only the Coriolis effect contributes since the fluid is not advecting relative to the reference frame.
Table 5.1 shows the results of theory versus observed shifts from the center frequencies.
n l m f (Hz) FEM f (Hz) observed % error
0 1 1 346.0 346.0 Calibrated
0 2 1 570.2 569.5 0.12
0 2 2 593.0 594.5 0.25
0 3 2 806.7 807.5 0.10
0 3 3 808.8 810.5 0.21
0 4 2 1008.9 1010.0 0.11
0 4 3 1013.6 1015.0 0.14
0 4 4 1013.7 1015.5 0.18
1 1 1 1045.4 1048.0 0.25
Table 5.1: A comparison of my FEM predictions with observed frequencies in the stationary
medium by targeted chirp, and errors are a small fraction of 1%. The error is likely even slightly
lower, since the soundspeed was calibrated using the lowest frequency, rather than the highest, or
some least squares fit to all. The measured acoustic data in the column for observed frequency
were performed by Matthew Adams.
67
∆1 FEM
∆5
5 observed error %
∆10
10 observed error %
∆15
15 error%
0.660 0.660 0.072 0.664 0.602 - -
0.214 0.213 0.346 0.212 0.580 0.205 4.089
0.728 0.730 0.249 0.725 0.438 0.707 2.956
0.459 0.455 0.792 0.460 0.298 0.462 0.662
0.704 0.700 0.580 0.700 0.580 0.708 0.604
0.315 0.310 1.438 0.317 0.946 0.315 0.151
0.502 0.460 8.366 0.503 0.100 0.500 0.398
0.670 0.650 2.995 0.673 0.363 0.687 2.477
0.135 0.155 15.164 0.150 11.449 0.132 2.173
Table 5.2: Comparison of the same mode set from the table 5.1 , at various rotation rates where
the inner and outer boundaries are spinning together, spun up to achieve uniform rotation of the
fluid. The ∆ values are the absolute shift of the modes from the center frequency in Hz. The
denominators on the observation ∆ values are the rate of medium’s rotation in Hz. The leftmost
column represents the expected shift from center frequency due to 1 Hz uniform fluid rotation.The
percent errors do not take any account of the observational uncertainties or modeling uncertainties.
The missing data in the first row is because measurement could not be reliably made due to that
portion of the spectrum at the high rotation rate of the experiment being very noisy. Under the
presumption that there should be no splitting at 0 Hz then, the linearity implies that the ratio of
the splitting to rotation rate should remain invariant during experiments and equal to the finite
element prediction. Note: This table does not give signs to the shifts, since signs were not measured
for the data, only the magnitude of the splitting was measured. Acoustic data measurements made
by Matthew Adams.
68
Table ?? shows a comparison of the linear theory of splitting with observed splittings from
the center frequency. The majority of measurements are within 1% error. The last row of data
corresponding to mode (1,1,1) is significantly higher in terms of percent error, though in terms
of absolute deviation is still quite small. Frequency resolution in the acquired data is comparable
to the amount of splitting being measured in those cases, and are not evidence of any theoretical
issue. Mode (1,1,1) at Ω = 10 Hz for instance has a shift of 1.50 Hz observed compared to FEM
prediction of 1.35 Hz. This translates to an error 11.5%. The small absolute error though, can be
explained by the observational uncertainty. Reading the splittings off the graph in such conditions
carries an uncertainty of roughly 0.4 Hz, translating to 0.2 Hz shift from center frequency. Errors
are also relatively high for the rotation rate of 15 Hz. Spinning the experiment at 15 Hz as opposed
to 10 Hz or lower produces significantly more noise, contaminating the spectrum.
5.3 Soundspeed Uniformity
The previous section assumed the soundspeed was uniform. Thermal stratification due
to the mechanics of the apparatus or in conjunction with velocity induced compression was not
incorporated in this study. The velocities were sufficiently low, |rΩ|2 <<< c2air, so the kernels did
not require weighting functions related to the soundspeed distribution.
Assuming incompressibility of liquid sodium, the deviations from unity weighting are at
most a small fraction of one percent for the experimental conditions analyzed here, much less for
the bulk of the cavity. The sodium experiment is 1.5 meters in radius and has a max outer rotation
rate of 4Hz, and the soundspeed of sodium is roughly 2.5 km/s.
Higher compressibility of air compared to sodium somewhat complicates matters, but, to
first order, the effects on soundspeed should be somewhat smaller for air, since the slight pressure
increase would induce slight density increase thereby preventing the soundspeed increase partially.
Corrections to the soundspeed distribution can be made using only the ms = 0 modes and a
traditional approach as used in seismology to infer the soundspeed profiles at long wavelengths:
the soundsdspeed squared can can be separated into a background model and a perturbation
69
quite accurately for small variations, and the Laplacian’s linearity can be employed to produce a
perturbation expression for use in a linear inversion. Given the uncertainties present in the data, it
is not possible to do so without conflating the effect with the noise and observational and modelling
uncertainties. As such, the matter has not been pursued herein.
Temperature variations were also considered. The uniform rotation of the air presumably
produces less mixing than in the data later used for inversions which were collected during differ-
ential rotation. Supporting that there were no significant effects due to that were occurring, the
ms = 0 modes were monitored. For the data used in chapter 8, 10 m=0 modes were monitored
across different rotation rates of the inner sphere, between -4 Hz and -40 Hz, and the m=0 mode
changes in frequency could be explained by a uniform temperature increase in the cavity. One
caveat is that slightly more elevated temperatures in regions with low acoustic sensitivity to flow
would not affect acoustics appreciably. Discussed further in chapter 9.
70
Chapter 6: Relating Fluid Rotation to Acoustics
6.1 Flow Characteristics
In this work mean fluid flow has been approximated as being purely azimuthal in the lab-
oratory frame of reference insofar as acoustic splittings are concerned, and flows where that ap-
proximation is valid are the interest herein. That is, centrifugal forces and poloidal flows are not
strong enough to sufficiently distort the mode shapes thereby influencing the splitting between
modes with equal and opposite azimuthal order ms. Poloidal flows can shift the center frequencies
though, and could be one of the explanations for the slight movements of spectral peaks seen in
the data depicted and discussed in chapter 8.
Acoustical modes sample the medium continuously with time, and so time varying flows
distort the idealized steady state oscillation of a mode. If the fluctuations are sufficiently small,
spectra acquired from sufficient time integration produce spectra with minimal defocussing of mode
peaks. The technique images the mean flow, where the duration over which the mean is taken
is the time interval over which the acoustic spectra are estimated. This type of mean flow is a
longer term average than simply averaging through the turbulence. The flow need not be steady
in the traditional sense; even if there is transient oscillation between nearby states via hysteresis
or inertial waves the acoustics record the mean over many seconds to as much as thirty seconds or
more.
6.2 Parametrization Style for the Flow
In this section, I describe the approach to parametrize flow used in most methods in this
work. Specific parametrization are addressed in chapters 8 and 9.
71
The flows are represented by vectors of coefficients with respect to a set of basis flows; a
flow being the superposition of the basis flows weighted by the coefficients. The linearity of the
frequency shifts due to flow then implies that the splitting can be calculated for the flow based on
those coefficients and the amount of splitting attributable to each basis flow. These coefficients
will be referred to as the model parameters, and the space in which they reside as the model space
M. If m ∈M, then it is associated to the flow Ω by
Ω =
dimM∑
j=1
mjΩj (6.1)
where {Ωj : j = 1, .., dimM} are the basis elements used for the flow parametrization. I adopt the
practical approach of forcing the dimension to be finite. Even for high resolution inversions the
expansion will be truncated at some index relative to the resolution needed or is reliable. Often
there is a limit imposed by the discretization as well, though well beyond what is resolvable in this
context.
Other styles of parametrization such as using model space parameters that are exponents
of basis flows or spatial variation wavenumbers were not used, as such styles of parametrization do
not as easily lend themselves to linear inversion. Since few modes are used in chapter 8 to perform
inversions, those alternate styles of parametrization are quite feasible in those settings, but this
work is intended for use beyond those demonstrative settings.
Parametrizing the flow well then means finding functions that can optimally capture large-
scale features of the mean flow when the flow is averaged over such time scales. Basis functions are
often remarked as being ”complete”, which means simply that for any member of a particular space
of functions the basis is complete with respect to, that function can be represented as a convergent
series expansion of those basis functions. When using an implicit regularization rather than a
traditional Tikhonov style damping strategy though, the effects of truncation of the series can be
more severe. For instance, even though products of zonal harmonics and spherical Bessel functions
can represent cylindrically symmetric functions well, i.e. they form a complete basis, it is much
harder to do so accurately when truncating to a dozen or so terms. For that reason, as discussed
72
here and in chapter 8, such basis functions were supplemented with vertically invariant cubic splines
in cylindrical radius r. Producing largely vertically invariant flows using many spherical based basis
functions for the meridional cross section requires far more basis functions than the flow requires
or the data can support. Using more basis functions, furthermore, results in a system of equations
far too underdetermined to be solved without imposing some form of explicit regularization, a
regularization that in this framework would be difficult to justify if it were constructed in a way
that allowed small length scale features not actually resolvable by the data.
To overcome this difficulty, a hybrid basis is primarily used in chapter 8, whereby a collec-
tion of column like splines that are a function of cylindrical radius r alone are used in addition to
the products of spherical Bessel functions with zonal harmonics; specifically, the Neumann modes
of the spherical annulus. The strategy is to allow the vertically invariant component to be largely
represented by the splines, while the superimposed deviations from cylindrical symmetry are repre-
sented by the other zonal harmonic related components. In order to assess that the reconstructed
flows were not simply mimicking the basis functions unrealistically, other basis functions such as
polynomials in r and z as well as 2D Cubic splines were used to gauge the effects of the basis in
biasing the results obtained in chapter 8 and discussed therein.
6.3 The Forward Operator
The forward operator C is a map from the set of admissible flows to the space of observations
of resulting frequency shifts. Within the linear regime, using a finite set of basis flows, C is a matrix
satisfying
Cij =
∂fi
∂mj
(6.2)
The range of C is the set of vectors of frequency shifts for the set of acoustic modes being imaged.
Thus Cij is the amount of frequency shift of the i
th mode due to basis flow j having coefficient 1.
The frequency shift is half of the splitting between the ±|ms| pair.
Computation of the forward operator allows for testing hypotheses about the flow without
73
explicit inversion. The forward operator can be calculated for even a very large basis for some inner
product space, and any flow can be projected the basis functions for that space, and the coefficients
can be formed into a model space vector that can have the related forward operator applied to it
to see what the acoustic effects are of that flow. In that regard, the forward operator can reduce
the computation time of applying the kernels to the flows directly, and utilize the coefficients of
the flow to rapidly compute the acoustic effects by matrix multiplication.
Another way of thinking of the forward operator, is that it is the matrix product KTB, the
set of inner products of the kernels with the flow basis functions; here the volume weighting by dV
implicitly absorbed into K where each column of K is a kernel for a specific mode appropriately
weighted at each element to allow the matrix product to carry out the volume integral.
6.4 Boundary Conditions of the Fluid Flow
No boundary conditions are imposed on the flow for inversions. Even if a no slip boundary
condition is approximately true, i.e. one where the flow is following the boundary at each location
of the boundary, the width of the boundary layer that decays to the boundary velocity at the wall
is suitably thin so as to prevent accurate reconstruction without greater amounts of data. The
viscous layer itself is so thin that possibly no amount of data is capable of accurately inferring that
acoustically. That discourages the use of a set of basis functions that have imposed boundary values
nor enforce a constraint that mandates their superpositions have prescribed boundary values; in
the absence of a prior model or data supporting those boundary values at least. Fortunately
though, this lack of resolution in that regard, is because such a thin boundary layer has little affect
on acoustics. The only boundary that would possibly even need included in error analysis is the
outer shell’s; yet for SCF in the hydrodynamic case the vast majority of the fluid at the outer shell
is following the outer shell rotation rate, in terms of azimuthal velocity component. So no such
error analysis is carried out in that regard. Perhaps if the fluid was not brought to a relatively
steady state before acoustics were acquired that would not be the case.
74
Chapter 7: Techniques of Inversion
7.1 Constrained Optimization
Relying only on the data and forward modeling of the splitting due to zonal flows as described
in chapter 5, leads to non-unique results. With limited data, the non-uniqueness is accentuated.
Using Normal Equations, i.e. minimizing weighted least squares, often mandates fewer basis
functions than data points. This is because the correlation and complementarity and coverage is
dissimilar from observing the flow at 15 directly sampled points spaced out in the domain; many
modes’ sensitivity pattern makes them more analogous to very closely sampled points with large
uncertainties on their functional value at those points. Introducing additional modes does not
necessarily increase the resolution of the inversion, but simply reduces the uncertainty in the
measurements the other modes that have similar sensitivity gauge.
Additional basis functions can only be added when using the Normal Equations if some type
of damping, truncated singular value decomposition, or other regularization method is employed.
Choosing the regularization parameters to do that can be difficult and highly subjective. Intead,
more information can be extracted from the data if it is supplemented with physical constraints
that lift the non-uniqueness of the inverse problem produced by a finer basis.
Consider the traditional weighted least squares problem
Φ(m) = ‖Wd(Cm− d)‖2
mest = m
∗
Φ(m∗) < Φ(m) ∀m 6= m∗
(7.1)
75
where Wd is some form of weighting of the data based on the confidence, or inverse data covariance
matrix. I solve the problem using techniques from nonlinear optimization theory. The model space
is decomposed into regions, some considered physically realistic and others not; here the model
space is simply the values that the coefficients of the basis functions reside in. Regions that are
considered physically realistic from the feasible set, and those considered physically unrealistic have
Φ(m) =∞. To minimize data misfit under this type of contraint multiple methods are employed.
The most robust method employed is the Barrier Function Method (BFM). That is, penalty terms
are added to the objective function that have negligible effect inside the feasible set develop large
and steep gradients as the model vectors m approach the boundary of the feasible set. For a set
of barrier functions {βj : j = 1, 2, 3, ...} the weighted least squared problem is reformulated with
the objective function
Φ(m, t) =
1
2
‖Wd(Cm− d)‖2 + γt
Nβ∑
j=1
βj(m) (7.2)
The barrier functions used herein this work are antilogarithms, and physical constraints can
be implemented in a variety of ways based on the prior model of the flow. The weights γt are
means by which the boundary of the feasible region is buffered by a sharper and sharper buffer
as the iterations proceed; in the limit the buffer retracts to the boundary and the gradient of the
buffer is infinity at the boundary. Taking γt =
1
t is most common; as t goes to infinity the penalty
goes to infinity at the boundary still, but the region in which the penalty is non trivial shrinks
as 1t tends to zero. See Byrd et al. [2000] for a discussion of the interior point method used that
utilizes the log barrier technique. To use these methods with a non-local basis, as opposed to one
that tiles the domain, requires some manipulation. When the basis does not tile the domain, one
cannot just enforce the sign of the model parameters (i.e. the flow basis coefficients) to have the
same sign, or employ some form of non-negative least squares or non-positive least squares method,
to bound the solution. To bound the solution even by a single value of rotation rate requires an
alternate formulation.
From a grid of points, sufficiently fine that the basis cannot produce local features sig-
76
nificantly finer than the interpolating surface through those points of minimal curvature, basis
function values at those grid points for each basis function are stored in a matrix A, where Aij
corresponds to the value of basis function indexed by j at grid point index by i. An upper bound
vector is implemented as a column vector ub, and thus
Am ≤ ub. (7.3)
where the inequality is component by component. The lower bound is similarly implemented by
−Am ≤ −ul. (7.4)
That the fluid cannot super-rotate the outer shell seems a fair approximation for these purely
hydrodynamic scenarios during counter rotation, but a perhaps better bound would be obtained by
bounding the average specific angular momentum. Derivatives and other functions can be expressed
in this matrix form as well. Experimentation with many types of bounds was performed, such as
enforcing monotonicity in rotation along profiles in r, and enforcing monotonicity in specific angular
momentum along profiles in r. Such constraints only requires a suitable banded matrix to express
the values at the grid points as a function of the model parameters. These matrix formulations
allow many complicated flow constraints, but for this study bounding the extreme values seemed
the more easily justifiable, and where used extensively in chapter 8.
Each row of this matrix inequality form is used to form a log barrier function, implicitly,
e.g. βj = [Am]j− [ub]j for the first set of equations, and βj = −[Am]j +[ul]j for the other. In the
log formalism the Jacobian and Hessian can be formed analytically using matrix expressions or
evaluated using finite differences, and prevented from penetrating the boundary by backtracking
line search; if the boundaries are violated, the step in the search direction can be lessened until it
is within or on the boundary of the feasible set.
The system is solved iteratively using the Matlab optimization toolbox’s interior point con-
vex method first by way of the Matlab function lsqlin, then checked with alternate methods such as
77
the fmincon sequential quadratic programming algorithm, and fmincon interior point method just
described that utilizes log barrier functions; no noticeable differences were observed between the
two algorithms’ results though for applications in chapter 8. The interior point convex algorithm is
much faster, and allows ensembles of inversions to be performed efficiently to assess the sensitivity
of the inversion to data uncertainties. However, it is non-obvious formulated problem is convex,
or whether local minima ensures global optimization. Therefore, the slower interior point method
with log barrier functions that allows an ensemble of starting points for the iterations to be chosen
was used as well, to assess whether the solution is relatively unique.
Interior point convex algorithm uses a combination of Newton steps and conjugate gradient
steps to perform the iterative optimization. The Newton steps attempt to solve the system first by
finding a solution satisfying the Kuhn-Tucker conditions; very few iterations are typically needed in
these inversions. The solution is directed along what is called the central path toward the solution.
Ensembles of solutions are formed by taking random samplings of the error intervals around the
data as well.
Finally, to support that the iterative solution is converging to a global minima, compar-
ison is made with solutions from an ensemble of initial conditions by using fmincon in Matlab,
implementing the log barrier method. They were only checked for the presented solutions though,
and similar behavior was assumed for the ensemble that samples the data in the error intervals
around the measurements. The only evidence for multimodal distributions resulting from that
form of ensemble creation based on differing initial conditions came when significantly more basis
functions were used than could be accommodated by the resolving power of the data.
The initial conditions used in chapter 8 are approximately 100 random starting points
formed by the spline components of the basis. Each starting point satisfies the constraints, and
also visually identical results were obtained by using the alternate algorithm sequential quadratic
programming algorithm. See Gill et al. [1991] for a discussion of the active set method employed by
the sequential quadratic programming method. Only including nonzero coefficients for the splines
only in the initial conditions, is sufficient at the resolutions in chapter 8 since the mean velocity
78
field is expected to be largely vertically invariant.
7.2 Minimizing Spread Functions
7.2.1 Basic Principle
An alternative means of supplementing the data to obtain well-posed inverse problems is to
place restrictions on regularizations by way of how the resolution matrix diffuses basis functions
for the flow. In this section, it is beneficial to provide some elementary abstraction of the linear
algebra involved in inverse problems, for more easily understanding what is meant by the control
of diffusion. For a given set of basis flows B = [B1...BN ] define the model spaceM associated with
it as the set of coefficient vectors {m} such that any rotation profile Ω in the span of B can be
written in the form Bm for some m in model space. Let R be the range of the forward operator
C such that
C :M 7→ R (7.5)
Inverting C can be done many ways, but for exposition consider regularization by a damping matrix
D weighted by the damping coefficient . Define an inverse operator Q by the following formula
Q() = (CTWC + DTD)−1CTW (7.6)
where W is some form of weighting matrix such as an inverse data covariance matrix, or simply a
diagonal matrix with the reciprocals of the squares of the error weights for each data observation.
The inverse in the first factor need not be a traditional inverse, but could also include a pseudo
inverse obtained by truncated singular value decomposition or another method. Also,  has been
left as a scalar for simplicity of exposition, it could be more compplicated,e.g. a collection of
parameters and the damping terms. The model resolution matrix is obtained by forming an
endomorphism on model space by the map
Rm :M 7→M (7.7)
79
Rm, for each choice of , is formed by the product
Rm() = Q()C (7.8)
Due to both data uncertainty implicit in Q() and also resolving power of the kernels used
in forming Q(), Rm is typically not the identity matrix. For a basis B such as a grid of points,
tiling of the domain by disjoint convex polyhedra, or some other localized basis, each basis function
may not be mapped purely to itself by the inversion process. It is said that Q diffuses ej when
Rm maps m · ej to a vector in M that is not ej . Ideally, Rm would be the trivial isomorphism,
the identify map on M. However, even in that case, if the actual flow producing the data does
not reside in B, the accuracy of the inverse generated by Q can be highly inaccurate. To reduce
the ambiguity of selecting  and reduce the subjectivity of the choice of epsilon, various functions
can be defined on Rm().
Finally, it should be noted that many of these techniques, at least in publicly available
literature, seem to be used somewhat exclusively by disciplines that cannot readily verify their
solutions, such as helioseismic inversion for rotation in the sun’s interior, gravity inversions for
planetary bodies, and seismology of the deep earth and inversions for other geophysical variables.
This is due partly to the fact that empirically validated results give rise to the ability to establish
prior models of the expected data of higher accuracy and also heuristically base choices of regular-
ization parameters. Since there are still so many parameters implicit in the formulation that can
be said to be subjective though, it is not even true minimization of spread functions do not suffer
from confirmation bias when inverting data; they lift non-uniqueness by placing the subjectivity in
other choices. At best, they lift non-uniqueness for one specific regularization approach. They do
not eliminate the non-uniqueness across different inversions from differing regularization methods.
80
7.2.2 Dirichlet Spread Function Minimization
A simple functional to define on Rm to for the purposes described in the previous material,
is the Dirichlet spread function (DSF).
JD() = ‖Rm()− I‖F , (7.9)
where the norm is usual Frobenius matrix norm, the standard vector norm formed by unwrapping
the matrix into a vector, and I is the identity matrix with the same dimensions as the reoslution
matrix. Frobenius matrix norms can also be formulated more abstractly as the square root of
the trace of the product of the operator argument with its transpose/dual. Minimizing the DSF
means making the resolution matrix as close to the identity as possible subject to a least squares
constraint on the element wise difference. As can be seen, minimizing DSF penalizes equally any
off-diagonal terms. In terms of a basis that tiles the domain by small disjoint basis functions, that
amounts to penalizing any diffusion, but once diffusion occurs it does not penalize diffusion into
distant tiles any more than diffusion into adjacent neighbors. This approach often better preserves
finer scale features at the expense of creating overshoot and oscillation in the reconstruction.
7.2.3 Backus and Gilbert Spread Function Minimization
Weighting the matrix elements in the argument to the DSF’s norm differently depending on
distance between the basis functions is by far more popular than the DSF minimization approach.
That is, relaxing the penalty function when the diffusion is from nearby neighboring regions of the
domain. The Backus and Gilbert spread function is given by
JBG() = ‖w  (Rm()− I)‖F (7.10)
where w is some weighting of the distance between basis functions, often the squared Euclidean
norm between their centroids or some other interior point, and  is the Hadamard product. The
81
Hadamard product is the element-wise matrix multiplication. When the distance between two
basis functions is greater, diffusion between them during inversion is penalized more. Because any
basis function is zero distance from itself, diagonal have no effect in this formulation. When the
inversion produces pathological behavior such as zeros or very large or negative numbers on the
diagonal of the norm’s argument, additional constraints can be imposed to prevent that behavior,
such as enforcing conservation of the rows of the resolution matrix, i.e. making the averaging
kernels be a partition of unity or have length 1 or some other constraint.
Since the inversions in this work is on a meridional cross section, the spread matrix can
have its terms weighted by the actual volume, i.e. by weighting by (2piridAi)(2pirjdAj). Including
the ri factors, however, can degrade focus in the smaller values of r though. An argument can be
made for using a weighting by the specific angular momentum as well, as it may be more desirable
to conserve that quantity. With high amounts of data that has suitable coverage, it is often the
case that the spread is localized and the coefficients, such as r, are approximately equal within the
neighborhood of diffusion and it makes less difference what choice is made.
The parameters that optimize the spread functions can be solved using constrained opti-
mization. In the case of a vertically invariant basis such as vertical strips, a single regularization
parameter can be used and interval searches work as well. Additionally, other imposed constraints
such as forcing the solution to be non-positive, the rows of the resolution matrix to sum to one
or some other function, can be added as nonlinear constraints, solved either using the Matlab
optimization toolbox or with Lagrange multiplier methods. It should be noted though, that the
constraints related to enforcing conservation in the resolution matrix are not a constraint that was
satisfied when using iterative procedures as discussed in the previous chapter. It is also questionable
what the best conservation conditions for the rows should be, since unlike in helioseismology and
other fields there will be non-localization to a greater degree and perhaps the angular momentum
should be taken into account in different basis functions for rotation.
Since a large portion of the cavity is at rest, conservation in terms of a sum along the rows
is troublesome anyway. For instance, if a row has many nonzero but small elements, typically
82
one would be inclined to interpret that as indicating only small diffusion into a basis function
coefficient. However, some large regions are typically at rest whereas others are more locally
varying. An averaging kernel though, only tells the relative weighting of the flow throughout the
domain that diffuses into a basis function, not what the flow is across the domain. Such matters
are somewhat ignorable when the resolution matrix simply performs a slight Gaussian blur on
point sources, but for limited data such as used in chapter 8, the resolution matrix often transform
point sources into objects that are exceptionally non-convex, delocalized, not positive definite, and
have many other unfortunate properties.
For these reasons, these approaches were not used extensively, and not used as one of the
methods in chapter 8; though one example is given. The noble goals of the methods were difficult
to meet with limited data. Spread function methods had too many subjective choices that allowed
tuning the formulation of the inverse problem to allow the inferred flows to take on characteristics
that though plausible, were largely sculpted by the subjective choices, while still being afflicted
by artifacts. However, they did serve as an excellent source for creating flows that fit the acoustic
data, which facilitated understanding what features were actually required by the data.
7.3 Adaptive Parametrizations
The methods discussed in this section were only briefly examined, since they are sluggish
computationally. They were later realized to be somewhat feasible since only 15 modes were usable
from the provided acoustic data. Since it is still not clear how many modes are realistic to expect
from liquid sodium data, it seems relevant to mention these findings.
7.3.1 Adaptive Interpolation
One way to overcome the issues with regularization that allows the parametrization to take
full advantage of the data without fitting noise or the downsides of the coverage problems, is to use
an adaptive parametrization. In this section, I introduce the idea of fitting the flow without using
the forward operator, but only the kernels and the adapted flow reconstruction. The basic idea is
83
to interpolate a solution from control points in the domain and on the boundary of the domain,
and instead of just performing constrained optimization on the values of rotational frequency at
the control points, also allow the control points to move in terms of location. This is a much
more computationally costly method. The objective function used is the same measure of misfit
as weighted least squares, but the model parameters are no longer simply coefficients of a basis.
To increase the speed, the fit is done in several stages. The initial state is a combination of 15
grid points just outside the boundary of the top half of the meridional cross section, roughly equally
spaced along the boundary, creating a set of points whose convex hull contains all of the domain
when they are replicated to the bottom half of the cross section; that is because the function is
being presumed to be even about z=0. These points do not move. An additional 8 points roughly
equally spaced out inside the top half of the cross section are used as free moving points. The
model parameters are the 23 values of rotation at each grid point, and the 8 x coordinates of the
moving points and the 8 y coordinates of the moving points. They are further constrained by
preventing positive rotation values, and the moving points are confined to being inside the box
[0, Ro] × [0, Ro], strictly inside the domain. Everything is replicated to it’s mirror image so as to
enforce even solution about z=0. The interior point method described earlier in this chapter is
used to enforce the constraints. The number of moving points on the interior of the domain was
chosen to be the maximum value that could be used without very small clusters of points being
created. Using more can often allow the moving points to form clusters so that the solution can
produced very high flow in a very small region that few modes are sensitive to, thereby fixing any
misfit their have with the bulk flow, while producing highly unrealistic features, especially at the
resolution expected form the data.
The first stage of adaptation produces a solution via nearest neighbor interpolation of the
values at the grid points, by a Voronoi tessellation of the domain. What is noteworthy, is that it
is much easier to fit the data using Voronoi tessellations than smooth interpolants, and many of
the tiled flows over the meridional cross sections, seem to recover plausible placements of the high
shear zone, as well as recovery of the jet instability where the inversions recovered it as well. It
84
suggests that perhaps high shear is more consistent with the data than smooth diffuse gradients
in angular velocity.
The computational speeds on a relatively fast desktop computer, can amount to an hour
or longer for an inversion if no prior is used and Hessians are computed numerically to find the
search directions; using Broyden-Fletcher-Shannon-Goldfarb algorithm or similar means. It can
be sped up by performing a preliminary inversion using one of the earlier described methods, then
writing with respect to the interpolant formulation as accurately as possible; for instance, selecting
points in the domain manually, perhaps where there is little posterior predictions in variance based
on data uncertainty. Otherwise part of the sluggishness comes from flow values changing as grid
positions change. Once the grid points stop moving as fast, the effects of changing the flow values
at each grid point, less backtracking is needed, and hence less function evaluations.
7.3.2 Non-coefficient Model Parameters
Often it more convenient to parametrize using means other than basis functions. For in-
stance, the more defining feature of a flow might be a region of shear and the gradient of that
shear. Consider a sigmoid function
S(r;A, rc, k) = A− A
1 + e−k(r−rc)
(7.11)
rearranging
S(r;A, rc, k) = A
e−k(r−rc)
1 + e−k(r−rc)
(7.12)
In this context, A is the amplitude of the lower rotation rate near the inner sphere and shaft, rc
is the cylindrically radial location of the shear zone, and k controls the gradient of the shear.
It is difficult to assign an uncertainty to this quantity, since the majority of the error stems
more so from the idealization of the flow rather than variations in inversions due to data uncertainty.
In chapter 8 this formulation is used by fixing k to be very large and estimating the location
of the shear zone. What has been found by experimentation with synthetics, is the the inversion
85
that is two essentially piecewise constant and disjoint basis functions who share a boundary at
some value r, will place the border where the shear zone is, since the data misfit incurred by not
fitting the larger volume piece by zero where it is relatively moving with the outer sphere is much
larger than averaging over the portion of the cavity at values of r less than the border value of r.
The data misfits are also quite small, having larger misfits where 2D features become prominent,
such as a jet-like coming off the inner sphere equator at roughly Ro=-4.67.
Such 2D features tend to defocus the shear when the value of k is allowed to vary, even
though the 2D inversion profile for rotation shows a relatively vertically invariant flow excluding
the jet-like feature and is no more diffuse in shear than other rotation rates; see chapter 8. This
is consistent with the adaptive parametrization approach using moving grid points. The grid
points do not actually spread out in the domain, but rather cluster along a the shear zone so that
they can create shear when interpolation is carried out on their adapted fluid rotation values. A
cylindrically symmetric basis was also used that constrained the specific angular momentum to
being monotonically increasing found similar results, consistent with the high shear zone. It was
not forced to taper to a constant rotation rate on values of r less than that at the shear zone, yet
recovered that feature as well.
86
Chapter 8: Applications to Ambient Noise Data
With only 15 usable modes, and not 15 specifically chosen modes, extracting useful robust
information about the flows is challenging; especially using ambient noise as the acoustic source.
In this chapter, inversions are performed, and discussion of quality of fit to data, and reliability of
the reconstructed flow are presented.
8.1 Data Processing
8.1.1 Array Configuration
The modes used for inversions had stationary frequencies up to approximately 2.1 kHz; see
table 8.1.1. The microphone configuration did not allow |ms| to be identified using array techniques,
since they were all separated by pi with respect to the azimuthal angle. The configuration could
not be easily changed, I was informed. Another microphone was offset in the polar direction to be
able to perhaps distinguish via l− |ms|, however only two microphones could be used at one time.
An additional microphone had been placed at an offset from another by roughly seven degrees
in the azimuthal direction, in the hope of being able to distinguish between values of |ms| and
increase the data set size. The small angle would have allowed many more modes to be identified,
but was perhaps overly ambitious in the small spacing; the small separation and combination of
the effects seen in the signals made difficulties in using the same cross correlation technique used
in Triana et al. [2014].
The microphones separated by 180 degrees in azimuth were used in the data. The micro-
phones were at 45 degrees in terms of colatitude. The explanation given, being that sum and
difference techniques were able to remove a great deal of the noise; specifically, odd valued |ms|
87
modes would be greatly amplified by taking the difference between the two microphone signals,
and since presumably digital noise and mechanical noise would be the same on each channel, would
be virtually removed in the result. Odd valued |ms| modes have equal magnitude but opposite sign
at locations separated by 180 degrees in azimuth; hence taking the difference doubles their signal
level while removing any features that are the same in both signals. The microphones were only
calibrated insofar as the same level was recorded at each location; the frequency response was not
calibrated. Calibration of the frequency response would likely have been of little help with ambi-
ent noise data, since the strength of the source term distribution is unknown. That is, knowing
the loudness of each mode at the microphone position is of limited value if the source terms are
not equal as well. Relative amplitudes should be equal at equal azimuthal angle for the mode as
well. These limitations of the array geometry made the FEM predictions all the more important,
since the array techniques available could only distinguish between evenness and oddness of |ms|,
placing additional emphasis on knowing the correct frequency range to expect modes to appear
in. The mode identifications were thus done by evenness and oddness and frequency location. In
all but one instance, the modes could be unambiguously identified based on that criteria. Mode
(1,3,2) was slightly ambiguous, but was distinguished based on its splitting behavior and consis-
tency with the flow suggested by the other modes observed. Although, the identification analysis
revealed the inversions were robust to removing this mode from the set, and hence the mode was
not providing any qualitative difference or visible addition of information. This would likely not
be the case if the residual frequency shifts were larger relative to the error bars associated with
the measurements of the spectra. It was tempting to try and identify other modes based on the
inversions and then increase resolution by including them. However, it seemed too likely to create
a gross error though, with no benefit. If a splitting was very accurately predicted by the inversion
with fewer modes used, then including the other modes that matched that inversion’s prediction
did not really have any added benefit. When the spectral location did not match the inversion
with fewer modes sufficiently for that to occur, and possibly presented new information, there was
still too much ambiguity to have any comfort in the accuracy of the identification. The benefits
88
index n l m f (Hz) δf6 (Hz)
1 0 1 1 346 4.0
2 0 2 1 570 1.3
3 0 3 1 774 0.4
4 0 3 2 807 2.8
5 0 4 1 970 -0.0
6 0 4 3 1014 3.0
7 1 1 1 1045 -0.8
8 0 5 1 1165 -0.2
9 0 6 1 1361 -0.2
10 1 3 2 1500 -0.2
11 1 4 1 1702 -0.7
12 2 2 1 1898 -0.6
13 1 5 1 1938 -0.5
14 2 3 1 2060 -0.8
15 2 3 2 2099 -0.5
Table 8.1: Mode wavenumbers and stationary frequency prediction at c=344 m/s. These are the
modes used to perform the inversions. Additionally, the frequency shift from center frequency due
to uniform rotation with the outer shell at 6Hz has been provided, rounded to 1 decimal place,
since it does not appear in the presented splitting data, it has been subtracted off. The negative
sign indicated that −|ms| has a higher frequency than +|ms|. Values rounded to nearest tenth of
a Hz. The splitting is thus approximately twice these listed values. The first column provides the
index used in plots of data fits in this chapter, for reference. They are always presented in this
order herein.
were too little to risk corrupting the solution with incorrectly identified modes.
Unfortunately, the majority of the data for 60cm setup is ambient noise based. That is,
the resonances are either excited by mechanical vibration or the flow of the air. The amplitude
of the signal in each frequency band was not constant across within every frequency band. This
means the spectrum in each band is convolved with the spectrum of the amplitude envelope of the
band. For instance, for a simple sinusoidal amplitude envelope in the time domain that amounts to
convolving with two delta spikes about the origin in the frequency domain, as discussed in chapter
1.
Signals were decimated by a factor of 6, and spectra were computed using prolate spheroidal
slepians with a time bandwidth product of 2.5; a window was chosen to smooth the clusters in
typical peaks that would be obtained by a steady resonance. The widths of the spectral peaks is
not due purely to defocussing of a single peak though, but also the smearing together of the cluster
presumably. Given the noise, the affects of the flowing air, turbulence, computing the splitting
was not as simple as locating two peaks and measuring the frequency separation, typically. For
89
each peak, I chose two points on both sides of the localization of frequency content contributing to
a mode. The mean location of the mode resonance inside that band was often within 1/2 radius
around the center of those measurements, but it could not be said that there were not fluctuations
beyond that interval for short instants. See figure 8.1. The shift from center frequency being half
of the splitting inferred from the measurements at corresponding |ms| pair. There is good reason
to suspect that to a significant extent, the pair moved in an either correlated or anti-correlated
fashion. As such, resolving power was not removed from a basis if the basis could fit the data
by 1/4 of the average error radius associated with the data. It does not mean the inversion is
over-fitting the mean angular velocity field, only possibly misfitting any one particular flow state
during the signal acquisition.
There is not a universal rule for calculating the splittings employed though, as all were not
as simple to measure as figure 8.1. The guiding principle was accurately gauging the mean splitting
to within as small an error radius as possible. If a splitting could be reliably constrained, it was
measured. For instance, often digital noise overlapped the spectra forcing the estimate of the flanks
of the spectra locations to be expanded, being conservative. In some instances, the mode peaks
were sufficiently clear, but the splitting was so small that not only was the overlap of the modes
taken under consideration in judging their mean spectral peak locations, but also the fact that it
was not completely obvious which mode in the mode pair |ms| was on the left and which was on
the right; that effectively doubled the error radius. That meant that often the average data misfit
was not a good indicator by itself, and the fitting across the mode set was very informative to
analyze.
90
Figure 8.1: Comparison of acoustic spectrum with different frequency resolutions. The black curve
has been created with a frequency resolution of 0.3 Hz whereas the red curve has a frequency
resolution of 0.05 Hz. Note: magnitudes are linear, not in decibels. The clusters of peaks could
be due to poloidal flow, noise envelopes, mechanical noise, changing mean flow state during time
integration, or some combination thereof.
91
As can be seen from the red curve in figure 8.1, it is useful to have error bars that encompass
values of splitting due to changes in flow as well. Though it could be the affects of noise envelopes
or transient poloidal flow changes, there could be changes in zonal flow state as well. The actual
estimate on the error of the mean is smaller than the error that encompasses possible changes
in the flow state. For the inversions this makes no significant different. Since a uniform scaling
only has effects on inversions when the data strongly dislike the constraints and trade-offs in data
misfit tolerances with constraint violation tolerances cause it to violate the constraints slightly.
A future improvement could include a hierarchical approach to different types of errors, but a
conservative approach was taken in this work since it is not completely clear what is affecting
the distribution inherent in those peaks, in terms of noise, different levels of excitation, different
microphone response, changing flow, quality factors of the modes, and turbulence.
To facilitate identifying the modes amidst the noise and other aspects of the signal, the
spectra were formed into an intensity plot with acoustic frequency as horizontal axis and inner
sphere rotation rate as vertical axis; see figure 8.2 for an example. In these spectrograms, the
dynamics of the peak locations could then be tracked over the range of rotation rates. Even with
the drifts due to thermal fluctuation as the experiments proceeded, there is a relatively continuous
trend in each peak. Thus alleviating any uncertainties about which peak is acoustic and which is
digital noise or another mode. The stratregy of measuring the splitting was then performed with
the signal plotted below these spectrograms, so that the signal could be analyzed in concert with
these spectrograms; see figure 8.5.
The first stage of processing the splitting data is assign a polarity to the splitting. The
measured splittings do not give the polarity of the splitting, whether the +|ms| or the −|ms| mode
is the higher or lower frequency. Sum and difference techniques were used to infer the even and
odd parity of the value of |ms|, but the polarity had to be inferred based on physical assumptions
of the flow.
For inner only rotating and the microphone in the outer shell rotation rates of the fluid with
the same sign as the inner sphere’s rotation rate can be presumed. Once the outer shell is moving
92
though, the effect of advection on the spectra in the majority of the cavity is greatly diminished,
and the majority of the volume of fluid is nearly at rest in the rotating frame for many Rossby
numbers; the splitting measured there is more so dominated by the Coriolis induced splitting. For
instance, if a splitting is measured to be 2 Hz in the rotating frame, and the prediction of uniform
rotation is a Coriolis induced splitting of -1.4 Hz, there is ambiguity if the splitting due to flow
relative to the rotating frame is an additional 3.4 Hz, or if it is only an additional 0.6 Hz, since
the original measurement could have actually corresponded to -2 Hz. Since data was acquired
at a broad range of rotation rates, for instance the outer spinning at 6 Hz while the inner was
spinning at -4, -6,...,-40 Hz, it is possible to track the splitting from a regime where the sign can be
safely presumed back to a regime where it would be more ambiguous. This places more emphasis
on the physical constraints used in the inversion, since the ill-posedness of the problem in their
absence would allow any polarity to be fit, at the expense of producing a physically unrealistic
reconstructed flow. The next stage of processing of the splitting data is to remove the effect of the
rotating frame.
The Coriolis induced splitting due to uniform rotation is subtracted from the data, and
then the inversion is done in the rotating frame as if it were at rest. This does not presume the
nominal value of the rotation rate is correct, but the flow rotation reconstructed is relative to that
assumption; the inversion is relative to a constant rotation frame at the nominal value of rotation.
For inversions that place upper bounds on the rotation rate, the expected average deviation from
the nominal rotation rate is not sufficiently large to warrant discussion of how a tiny fraction of
a Hz error int he bound affects the inversions, since they are not at an accuracy level it would be
relevant. There is a tolerance of a several hundredths of a Hz allowed also, and testing was done
to show that small fluctuations in the upper bound did produce a visually dissimilar answer; often
that bound could be substantially raised by several Hz while often only incurring small oscillations
above the upper rotation rate limit int he profiles. If the flow was sub-rotating relative to the
outer shell slightly, this work does address whether that is due to the fluid rotation rate being
slower, or if the nominal value of rotation is slightly inaccurate; a precision motor was used, and
93
the deviation is claimed to be negligible for these pruposes [Lathrop, personal communication].
Figure 8.2 shows an intensity plot for power at each frequency with axis representing acoustic
frequency and inner rotation rate frequency for the mode (0,1,1), as recorded in the rotating frame.
Microphones at φ = 0 and φ = 180 degrees and both at θ = 45 degrees were used, and their
difference taken so as to avoid noise and modes with even m. The mode (0,1,0) is absent in
the plot as a result of the difference being taken of the two microphones. The sum of the two
microphones is plotted in figure 8.3.
94
Figure 8.2: Power spectra plot of acoustic frequency(horizontal axis) vs. inner core rotation rate
( vertical axis ) near the modes (0,1,1) and (0,1,-1). Color is relative amplitude in log scale. The
spectra is computed by taking the difference at two microphones at equal latitude spaced apart
by 180 degrees in terms of azimuth. This shows drift as the temperature increases inducing center
frequency shift, and shows that the dominant effect on splitting is the outer spinning at 6 Hz rather
then the greater splitting that occurs due to additional rotation of the inner core. There could
also be some effect of poloidal flow. In this regime centrifugal force is not expected to have any
significant effects.
95
Figure 8.3: Similar figure and axes as in figure 8.2. Power spectra plot of power spectral density
vs. inner core rotation rate near the mode(0,1,0). This shows drift as the temperature increases
inducing center frequency shift. Color is relative amplitude in log scale. The spectra is computed
by taking the sum of two microphones at equal latitude spaced apart by 180 degrees in terms of
azimuth.
96
Figure 8.4: Power spectra plot of frequency vs. inner core rotation rate near the mode (2,3,1).
Intensity of the plot is power spectral density, relative amplitude in log-scale. Note splitting
changes with rotation rate, and also the drift of the center frequency.
97
(0,1,1) and (0,1,-1) are the lowest frequency acoustic modes of the resonant cavity, and their
splitting is dominated by the Coriolis effect due to being in the rotating frame in the regimes
of these experiments. The additional shifts due to the shear flows are apparent, but much less
pronounced; the drift due to increasing temperature associated with spinning the inner core faster
and the fact that the experiment had been running longer is more obvious. Plots such figure
8.5 allow the trend of the mode peaks to be seen in the spectra, since when looking at a single
inner core rotation rate, ambiguity can arise as to what the actual acoustic peak is and what is a
mechanical noise or resonance.
For instance in figure 8.4, which shows the difference power spectra between two micro-
phones, centered near (2,3,—1—) in the spectrum, the signal to noise ratio is not as high and noise
is more apparent.
Figure 8.4 also shows that the effect of the inner core is more significant for (2,3,1). The
splitting is nearly 0 when the Rossby number is low; it is offsetting the effect induced by the
Coriolis force to split the pair negatively. The error bars are large for this mode when the Rossby
number is nearly zero, since the peaks cannot be separated well. The splitting can be tracked in
figure 8.4 back to where that splitting of zero occurs. For a inner core rotation rate of +6 Hz, the
expected splitting is that of uniform rotation, but between inner core rotation rates of +6 Hz and
0 Hz, without making assumptions about the flow pattern, it is impossible to judge the sign of the
splitting. This sign ambiguity was tested by assuming a geostrophic profile and a monotonic profile
between the axle and where the expected region that is executing solid body rotation occurs, and
found the sign was too variable in regard to profile input.
The plotting of the profiles for each rotation rate was very useful in more ambiguous cases
such as the one illustrated in figure 8.5, which shows power spectra in close vicinity to digital noise
peaks in the spectrum.
98
Figure 8.5: Power spectra plot (difference between microphones) of power spectral density vs.
inner core rotation rate near the modes (0,2,1) and (0,2,-1). This shows drift as the temperature
increases inducing center frequency shift, and shows that the dominant effect on splitting is the
outer spinning at 6 Hz rather then the greater splitting that occurs due to additional rotation
of the inner core. The lower plot is a power spectral density estimate placed on a linear scale,
for inner core rotation rate of -40 Hz. It is a relative amplitude, not calibrated to any particular
loudness level.
99
Figure 8.5 shows a double peak around the second mode of the pair. The spectral plot
shows a periodic noise feature in the range of the spectra that overlaps the two modes’ spectra; 4
columns can be seen in the plot evenly spaced, due to noise. The acoustic peaks have a different
characteristic shape, but the plots make it much more obvious what is the acoustic signal and what
is the electronic noise. It also shows how in the left peak, there is overlap, which causes a slightly
increased amount of uncertainty placed on the peak location. It is not clear that the peak is the
best estimate of the mean location of the acoustic spectral energy.
To ensure that visible drift of the center frequency is explainable entirely by thermal in-
duced uniform soundspeed increases associated with the experiment running, the frequencies of the
ms = 0 modes were calibrated to the highest observed ms = 0 mode (1,6,0) to test their relative
invariance; each observed ms = 0 mode’s observed frequency was multiplied by the ratio of the
original starting center frequency associated with (1,6,0) to the observed (1,6,0) frequency at the
time of measurement. This follows from assuming the wavenumber k such that ω = kc is invariant
then the calibrated frequency
fcalibrated = ft · f
(1,6,0)
o
f
(1,6,0)
t
(8.1)
should remain invariant as well, where fo is the frequency at some reference state soundspeed, and
ft is the frequency at the current observation. This is equivalent to
ft = fo
ct
co
(8.2)
where c is the sound-speed defined with the same subscript notation as for frequency.
Once normalized this way, the drift of the ms = 0 modes all had a standard deviation of
0.33 Hz or less, with only a few outliers that were easily explained by the data uncertainty.
100
8.2 Parametrization
8.2.1 Approach to Parametrizing Fluid Rotation
The rotational profile is parametrized two different ways and inversions are performed. One
method utilizes a basis that caters to expected flows, the other being more generic and more
localized.
Also, the constrained optimization approach bounds the solution, which can strain the side
lobes of modal type basis functions that prefer to oscillate around the answer rather than be
bounded by it; in regions of the cavity near uniform rotation with the outer shell, the ripple being
bounded below that rotation rate could distort the answer slightly more than necessary. The
bound related to the inner sphere did not seem to have significant effect so long as a reasonable
rate was chosen, between one half and the rate of the inner sphere. The locations that seem to be
more affected by the specific bound, are inside the missing pieces of the shaft; a region that little
significance should be attached to with respect to the reliability of the inversions already.
Additionally, the inclusion of cubic splines that are only a function of cylindrical radius
strongly supports flows that have rotation profiles that are strongly cylindrically symmetric; i.e.
vertically invariant to a marked degree. However, when the vertical heterogeneity of the profile
is highly localized, smearing into the splines occurs as a result of them being able to fit the data
better than the modal basis components; many of the regions of the cavity have little sensitivity.
Comparison with a more localized basis offers the chance to see if the inversion produced that way
has contiguous vertically invariant components or if they are possibly an illusion due to smearing
inherent with the lossy inversion.
The most difficult aspect of using such a small data set, is that it is quite easy to fit the
data with any number of profiles usually; more advanced techniques to regularize the problem
are difficult without introducing subjective parameters and ad hoc justifications. The approach
adopted with the limited data, is to provide a tool for laboratory experiments which generate
inversions that robust features can be assessed from, and combined with other observations provide
101
a tool for improving methods to explore dynamo theory; not necessarily provide rigorous evidence
for the recovered profiles’ accuracy. With So few spectra and difficult to interpret error bars
incurred by being relegated to ambient noise data, providing quantification to the recovered profiles
for very localized regions of flow is not recommended.
The resolution chosen, i.e. the basis chosen, was also to prevent over-fitting the data. The
goal of having misfits between 1/4 and 1 as the average misfit was chosen for two reasons: 1) it
prevents outliers with respect to the error intervals assigned, and 2) the actual mean flow splitting
is likely to be as low as approximately one half or less of the estimate, since the majority of the
error comes from clustering of peaks around each spectra; the +|m| versus −|m| pair are highly
correlated in their relative motion in the spectrum. An error of 1 would correspond to completely
uncorrelated errors in +|ms| and −|ms| modes, 1/2 would correspond to assuming correlation
of their mean location, and 1/4 would correspond to further assuming that their mean is more
accurate than the actual measurements. What is meant by the last statement, is that the way the
way the splittings were originally measured attempted to capture the variation about the mean,
but the mean is much more accurate. As a demonstrative example, suppose there was a 2-peak
cluster at 1001 Hz and 1003 Hz and it’s mode pair of the same (n,l,—m—) had a corresponding
2-peak cluster at 1011 Hz and 1013 Hz. The measurement would grab the the point on the left
of 1001 Hz at FWHM of the 1001 Hz peak, and at FWHM to the right of the peak at 1003 Hz.
The location of the average is then 1002±FWHM, and similarly for the right 2-peak cluster. The
shift due to azimuthal flow is half the splitting.
Even if changing flow state anticorrelates which peaks correspond to each other from cluster
to cluster, their mean splitting over the time acquisition is still approximated by the mean sepa-
ration of the cluster centers; though it could be a splitting value the actual flow does not produce
for any significant amount of time, but the mean flow over a large time interval might not ever
be realized by the fluid as a mean flow over shorter time interval either. Most often, as was the
case for the example in figure 8.5, the mean of the peak cluster was more accurate than the actual
distribution could be approximated by it’s mean. The error on the mean was due to them not
102
being delta spikes, not because the mean was not known more accurately. Since there were so many
different types of errors though, it seemed more prudent to present them using this scaling, rather
than rescale to this interpretation. For instance, at some of the low rotation rates, the splitting
was so small that it is difficult to judge the polarity of the modes, i.e. which member of (n, l, |ms|)
is on the left and which is on the right in the spectrum. The error bar had to be larger not because
it was difficult to measure the small splitting δf , but because the actual splitting could be −δf .
8.2.2 Cylindrically Radial Cubic B-Splines and Neumann Modes
The flow is parametrized by a combination of cubic splines in the cylindrical radius param-
eter and Neumann modes of the spherical annulus, the total number of basis functions numbering
15, the total number of splittings used; the modes had stationary frequencies up to approximately
2.1 kHz.
Six equally spaced cubic splines with two on the endpoints of the range of cylindrical radius
values having their peaks at the outer boundary and at the shaft radius were used. The cubic
splines help to parametrize the cylindrically symmetric component of the flow. They are essential
for slow flows relative to the outer sphere’s rotation rate, which are approximately geostrophic in
much of the cavity, and for the capturing the cylindrically symmetric component of the flow at
higher flow rates. The splines in cylindrical radius r cannot create a very sharp shear at the tangent
cylinder, it must be somewhat defocused. Using more splines so that such an effect is achievable,
runs the risk of over fitting the data and decreasing robustness in favor of higher resolution. It is
however possible, to simply add a vertical rectangle function there to accommodate shear at the
single region, but in this study care has been taken to prevent as much as possible the creation
of features in the inversion that are idiosyncrasies of the particular basis chosen, and so that has
not been done herein; the way it is done here, may not be ideal for inverting when there is very
strong shear at the tangent cylinder, but when the inversion stratifies the flow into obvious regimes
inside and out of the tangent cylinder, it is less likely an artifact of the parametrization; and the
interpretation of defocussing of the shear at the tangent cylinder by the inversion requires user
103
interpretation and judgment when making decisions.
Deviations from vertical invariance were parametrized by the products of zonal harmonics
with the radial wave functions of the acoustic spherical annulus; the spherical Neumann modes
of the spherical annulus. These Neumann modes were all the combinations of the first 3 values
of n for even values of l=0,2,4. They were computed using the same method as described in the
Numeric Methods section. The radial functions were computed by FEM and the polar component
by using zonal harmonics. The actual modes of the cavity were not used since those are not strictly
even about θ = 0 due to the asymmetry of the missing pieces; though I did use that approach
in earlier work. Even in absence of the asymmetry, though, vanishing Neumann conditions are
not ideal along the shaft anyway, and any flatness along the shaft of the 2D profile is well com-
pensated by the cylindrical spline overlapping it. For much higher resolution inversions, using the
modes of a modified cavity whereby the missing pieces on the bottom and top are made equal to
ensure symmetry could be preferable. Odd values of degree l were not included, since they yield
components that are approximately in the nullspace, and slight deviations due to the asymmetries
due to the missing pieces seemed unlikely to be captured accurately with low-resolution inversions
such as these. This means that the reconstructed flow around the missing piece in the bottom
hemisphere may not be very reliable at all, since it is a composite effect of assuming the flow is even
about z = 0 when the missing pieces are slightly different size in the top vs. bottom hemispheres.
All inversions except where otherwise stated, use this parametrization. Other bases consisting of
splines or polynomials or indicator functions have been used to assess bias introduced by the bias,
and are noted directly throughout. The bias introduced by a basis does not only stem from its
span, but also how that span relates to the sensitivity kernels of the modes used to perform the
inversions.
It is important to note that using low-resolution basis increases the robustness of the inver-
sion to data uncertainty, but it does not actually increase the confidence in the inversion unless one
presumes that the flow is virtually devoid of variations at length scales that cannot be resolved by
the flow basis. Using a low-resolution basis actually magnifies the errors with respect to small-scale
104
variations in favor of capturing broad averages more reliably, presumably. To prevent bias caused
by the parametrization, alternate parametrization are presented as well, to show the manner in
which low-resolution reconstructions can defocus flow differently, as a means of gauging what are
the actually robust and persistent features across inversions.
8.2.3 Choice of Reference Frame
The inversion was performed in the rotating frame by removing the effect of uniform rotation
Coriolis induced splitting from the data using first order perturbation theory prediction, then
performing the inversion as if at rest. The outer sphere is rotating at 6 Hz; the inner sphere is
rotating at -4 Hz for the low rotation rate inversion and at -22 Hz for the next rotation rate focused
on, and -40 Hz for the highest magnitude inner rotation rate inversion; intermediate values are also
used for inversions and analysis. The solution was solved by constraining the values of rotation
at any location to be between the inner and outer rotation rates. Regularization is implicit due
to the low-resolution basis used to accommodate the low amount of usable data and resolving
power of the spectra observed in terms of sensitivity. Part of this work is developing means of
implementing such types of constraints and assessing their influence on solutions, more so than
the specific applications in this chapter.
8.3 Processed Data
In this section I present the data used throughout the rest of the chapter. Three tables, 8.2
8.3 and 8.4 give rounded values for the data used in this chapter.
105
(n, l, |ms|) -4 -6 -8 -10 -12 -14 -16 -18 20 -22 -24 -26 -28 -30 -32 -34 -36 -38 -40
(0,1,1) 0.3 0.4 0.4 0.5 0.6 0.6 0.6 0.5 0.5 0.7 0.6 0.8 1.0 1.1 1.0 1.1 0.7 1.2 1.3
(0,2,1) 0.5 0.7 0.6 0.6 0.7 0.7 0.7 1.0 1.2 1.2 1.5 1.6 1.6 1.7 1.8 1.7 2.1 2.4 2.3
(0,3,1) 0.7 0.5 1.1 1.1 0.8 1.3 1.5 1.7 1.6 1.4 2.0 2.4 2.5 2.6 2.7 2.6 2.8 2.6 2.7
(0,3,2) 0.3 0.6 0.7 0.9 1.1 1.3 1.4 1.8 1.7 2.0 2.2 2.5 2.8 2.9 3.1 3.4 3.7 3.9 4.4
(0,4,1) 1.0 1.2 1.1 1.4 1.6 1.7 1.9 2.2 2.1 2.5 2.7 2.8 3.0 3.2 3.4 3.5 3.4 3.4 3.6
(0,4,3) 0.2 0.3 0.5 0.7 1.0 1.2 1.4 1.6 1.8 2.3 2.7 3.1 3.5 3.9 4.2 4.6 5.1 5.5 5.9
(1,1,1) 0.8 0.8 0.8 0.9 1.3 1.5 1.8 2.0 2.4 2.6 2.9 3.3 3.5 3.7 4.0 4.2 4.3 4.5 4.7
(0,5,1) 0.5 0.6 1.2 1.2 1.6 2.0 1.8 2.1 2.0 2.1 2.8 2.6 2.6 2.8 3.0 3.2 3.5 3.6 3.5
(0,6,1) 0.8 1.1 1.6 1.2 1.4 1.8 1.8 2.0 2.4 2.5 2.6 2.7 2.9 2.9 3.0 3.1 2.9 3.2 3.3
(1,3,2) 1.2 1.9 2.4 2.5 2.9 3.5 4.0 4.8 5.2 5.7 6.0 6.5 7.0 7.4 7.4 7.9 8.1 8.3 8.5
(1,4,1) 0.8 2.5 2.8 2.6 3.0 3.1 3.6 3.6 3.8 3.8 4.3 4.2 4.8 4.7 4.6 4.3 4.9 5.3 5.1
(2,2,1) 1.2 1.5 1.7 2.2 2.5 2.6 2.9 3.4 3.6 3.6 4.3 3.9 4.9 4.5 4.6 4.5 4.9 4.3 4.4
(1,5,1) 1.5 2.2 1.8 1.9 2.0 2.0 3.0 3.5 3.4 4.1 3.8 3.9 4.3 4.2 4.6 4.8 4.4 4.1 3.7
(2,3,1) 2.2 2.1 2.5 2.7 2.9 3.2 3.6 3.8 4.1 4.3 4.5 4.5 4.5 5.0 5.2 5.2 5.1 5.1 5.1
(2,3,2) 1.4 1.9 2.6 3.2 3.5 4.2 5.5 5.2 6.1 6.7 6.9 7.3 7.7 8.3 8.2 8.4 8.6 9.2 9.6
Table 8.2: Frequency shift from center frequency in Hz. Top row is Fi in Hz, and each row is a mode pair given by the first column. Measurements
are in Hz. Outer Sphere is rotating at 6 Hz and the inner sphere rotation in Hz is given at the top of each column. Each row is for one particular
mode. The data has had the correction due to being in a rotating frame removed from it. For instance, the first row shows (0,1,1) has relatively
small shifts, but its shift due to the Coriolis effect of being in a rotating frame is almost 4 Hz. The value stated is the residual after that amount was
subtracted off. Refer to table 8.1 for the approximate Coriolis correction of the rotating frame that has been removed.
106
(n,l,|m|) -4 -6 -8 -10 -12 -14 -16 -18 20 -22 -24 -26 -28 -30 -32 -34 -36 -38 -40
(0,1,1) 0.3 0.2 0.3 0.3 0.2 0.4 0.4 0.2 0.2 0.2 0.3 0.2 0.2 0.3 0.2 0.2 0.3 0.3 0.2
(0,2,1) 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.4 0.3 0.2 0.2 0.1 0.1 0.2
(0,3,1) 0.2 0.2 0.2 0.2 0.2 0.2 0.3 0.2 0.2 0.2 0.2 0.3 0.4 0.3 0.2 0.2 0.2 0.2 0.2
(0,3,2) 0.3 0.5 0.4 0.3 0.3 0.4 0.4 0.4 0.3 0.3 0.3 0.3 0.3 0.5 0.4 0.4 0.3 0.4 0.2
(0,4,1) 0.3 0.3 0.2 0.3 0.4 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
(0,4,3) 0.3 0.3 0.3 0.3 0.3 0.2 0.2 0.3 0.3 0.3 0.3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
(1,1,1) 0.3 0.2 0.2 0.2 0.2 0.1 0.1 0.2 0.2 0.3 0.2 0.2 0.2 0.2 0.2 0.3 0.3 0.3 0.4
(0,5,1) 0.6 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.3 0.3 0.5 0.5 0.5 0.6 0.5 0.5 0.4 0.5 0.4
(0,6,1) 0.5 0.4 0.5 0.4 0.5 0.5 0.5 0.6 0.5 0.5 0.5 0.6 0.6 0.6 0.7 0.6 0.6 0.6 0.5
(1,3,2) 0.3 0.3 0.3 0.3 0.3 0.2 0.3 0.3 0.3 0.3 0.3 0.5 0.4 0.4 0.5 0.5 0.5 0.5 0.4
(1,4,1) 0.7 0.4 0.4 0.4 0.5 0.5 0.5 0.5 0.6 0.7 0.6 0.8 0.9 0.8 0.9 0.7 0.6 0.7 0.6
(2,2,1) 0.4 0.3 0.3 0.4 0.3 0.6 0.3 0.3 0.3 0.4 0.5 0.3 0.3 0.2 0.3 0.2 0.3 0.4 0.4
(1,5,1) 0.5 0.8 0.7 0.8 0.6 0.7 0.6 0.6 0.5 0.7 0.7 0.5 0.4 0.7 0.8 0.6 0.6 0.5 0.3
(2,3,1) 0.5 0.5 0.5 0.5 0.6 0.5 0.5 0.6 0.5 0.5 0.4 0.5 0.5 0.5 0.6 0.5 0.8 0.5 0.6
(2,3,2) 0.4 0.4 0.4 0.3 0.7 0.5 0.7 0.6 0.5 0.4 0.4 0.5 0.4 0.5 0.4 0.6 0.6 0.7 0.6
Table 8.3: Errors (Hz) relative to the mean that are used for performing inversions. Table layout presented in the same form as table 8.2. Values are
rounded off to nearest tenth in the table.
107
(n,l,|m|) -4 -6 -8 -10 -12 -14 -16 -18 20 -22 -24 -26 -28 -30 -32 -34 -36 -38 -40
(0,1,1) 111 56 60 62 38 64 57 44 44 29 49 20 17 23 18 22 44 22 18
(0,2,1) 33 25 27 37 23 23 27 19 18 17 15 16 23 16 12 10 5 5 7
(0,3,1) 32 35 22 22 26 17 18 12 14 14 11 11 15 10 9 8 7 8 8
(0,3,2) 98 81 52 29 27 32 25 20 16 16 12 12 11 16 11 10 9 10 5
(0,4,1) 28 21 21 19 23 21 16 15 13 16 17 15 13 13 11 12 13 12 12
(0,4,3) 202 111 56 45 25 21 18 19 17 12 12 12 12 10 9 9 9 7 6
(1,1,1) 35 28 24 26 12 8 8 11 8 11 8 7 7 6 6 7 7 7 7
(0,5,1) 115 67 31 32 25 19 20 17 14 16 18 18 18 21 18 15 12 14 13
(0,6,1) 61 39 30 36 33 27 27 33 20 21 21 23 22 20 23 18 22 20 17
(1,3,2) 26 17 14 13 11 7 7 7 6 6 5 7 6 5 7 6 6 6 5
(1,4,1) 83 16 16 15 16 17 14 15 16 17 15 19 19 16 20 15 13 14 12
(2,2,1) 31 23 16 16 12 21 12 8 9 11 12 9 7 5 7 5 6 9 8
(1,5,1) 34 36 40 39 30 32 18 18 15 18 18 13 10 15 18 12 14 13 7
(2,3,1) 21 22 21 17 22 15 13 16 12 11 10 12 10 11 12 9 15 9 13
(2,3,2) 26 23 15 10 19 11 13 11 8 6 6 7 5 6 5 7 7 7 6
Table 8.4: Fraction of measurement that the error radius represents, as a percent. Often the uncertainty is small relative to the absolute shift that
includes the shift induced by the rotating frame, but much larger relative to the residual shift. Table layout same as in table 8.2.
108
8.4 Inversions for Rotation
8.4.1 Constrained Optimization Inversions
The ability to predict the acoustic data with the previously described basis,cylindrically
radial cubic B-Splines and Neumann modes, is shown in figure 8.6, where average misfit is the
average misfit is defined as the difference between the prediction form the inversion and observed
value, denominated in the error radius.The error radius reflects both possible variability due to
flow, as well as decreased confidence associated with measurements described earlier throughout
this chapter.
Figure 8.6: RMS misfit of the ratio of the discrepancy between predicted frequencies from inversions
and the measured values of frequencies to the error weights. The rotation rate of the outer sphere
was 6 Hz, and the inner sphere’s rotation rates were -4,-6,...,-40 Hz and shown in the horizontal
axis. For the purposes of assessing overfitting vs. underfitting the data, the vertical axis values
should be scaled by a factor of between 2 and 4 for converting to uncertainty on the mean. Any
choice of a particular value for scaling would be somewhat arbitrary.
109
One particularly useful idea, is that often velocimetry is a means to an end. It is not actually
the velocity at any one particular point that is of critical importance, it is that the velocity field
can be used to compute quantities that inform one about related phenomena. In that regard,
many issues of finer scale features in inversions do not contribute substantially to variation in
those quantities. So, to first motivate discussion of the reconstructed profiles for different Rossby
numbers, first I will discuss many functionals that can be defined on the recovered profiles, since
they are highly informative and partially motivate interpretations of the inversions.
Functionals defined on the inversions were most revealing, such as total azimuthal angular
momentum, average azimuthal velocity, average angular velocity, and for ωeffect. Finer details of
the reconstructions can often be difficult to interpret the veracity of, especially for high Rossby
numbers, but they typically occur in regions that do not produce pronounced affects on the global
characteristic features evaluated by functionals defined on the profiles.
Consider the notation
%X6 := 100 ·
∫
V
XdV∫
V
Xsbr6 dV
(8.3)
where Xsbr6 means the quantity evaluated when the fludi is roating uniformly at 6 Hz with the
outer shell, in “solid body rotation”.
The inferred azimuthal angular momentum decreases with the square of the Rossby number,
as the fluid is slowed down relative to the outer shell by the inner sphere, eventually rotating in the
opposite direction. The fit is depicted in figure ?? along with the data and variation in the values
found by using an ensemble of 500 inversions for each Rossby number. Each inversion is performed
by choosing values in the uncertainty radius of the measured data, and then inverting with uniform
and small error radii to create the weights in the optimization data misfit functions.The ensem-
ble produced can yield a posterior estimate of a functional Φ(m) defined on the model space, by
computing statistics on {Φ(m) : m = Inverse(dj), j = 1, ..., N}. This is necessary since iterative
methods with penalty functions cannot have their posterior error distributions formed by convo-
lution of distribution functions, as is the case for a linear inverse operators. Using the ensemble
approach, variations due to data uncertainty were thus able to be incorporate to weighting the
110
data.
To determine how large of an ensemble was required, several tests were performed using 5000
inversions sampling the error distributions. See figure 8.7 for an example. Though 500 samples
did not produce as smooth a distribution, the mean and standard deviation converged sufficiently
by 500 samples, and to the same values as the larger sample size, that more samples were not
warranted given the accuracy being employed. Profiles using 30,000 or more were also visually
inspected in earlier work and found to behave similarly.
Figure 8.7: Histogram with Gaussian fit for an ensemble of 5000 inversions for Ro=-4.67, for the
angular momentum quantity described in figure 8.8. The distribution was sufficiently Gaussian
and similar in standard deviation and mean for smaller ensembles that it allowed for sufficient
accuracy in the weighting of the data when performing fits.
111
For net azimuthal angular expressed by %L6 the relationship is
%L6 ≈ 100− pi
5
Ro2 (8.4)
Figure 8.8: The azimuthal angular momentum as a percent of the uniform rotation case at 6Hz,
depicted by triangles. The variation with respect to sampling the data within the uncertainty
interval around presented as error bars, and the curve depicting a purely qudratic relationship
with Rossby number after offsetting by 100% at uniform rotation, i.e. Ro=0. The coefficient in
the quadratic term is ≈ −0.62
.
112
The most interesting aspect of the recovered parameter trends in this section, such as that
shown in figure 8.8, is how well they can be fit by simple functions such as a pure quadratic, a
square root, or a line. In most cases, they intercept intuitive values of the parameter at Ro=0.
Also, in the majority of cases, the simple functional forms fit well enough that the values can
used to invert for the Rossby number if the Rossby number is unknown. This is really significant,
since even if the predicted flows turned out to be innacurate, or the functionals were nto gauging
the quantity they were designed for, the simple relationship between acoustics and the Rossby
number contained in these functions gives the ability to infer the Rossby number directly from
the acoustics. Preliminary analysis has revealed that even smaller basis function sets with less
resolving power and fewer modes in the data set, also fit these profiles. Though those less accurate
inversions imply values of the functionals that scatter around the fitting curves more, the fitting
parameters are still inferred to within a few percent. However, the analysis has also revealed that
the primary contribution to many of these profiles comes from outside the tangent cylinder. That
is, the relationships may not be representative of flows inside the tangent cylinder.
The average velocity can be expressed similarly:
%V6 = 100 + 1.77 ·Ro− 0.54 ·Ro2 (8.5)
The fit is depicted in figure 8.9. A line could also be used for a portion of the data, but it wouldn’t
predict the value of 100% at Ro=0, nor fit the data as tightly over a broad region as in figure 8.9.
The quadratic term in Ro was essential for fitting the profile globally and assuming the value of
100% at Ro=0.
113
Figure 8.9: The fit of %V6 is a second order polynomial in Ro with an intercept of 100% at Ro=0.
.
The average rotation rate as a percent of the rotation rate associated with uniform rotation,
is somewhat interesting as well. It yielded the following relation quite tightly fit by the equation
%Ω6 = 105.4 + 7.66 ·Ro (8.6)
and the fit is depicted in figure 8.10.
114
Figure 8.10: The fit of %Ω6 is linear but does not go through the 100% at Ro=0.
.
These fits are tight as can be seen in figure 8.8, figure 8.9, and figure 8.10. They are
perhaps far more accurate than than any specific location in the cavity reconstructed from an
inversion. Especially since one of the most difficult places to resolve the flow well from inversions
is inside the tangent cylinder, but a volume integral is biased toward values of r that are larger
than they are inside the tangent cylinder. With the exception of the last two data points, it could
be possible that the majority of the variation about the fitting curves is due to variations inside
the tangent cylinder, or the affects on the flow of the shaft and other non-ideal aspects of the
cavity geometry; it also possible though that since the flows are not purely azimuthal and there
is turbulence, there is conversion occurring that is causing slight deviations. Since these are not
high resolution inversions though, it could also be a case of lack of well resolved features in the
azimuthal component. Decomposing into an inside and outside of the tangent cylinder, the outside
was still linear whereas the inside did not have a simple trend, when evaluating the same quantity
115
in the figure 8.10.
The most interesting functional is perhaps the %Ω6 profile. It does not go through 100%
at Ro=0. Trying to force the fit to do so, drastically reduces the accuracy of it’s prediction of the
data. The deviations from linearity in the %L6 and %V6 profiles already imply that the flow is not
scaling with Rossby number, after removing the offset due to the rotating frame of reference.
For net ωeffect increases as the inner sphere spins faster in the opposite direction of the
outer sphere. The non-squared average can also be fit by a simple linear model in Rossby that goes
through Ro=0, but the fit is not as tight for higher magnitude rossby numbers. This is because as
the flow became substantial at higher values of r, the inversions produced what looks like artifacts
of the low resolution basis within the tangent cylinder, and even slightly beyond (see figure 8.14).
This is because to fit the rotational profile accurately where it affects the acoustics more so, does so
at the expense of creating artifacts where the acoustic modes are less sensitive. The side lobes and
non local effects in the smaller values of r regions affects the ω more so since it contains gradients
in r of the rotation profile. In figure 8.11 what can be seen is a somewhat linear trend through 0
at Ro=0, yet by Ro=5 there starts to become not only deviations from the trend, but also larger
variations seen across the ensembles of solutions.
Figure 8.11: Inferred ωeffect, described in chapter 1, as a function of Ro.
116
Figure 8.11 shows the ωeffect as having slope ≈ −0.7 trend with Rossby, once normalized
by the outer rotation rate; the actual value of the trend is likely slightly higher, since the higher
Rossby value data should perhaps not be included in the fit. The normalization could also include
a volume factor, to help nondimensionalize the quantity. The volume is just over 0.1 cubic meters,
and so the slope is between -6 and -7 once nondimensionalized. Since there is presumably zero ω
effect for uniform rotation, the normalization cannot be performed as in previous quantities.
Other complications are that if the shaft is actually slowing down the rotation rate near it,
it would cause gradients that actually decreased the net ω effect. For that reason and others, of
all the quantities examined, the ω effect is probably the least accurate, since it is corrupted more
so by the tangent cylinder interior’s fluid rotation and the gradients of that rotation with respect
to r; the inversions are likely to be less accurate there than average over the volume of the cavity.
By using a steep sigmoid function, approximating an infinite shearing of the fluid at a
particular value of r, parametrized so that it was negative near the shaft and zero on the other
side of the shear zone at higher values of r, constrained optimization over A and rc was performed
to fit the most likely location for the shear zone.
117
Figure 8.12: Estimating the location of high shear. The horizontal axis represents the location of
the vertical line separating the two piecewise constant basis functions, as inferred by minimizing
the misfit to the acoustic data; described in text. Though a line is also possible to fit the data
nearly as well, the square root function looked like a better fit and also goes through zero, which
is consistent with solid body rotation of the fluid. The vertical axis has been non-dimensionalized
by dividing by the radius of the outer sphere.
Figure 8.13: Data misfit for the shear zone parametrization by a sigmoid, using similar layout as
figure 8.6. This corresponds to the data plotted in figure 8.13.
118
The profile suggested by figure 8.13 was chosen to be a square root function of the magnitude
of Ro. The fit for a line is also convincing, but not only does the square root fit better, but it also
goes through 0 at Ro=0. Perhaps it should go through the axle radius or shaft radius at Ro=0
but 0 was an approximately valid location for the origin, where the fluid is rotating uniformly and
has no shear.
This is also in agreement with the previous figures in this section that used the higher
resolution basis, since it is showing the affect on the flow profile in the rotating frame is not a
simple scaling with the inner sphere rotation rate or Rossby number. Unlike the case of the inner
sphere spinning only, the Coriolis force seems to be preventing the shear zone from diffusing out
and the angular momentum from mixing as much.
Other features present, include the rotation of fluid following the shaft somewhat near it,
rather than the inner sphere. This could be to some extent due to the lack of resolution near the
shaft, which was investigated with synthetics and found that there is occasionally drops in the
rotation rate near the shaft not due to actual reduction in the input profile. On occasion the drop
is more considerable though, and it is possible that the large surface area of the shaft spinning
with the outer sphere is reducing the rotation rate as well.
119
Figure 8.14: Reconstructed inversions at each rotation rate. The rotation rate of the inner sphere
has been placed on top of each profile, but all other labelling has been omitted to reduce clutter.
The outer sphere is spinning at 6 Hz.
Broadly speaking, what the inversions reveal is that as the Rossby number increases as a
shear zone moves out from the tangent cylinder towards the outer sphere’s equator.
In figure 8.15 there appears to be a jet instability forming on the inner sphere’s equator.
The rotation rate appears to be less negative in the regions of the column immediately above and
below it than at higher |z| values. It was thought that perhaps this was an affect of the sidelobes
of the feature forming at more distant portions of the cavity, but the variation in the column was
also recovered, and more pronounced using 2D cubic splines. It seems that the undulations in the
profile at higher values of r may actually be the artifact, related to the jet feature.
120
Figure 8.15: Reconstructed flow profile for Ro=-4.67, i.e. Fi = −22 Hz and Fo = 6 Hz. There
appears to be a jet forming on the equator of the inner sphere.
121
Fr
eq
ue
nc
y 
Sh
ift
 (H
z)
Figure 8.16: Data prediction from reconstructed flow profile for Ro = -4.67, i.e. Fi = −22 Hz and
Fo = 6 Hz. The data is plotted in red with the error weights used in the inversion.
122
Tools of assessment for these types of inversion included assessing variability across ensem-
bles, as done for the functions. In figure 8.17 the relative variation
The jet is likely always present, but perhaps less substantial at other rotation rates. It is
difficult to judge the actual size with a low resolution basis. One possible reason for it’s emergence,
is that as can be seen, a large portion of the fluid is coming to be approximately at rest in the
lab frame. The jet is likely much faster in the negative direction, yet the defocussing smears that
velocity feature over a larger volume of fluid in the Reconstruction.
Figure 8.17: Relative variability profile for Ro=-4.67, i.e. Fi = −22 Hz and Fo = 6 Hz. The units
of Hz is consistent with the mean flow error radius. This is only variation across an ensemble
though, not necessarily reflective of the uncertainties related to how well the basis can match the
flow in very localized regions.
123
Figure 8.18 shows the profile prior to the formation of that feature, when the outer was
spinning only -20 Hz. The data misfit is nearly identical, except mode index 5 was slightly better
fit than for Fi = −20 Hz, but investigation found that not be what is causing the jet like appearance
in figure 8.15.
Figure 8.18: Reconstructed flow profile for Ro=-4.33, i.e. Fi = −20 Hz and Fo = 6 Hz. Any jet is
less pronounced in this inversion than in the one for Fi = −22 Hz.
124
Beyond these regimes, i.e. beyond Ro=5, the profiles start developing vertical variations
that are difficult to explain, at least at lower values of r. It could be related to fluid coming to
rest in the lab frame, the shaft, other factors left out of the formulation, or simply that fitting the
profiles well where the acoustic modes are sensitive trades-off errors at smaller r values in a very
disproportionate way.
The last two rotation rates, Fi = −38 Hz and Fi = −40 Hz seem to be outliers in many
analyses such as total azimuthal angular momentum. Several possibilities seem likely. First of all,
there are electronics at 45 degrees θ and 135 degrees θ colatitude at several azimuths, and the fluid
appears to be following the outer shell; the electronics seem too small in volume for a column to
cause that though. It seems to always be present possibly, yet it is not until the fluid begins sub
rotating the outer sphere beyond the column that it is noticeable.
Another strong possibility, is that the bound of not super-rotating the outer has a different
affect on the inversion once the most of the fluid inside is no longer following the outer shell.
Also,, using fewer basis functions so that a traditional weighted least squares inversion could be
performed, recovered a rotation rate there that was super-rotating the outer shell. It could actually
be the constraint causing the issue. It could be possible that simply the specific angular momentum
is non decreasing with r and not velocity necessarily, similar to a conical pendulum speeding up
as the level length is decreased; analogous to fluid falling to a shorter r radius.
8.4.2 Consistency with Alternate Basis Functions
To gauge that features observed in the inversions are not artifacts of the basis chosen,
alternate basis functions were used, as discussed in the previous section. For instance, the jet-like
feature that appears when the inner sphere is rotating at -22 Hz is also recovered by 2D cubic
splines, and also by a coarse tiling by Voronoi cells. Since such basis functions are less prone to
nonlocal effects that does not stem from the data itself, there is confidence that these features are
not artifacts of the parametrization. The concentration of negative rotation in the top and bottom
hemispheres also appears in those alternate parametrizations.
125
Fitting the Voronoi cells iteratively using Quasi-newton methods, and monitoring the inver-
sion, was informative. The jet appears early on, even by the time that the misfit RMS is roughly
1 error weight. That indicates that it is a feature robust across the variations of flow occurring.
By roughly an RMS misfit of 0.6 the stationary region is apparent, and only by the time the misfit
reaches 0.5 or 0.4 does the stationary region begin to become more varying and lobes present in
the hemispheres. Likewise, the rotation rate matching the shaft near the shaft also only occurs as
the misfit decreases similarly.
Polynomials in (r,z) also recover similar features as the basis used in the last section, but
also produce many artifacts analogous to Runge phenomena, and also oscillate wildly where there
is poor coverage of the domain with the sensitivity kernels of the acoustic modes.
A very minimal basis was used, approximating the flow as being vertically invariant and
only using 4 splines from the primary basis used in the previous section. Leaving out the 2 splines
concentrated at the highest values of r. So, the basis tapers to zero at ro. The profiles are shown
in figure 8.20 and the fitting analysis is shown in figure 8.19. The poor fit of the data for the 2
highest Rossby numbers is due to the constraints possibly. Weighted least squares formulations
with smaller basis functions that made no requirements on flow rotation bounds, fit those profiles
well but using values of rotation that super rotated the outer sphere. Another large misfit occurs
at Fi = −22 Hz, which is to be expected since the 2D basis revealed a jet like feature that cannot
be expressed in a vertically invariant basis well. The large misfit at Fi = −6 could be due the
fact that the splittings are small there; such is the case for Fi = −4 Hz as well, but it is perhaps
compensated by the error bars being much larger as well, since there was at times issues of polarity
of the splitting being difficult to determine.
126
Figure 8.19: Reconstructed Profiles using only 4 cubic splines, as discussed in text. Recovered
profiles look similar to how the profiles in figure 8.14 would look if they were averaged vertically
and homogenized in the z variable based on that average. Figure similar labelling convention as
figure 8.14.
127
That figure 8.19 suggests that the vertically invariant basis produces inversions similar to
vertically averaging 2D reconstructions is reassuring. Often in inversions, projection operators and
averaging operators do not commute well wiht the resolution matrix, meaning that basis that is
a subset of larger basis that can express the flow well, will not necessarily recover the projection
of the flow onto that subset of the basis vectors. Again, the absence of relative rotation along the
shaft that begins to appear at Fi = −12 Hz, may be tempting to attribute to the fact that the
shaft spins with the outer sphere, but resolution tests have shown it is also at least spatially due
to how defocussing of shear zones trades off with the error in flow along the shaft.
Figure 8.20: Data misfit analysis for basis consisting of 4 cubic splines, discussed in text. Figure
similar to figure 8.6.
128
Figure 8.20 suggests, is that the shear location selected by the sigmoid basis, is not neces-
sarily the region of highest shear, but the region of shear nearest the portion of the cavity following
the outer sphere.
8.5 Alternate Methods of Inversions
8.5.1 Minimization of Spread Functions
Using 192 2D cubic splines in (r,z) the minimization of Dirichlet spread function approach
was performed. Splines, though overlapping their neighbors, were found to produce less artifacts
than using disjoint tilings. Curvature damping in 2 variables, in this case r and z, was used to model
the problem, and the Dirichlet spread function was minimized to choose the damping coefficient
parameters.
129
Hz
Rotation Prole Data Fit Analysis
Dirichlet Spread Function Minimization
 Inner Rotation Rate of - 22 Hz
Figure 8.21: Result of minimizing Dirichlet spread function for 2D cubic spline basis. It has many
realistic features, but also has some overshoot in the sense that it super-rotates the outer shell in
a way that does not seem plausible.
130
Figure 8.21 has many plausible features and likely is capturing true aspect of the flow. The
overshoot artifacts are unfortunate, since they do not appear of a form that physically plausi-
ble. Their location looks similar to early iterations of iterative methods that used Quasi-newton
methods, suggesting they may a mutual region of sensitivity for two or several modes. It was not
investigated further though, since these results proved difficult replicate with different weighting
functions, bounds set on the damping parameters, and weighting of the spread function. The
trade-offs were also apparent. For instance, at an inner rotation rate of 20 Hz it did show there
was less flow at the inner sphere equator, which is consistent with the other inversions, but also
showed elevated rotation on the top and bottom of the inner sphere, yet less overshoot in the
regions super-rotating the outer shell. It would seem that the variation in the figures could be as
easily attributed to how much it is trading off with other parts of the cavity. It would be hard to
draw anything conclusive from the inversions. Even though a jet is not only possible but likely, if
there are nonphysical artifacts in the inversion, it is not wise to selectively choose which features
to lend credence to and which to ignore without good reason.
All issues of parameter choices aside, even profiles such as in figure 8.21 did a remarkably
good job at recovering that the flow is closely following the outer sphere, and the spread mini-
mization formulation did not even impose an upper bound on the rotation rate of the fluid in the
inversion.
Combining the constraint that the solution must be non-positive did not improve the plausi-
bility, or lift the subjectivity; either by imposing condition for the inverse operator when choosing
damping parameters, or by choosing damping parameters and then imposing the restriction by
solving the inverse problem constrained.
The Backus and Gilbert approach had even more issues of subjectivity. Unlike the Dirichlet
spread function, the weighting of the spread function was much more important. Constraining
the rows to sum to 1 of the resolution matrix was not ever close to being obtained, and weighting
differently whether by volume or angular momentum did not seem to make much difference. Also,
using disjoint tiles instead of splines, in the (R, θ) plane or the (r, z) plane always had issues at
131
the boundary, since neither could match the boundary everywhere. The splines could have their
second derivatives evaluated by interpolating their gradients at the triangle centers back to the
nodes and then differentiating at the triangle centers again. Dealing with the fractional areas at
the border of the domain did not seem to have any issues with the splines, yet for tiles it would
often produce zero as the solution near the boundary even when weighted in many different ways
to make up for the fact that the area or volume contribution was different. The same issues were
found for both curvature damping in r and z, as for R and θ, using tilings. The use of these methods
was a constant search for parameter combinations that yielded plausible results, and so at best
seem like confirmation bias. At least theoretically, these methods are very appealing, and would
no doubt perform quite well with sufficient data. They also allow for the rapid propagation of
uncertainties via convolution of the distribution functions associated with the uncertainties on the
measurements; though given the types of errors, it is not completely clear if a standard probability
interpretation of the posterior distribution for the rotational profile is fully logical.
Lastly, judging a minimizer by a convergence of the optimization method, here again using
interior point method, often produced multiple minimizers depending on the initial conditions, i.e.
the damping parameters the iterates started at, though typically only one solution was plausible,
or not at the boundary of the constraints. The other would have many localized lobes or be almost
completely flat and not fit the data at all when used for inversions.
8.5.2 Interpolating on a Fixed Grid
The domain can be decomposed via constrained Delaunay triangulation, and the flow can
be reconstructed by interpolating on the nodal values. Instead of using a forward operator, the
kernels are used directly to predict the acoustic effects, and values at the nodes can be determined
by traditional optimization methods, such as Broyden-Fletcher-Shannon-Goldfarb algorithm. Dif-
ferent methods of interpolation are possible. The simplest is the nearest neighbor, which tessellates
the meridional cross section into Voronoi cells each cell having uniform rotation value. Another
simple method is simply linearly interpolates between the nodal values. Finally, natural neighbor
132
interpolation produces a smooth interpolant except at nodes.
One downside to the slowness of these methods, is it makes propagating uncertainty asso-
ciated with the data more cumbersome computationally expensive. Furthermore, these inversions
were not checked with a plethora of initial conditions to ensure the solution is not a local mini-
mizer of the objective function. However, since alternate means have been presented already that
do perform these checks of the answer, when these means produce inversions that are consistent
with those alternate means, it increases reliability. When the inversions differ in some feature
from method to method, what it can be quite telling. If the data misfit they predict are equally
plausible, it says that a feature is not in any way required by the data and it’s reliability is low.
That may not even be apparent by propagating the data uncertainty, since a single parametrization
could have biases not stemming form the data, but from the span of the basis functions possible.
When the inversions disagree on some feature, yet have significant differences in how they predict
the data accurately, analysis can be very informative. If a single mode’s splitting is not predicted
well by one, it could allow reassessment of the nature of that data uncertainty or upon examination
of the sensitivity kernel determine that the feature in the inversions in question is related to that
specific mode splitting that the flows predict differently.
Another aspect of variability is the choice of points to use to define the grid. When us-
ing many dozen of points somewhat uniformly sample the domain this type of variation is less
pronounced, but small amounts of data discourage using too fine a mesh for the node locations.
8.6 Mean Flow Approximation Issues
Since the flow is varying about the mean flow, even when the mean flow is steady, it means
that there are accelerations occurring. The Euler forces associated with Ω˙ is ignored, approximating
them as zero in terms of their effect on the low frequency modes used. It is expected that those
effects are quite small, and the linear manner in which they add to the force implies they will
average out, since the average value must be zero; or at least very small if the flow is close to being
steady.
133
The mean flow is assumed to be steady, which may not always be warranted, especially near
a bifurcation value for inner rotation rate. Hysteresis can occur and other forms of oscillations
between states, as well as instabilities that obscure what the mean flow represents. Such phenomena
are easy to conflate with peak clusters created by the noise envelopes though. The multiple peaks
clusters occur around the m = 0 modes as well, which strongly suggests noise envelope effects,
though it is also possible to some extent the phenomena is due to other forms of coupling such as
described in Roth and Stix [1999].
8.7 Alternate Basis Discussion
Polynomials produced nearly identical results, using r and z up to order 5 or greater, but
keeping only even orders of z, except would produce far more artifacts such as oscillating at the
boundary or regions where sensitivity was low; greater maximum order producing more artifacts.
The same features were produced using polynomials by synthetic flows that did not have those
oscillations in the flow pattern, which supported the contention that is was just an artifact. Polyno-
mials tend to produce Runge phenomenon that behaves far worse in conjunction with constraints
than the basis I used herein. That is, polynomials of low order will prefer to oscillate around the
solution, and a large portion of the solution likes to be right on the boundary of the constraint.
Relaxing the constraints can remove the problem, but also allows a different type of artifact. The
Neumann Modes that have zonal harmonic based factors do not seem to suffer the analogous Gibbs
like effect to as marked a degree, though it is because they are supplemented by cubic splines in r.
There is one inversion where a pronounced feature forms around where the jet instability is located,
which does appear to produce what is likely sidelobe artifacts due to the modal basis components,
since it is difficult for the cubic splines to adequately support that high vertical gradient at the
equator of the inner sphere.
2D Cubic splines, coupled with their mirror image about z=0 to produce a basis functions
even about z=0, also produced markedly similar effects. That was especially important in the case
of the jet like instability that appeared in the rotation rate of Fi = −22 Hz with Fo = 6 Hz. Since
134
only even l zonal harmonics were used, all the modal components of the flow basis had lobes in
that region. Using least squares projections with synthetics showed that artifacts were never that
pronounced as in that inversion, and the fact that the splines recovered it as well showed that it
was not an artifact of the basis chosen; it could only be a basis given by the acoustic modes used for
data. The cubic splines do not have direct nonlocal effects except at neighboring regions, and so it
also lent weight to the idea that there was some artificial sidelobes created by the low order basis
having to represent that feature, since they are absent in the spline version. The spline basis did not
have the vertically invariant basis components included in their basis, and tended to have slightly
more variance in the z direction. While it is possible there is more variation in z direction, since
the forward operator was horribly conditioned if the inversion was tried int he traditional manner
rather than using constrained optimization, even with as many basis functions as data points, and
also the error analysis revealed the regions of higher variation in the z direction were also regions of
higher uncertainty, the basis was not preferable. As a diagnostic device for laboratory experiments
and assessing if flow changes are occurring concomitantly with magnetic phenomena, I recommend
using multiple bases in order to become more fully aware of these aspects. The reason for this
commentary, is to convey that a great deal of effort was put into verifying that any feature in the
flows was relatively robust to the basis chosen; any bias was due to the limited number of modes
available as data. Gauging the effects of coverage issues, synthetics were used to understand what
kinds of defocusing were possible. The basis used here was found to perform best and produce the
least amount of artifacts across the different differential rotation rates. This is because it has less
resolving power perhaps, but the data has little resolving power relative data sets such as those
used in helioseismology.
135
Chapter 9: Concluding Discussion
9.1 Conclusions
9.1.1 Application of the Finite Element Methods
Flow and Acoustics cannot easily be jointly parametrized with a basis that would allow
an analytic construction of the forward operator matrix. Even when the modes and frequencies
can be well approximated by the spherical annulus approximation, discretization typically must
occur in order to calculate the forward operator. An earlier incarnation of my software used
the acoustic modes of the actual geometry to parametrize the flow, but found to be suboptimal
since the flows of interested tend to exhibit more cylindrical symmetry and variation with r than
spherical symmetry or variation easily characterized by varations in θ. A numeric approach must
be employed then to perform the inversions, and little is gained by using analytic expressions
for the acoustic modes in these contexts. Many methods in helioseismology benefit from those
expressions since they allow for asymptotic formulations and more advanced forms of inversions
for high frequency modes and amplitude analysis. For these types of laboratory experiments, the
finite element method is optimal.
9.1.2 Bulk Fluid Rotation Recovery
The greatest success of the inversions, has perhaps been ignored due to the context of the
study: across many different methods and parametrization, the recovery of the bulk fluid rotation
was quite robust. This was perhaps less impressive in these experiments, since the fact that the
majority of the fluid would follow the outer shell rotation rate was something that could have been
136
surmised without the inversions.
Optimizing the resolution matrix often produced super rotation artifacts or implausible
rotation rates inside the tangent cylinder. However, these undesirable features were only large
relative to the flow in the rotating reference frame. They were quite small errors relative to the
bulk of the fluid rotating at approximately 6 Hz, and the methods all recovered that the fluid was
largely following the outer sphere reasonably well throughout a substantial portion of the cavity.
Even the constrained inversions using the primary method used in chapter 8, when the constraints
were greatly relaxed, also recovered the overall following of the outer sphere.
The functional that produced simple formulas for their values as a function of Rossby num-
ber, allow the Rossby number of the flow to inferred quite accurately, if that number was not
known. It is of little use in an experiment where the fluid has been brought to a relatively steady
state and the Rossby number given by the spheres’ rotation rate characterizes the flow, but in
many other contexts that is a very useful ability to have. It provide a direct way to measure bulk
fluid properties form the acoustics, and will work whether the functionals accurately describe the
variables they are formulated as or not.
9.2 Looking Forward
9.2.1 Liquid Sodium Experiments
The three meter setup at University of Maryland has pressure probe ports at equally spaced
azimuths, and should in theory facilitate identification of modes by azimuthal order ms by tra-
ditional array techniques. There have been multiple difficulties identifying the modes in pressure
probe data in the sodium experiments. This engineering aspect of the acoustic velocimetry system
will hopefully be resolved in the not too distant future.
9.2.2 Poloidal Flow
If poloidal flow were inducing more shifts than the small ones discussed in chapter 8, it would
have corrupted the inference that temperature was causing a uniform shift. The ms = 0 modes
137
have different sensitivities to poloidal flow, so if the mode frequencies drifted as if by a uniform
temperature shift, it seems to suggest that the poloidal flow effects were only causing small drifts,
possibly the cause of the clusters of peaks around the mode spectra discussed in chapter 8.
It seems unlikely that temperature and poloidal flow could trade-off in concert so as to create
this illusion of the overall drifts being explainable as a uniform temperature shift. More accurate
temperature readings could possibly support this better. All evidence indicates an approximately
uniform soundspeed to the degree of accuracy needed in this study though.
9.3 Recommendations for Future Improvements
9.3.1 Inversion Formulation
This work developed many of the methods under the assumption of a larger data set. That
is still the hope for the data sets once liquid sodium acoustics are functioning properly. However,
if only such a meager data set is to be used, there is perhaps greater utility in actually computing
the predictions of acoustic splittings without using perturbation theory.
As a test of idea, I coded the acoustic wave equation in the material displacement domain
for an infinitely tall Taylor-Couette flow using a FDM. It can also be solved using quadratic
eigenvalue problems, but it was simple enough to use an iterative method. I begin by assuming the
stationary cavity frequencies, insert them into the convective component of the wave equation in the
presence of rotation, solve for the perturbed frequencies, then update the equation and repeat until
convergence is attained. This method of successive approximations allows the shifted frequencies
to found without using perturbation theory; though it requires them to be found one at a time
since the updated forcing terms are only valid for a single acoustic frequency. Furthermore, the
flow can be parametrized in ways other than coefficients of basis functions. The model parameters
can be exponents or other parameters, and the inversion can be carried by Quasi-Newton methods
or other techniques in nonlinear optimization. It is somewhat more computationally costly though,
especially if error propagation is sophisticated.
138
9.3.2 Poloidal Flow
Using a very targeted chirp with high amplitude, it should be possible to assume frequency
fluctuations in the acoustics come primarily from flow or thermal effects. If the thermal effects can
be inferred to be uniform from, then the movements of center frequencies can be attributed to Non-
azimuthal flow. In the case of liquid sodium, there could be corresponding magnetic phenomena
associated with poloidal type flows. The kernels associated with them can be easily computed
from the convective operator applied to the acoustic perturbation.
9.3.3 Modeling Geometry
The missing pieces of the shaft seemed to be of less importance, and actually a source of
nuisance. They affected the frequencies predicted negligibly, and instead of preventing corruption
to the flow, possibly increased it.
There was also an apparent artifact that appeared near the corner of the shaft inside the
missing piece in the upper hemisphere in a few of the sensitivity kernels. This could be possibly
due to Neumann condition being somewhat poorly defined at a corner, since at a corner the normal
derivative is singular technically. This artifact did not appear in the modes, and seemed to be
due to computing gradients on the triangular mesh. Different methods were used to compute the
gradients to remove the artifact to no avail, and it was present when computing the kernels from
modes I computed using the Matlab PDE toolbox as well. It only occurred over a very small
portion of area, and was not significantly large in amplitude, so after concluding that zeroing out
the sensitivity there and also smoothing the mode or kernel did not produce visible differences in
inversions, it was ignored. It is likely due to the fact that despite many modes being nearly zero
there, the high stiffness encountered at the corner caused some lack of smoothness. Ultimately
though, the region’s geometry was highly idealized, and the model’s deviation from reality was likely
far greater than any numerical aspects such as this, but if I had to do it again I would smooth out
all corners, for cosmetic reasons and to prevent having to test the effects of the possible artifact.
The issue did not occur in my models of the 30-cm setup or the 3-m setup since they had only an
139
axle though; it only seemed to occur on corners that protruded into the cavity.
9.3.4 Meshing Geometry
A great deal of exploration with meshing techniques was performed. The methods set forth
in Persson and Strang [2004] were by far the best approach. The technique is improved mesh
quality and the regularity of the triangles in the mesh outperformed anything I was able to do
easily with other other available software. The method was also easy to understand, so coding it
with the specifications of the particular application did not require much time investment.
Adaptive meshing and nonuniform meshs were also explored, but for solving for modes, the
linear algebra was simpler to assess the numerics of when a very regular mesh was used and the
matrix entries of the FEM operators had relatively equal magnitudes.
9.3.5 Test Functions and Shape Functions
Using a much finer basis during experimentation with minimizing spread functions and with
more complicated constraint equations often involved computing second derivatives. One natural
improvement would be to abandon the piecewise linear formulation so that the second derivatives
could be computed inside each element, and the shape functions coudl be used as the bass func-
tions. Gradients can be computed naturally for piecewise linear elements, but second derivative
calculations often resembled finite difference computed derivatives as a result of not having cur-
vature in the test functions locally; I used interpolating functions and computed the answer via
difference methods. When spline functions or polynomials were used, since the derivatives could
be computed analytically then evaluated on the mesh, but the freedom to compute the second
derivatives as straightforwardly as the gradients would be helpful when using the shape functions
as basis functions.
140
9.3.6 Inclusion of Magnetic Data in Constraints
Possibly by including magnetic data based inferences the constraints can be improved. The
more information about the flow that can supplement the data the better. It is possible electro-
magnetically to determine the fluid velocity near the outer shell. In the case of liquid sodium, it
also possible to tell how much torque is being required to maintain the rotation of the inner sphere.
All such information can be incorporated to augment limited acoustic data, thereby increasing the
amount of information the acoustic data can convey. It has also been suggested in [Adams, 2016]
that selection rules could be used to suggest the possible flows inside, though not unambiguously.
In which case it would only require simply testing the different possibilities using the forward op-
erator for the acoustics. If the different flow states have sufficiently different acoustic signatures,
the flow could be inferred without carrying out an explicit inversion.
9.3.7 Adaptive Interpolation
Allowing the grid points to move as well, introduces more flexibility into the inversion. I
ran some preliminary tests, and found that the grid points tend to cluster so as to be able to form
regions of high shear. Perhaps the method would be ideal to capturing high shear without having
to use a very fine set of basis functions.
141
Bibliography
Santiago Andre´s Triana, Daniel S Zimmerman, Henri-Claude Nataf, Aure´lien Thorette, Vedran
Lekic, and Daniel P Lathrop. Helioseismology in a bottle: modal acoustic velocimetry. New
Journal of Physics, 16(11):113005, 2014.
Daniel P Lathrop and Cary B Forest. Magnetic dynamos in the lab. Physics Today, 64(7):40–45,
2011.
Emmanuel Dormy, Jean-Pierre Valet, and Vincent Courtillot. Numerical models of the geodynamo
and observational constraints. Geochemistry, geophysics, geosystems, 1(10), 2000.
Gary A Glatzmaier and Paul H Roberts. A three-dimensional convective dynamo solution with
rotating and finitely conducting inner core and mantle. Physics of the Earth and Planetary
Interiors, 91(1):63–75, 1995.
Gary A Glatzmaier, Darcy E Ogden, and Thomas L Clune. Modeling the earth’s dynamo. The
State of the Planet: Frontiers and Challenges in Geophysics, pages 13–24, 2004.
Ulrich R Christensen, Julien Aubert, and Gauthier Hulot. Conditions for earth-like geodynamo
models. Earth and Planetary Science Letters, 296(3):487–496, 2010.
RA Secco. Viscosity of the outer core. Mineral Physics and Crystallography, A Handbook of
Physical Constants, pages 218–226, 1995.
Eugene N Parker. Hydromagnetic dynamo models. The Astrophysical Journal, 122:293, 1955.
JB Taylor. The magneto-hydrodynamics of a rotating fluid and the earth’s dynamo problem. In
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences,
volume 274, pages 274–283. The Royal Society, 1963.
Matthew M Adams, Douglas R Stone, Daniel S Zimmerman, and Daniel P Lathrop. Liquid sodium
models of the earths core. Progress in Earth and Planetary Science, 2(1):1–18, 2015.
Peter Olson. Laboratory experiments on the dynamics of the core. Physics of the Earth and
Planetary Interiors, 187(1):1–18, 2011.
Y Takeda. Development of an ultrasound velocity profile monitor. Nuclear Engineering and Design,
126(2):277–284, 1991.
H-C Nataf, Thierry Alboussiere, Daniel Brito, Philippe Cardin, Nadege Gagniere, Dominique
Jault, and Denys Schmitt. Rapidly rotating spherical couette flow in a dipolar magnetic field:
an experimental study of the mean axisymmetric flow. Physics of the Earth and Planetary
Interiors, 170(1):60–72, 2008.
Daniel R Sisan, Nicola´s Mujica, W Andrew Tillotson, Yi-Min Huang, William Dorland, Adil B
Hassam, Thomas M Antonsen, and Daniel P Lathrop. Experimental observation and character-
ization of the magnetorotational instability. Physical Review Letters, 93(11):114502, 2004.
142
Michael J Thompson, Jørgen Christensen-Dalsgaard, Mark S Miesch, and Juri Toomre. The
internal rotation of the sun. Annual Review of Astronomy and Astrophysics, 41(1):599–643,
2003.
Matthew Adams. Magnetic and acoustic investigations of turbulent spherical couette flows. Dis-
sertation, 2016.
Per-Olof Persson and Gilbert Strang. A simple mesh generator in matlab. SIAM review, 46(2):
329–345, 2004.
Yousef Saad. Variations on arnoldi’s method for computing eigenelements of large unsymmetric
matrices. Linear algebra and its applications, 34:269–295, 1980.
Richard B Lehoucq, Danny C Sorensen, and Chao Yang. ARPACK users’ guide: solution of
large-scale eigenvalue problems with implicitly restarted Arnoldi methods, volume 6. Siam, 1998.
Richard H Byrd, Jean Charles Gilbert, and Jorge Nocedal. A trust region method based on interior
point techniques for nonlinear programming. Mathematical Programming, 89(1):149–185, 2000.
Philip E Gill, Walter Murray, Margaret H Wright, et al. Numerical linear algebra and optimization,
volume 1. Addison-Wesley Redwood City, CA, 1991.
M Roth and M Stix. Coupling of solar p modes: quasi-degenerate perturbation theory. Astronomy
and Astrophysics, 351:1133–1138, 1999.
143