ID 
Date 
Author 
Type 
Category 
Subject 
121

Tue Dec 5 10:50:54 2017 
aaron  Mechanics  PonderSqueeze  FEA on optimus  I had some trouble running this on optimus.
Optimus has COMSOL 5.1 installed, but I made these files in 5.3. I downloaded the comsol 5.3 dvd.iso file last night, but on install I'm now getting the error "No locks available." I wasn't sure if this is a file permissions issue (sometimes the file has been 'locked' when it is saved and I open it in the GUI), an svn issue (most of the information I find online is from people running into some problem with lockd or their NFS), or a licensing issue (I don't think this is likely, COMSOL on my laptop shows more licenses available).
Has anyone seen this? Do I need to do something else with /cvs after install? 
122

Tue Dec 5 19:50:47 2017 
aaron  Mechanics  PonderSqueeze  FEA on optimus  Gautam advised me against trying to install version 5.3, lest it break version 5.1I had already gone through the install, but looking at the install manual it says it shouldn't affect previous installs except that the default behavior when double clicking .mph files will always choose the latest version of COMSOL. Since we mostly run on the command line we should be fine. That said I haven't tested files with COMSOL 5.1.
Zach pointed me to sandbox1.ligo.caltech.edu at the lunch meeting, as well as these notes from Rana about running comsol on a remote server. I couldn't run the file myself on sandbox1 because I don't yet have a home directory there, but I asked Larry if he could set one up for me so I can use the clusters. Gautam helped me run them on sandbox1 by having me scp the file onto the shared 40m drive, then he moved it in to his home directory on sandbox1 and ran it. A version of the model asking for an analysis at 60Hz and 600Hz ran in 1020 min, and looked good (I was able to scp the output to my laptop and open it in the gui, though I have a script not quite to the point where Matlab can do this for me, it's not ready yet and anyway I wanted to directly see that the output was normal). I modified the study to take 100 steps between 60Hz and 5kHz, and Gautam has now started that run, which if it scales linearly (hopefully better than linear, since it won't have to remesh and etc) will take about 816 hours.
Thank you Gautam for the extended mattermost session helping me run these! 
123

Tue Dec 12 11:50:12 2017 
aaron  Mechanics  PonderSqueeze  FEA on optimus  Simulation results
First run
Gautam ran the COMSOL model on sandbox1 since we were trying to run it before I had a home directory there to run from my login. Since that first run, Larry set me up on sandbox1 so I was able to run a few more times with some tweaked model parameters.
Here are the results from the first run, which uses the nominal aLIGO test mass parameters:
I don't think I was expecting such a regular pattern here. Are the modes really that closely and regularly spaced? This might also be indicating some problem with my meshing, I could imagine a few scenarios. To go from this chart to something that looks like the suspension thermal noise plot, I'll have to scale this by the frequencydependent loss angle and a strain sensitivity TF probably? I need to remind myself for another few minutes on this, maybe chat about it.
Relevant properties of the file that might need to be tweaked:
 F_load
 The applied force may have been resulting in large displacements of the test mass (and thus large curvature in the fibers). Since the model contains geometric nonlinearities (I think it has to, no?), large displacements may not be able to be linearly extrapolated back to the small displacements we expect. Since I am running a frequency domain study, COMSOL expects a load force, which means I am specifying a force rather than a displacement. In future runs, I wanted to reduce the displacement by scaling down the force... I'll discuss my first attempt at this later.
 skin_depth
 Gabriele mentioned that when he models surface losses for wafers he does not explicitly mesh a boundary layer, but rather asks COMSOL for a boundary integral (so just the 2D integral on the outer surface of the wafer), and then presumably makes some assumption about the thickness of this layer to go from J/m to J if need be. In contrast, my suspension thermal noise model defines the 'surface' of the fiber as the 15um outer layer, and I explicitly mesh that surface volume and do a 3D (volume) integral over the surface domain.
 Using a 2D integral would significantly simplify the model and probably run faster and result in fewer bugs when the scale of the geometry changes the necessary fineness of the mesh. However, I am a bit uncertain about whether this is what we want. At least in my understanding, we are interested in asking about what happens when the surface layer depth becomes comparable to the radius of the fiber. Since the 2D integral would always represent an infinitesimal volume, and the geometry of an extended (deep) surface layer would differ significantly from that of an infinitesimal shell, so I am concerned that the strain energy density would vary significantly with the radial postition in the fiber. Nonetheless, even the mesh I am using now has only a few (24) mesh elements along the radial direction, so I'm not sure I'm doing much better even with the additional computation time. Therefore, in future runs I'd like to try Gabriele's suggestion.
 Frequencies
 I was using 50Hz spacing on the frequencies from 60Hz to 5kHz. I haven't been very systematic about this, but I'm getting some convergence issues going to lower frequenciesmaybe by taking the suggestion for eliminating the skin_depth I can make a finer mesh and go to lower frequencies. This frequency sweep is too coarse to find welldefined resonances, but they are suggested in the plot below.
More recent runs
I wanted to automate this loop over mass parameters in a matlab script, so I set up the Livelink handshake so Matlab would send the model to sandbox1 for solving, then MATLAB on my machine would work with the solved model to extract results. I realized later during running that this might not be optimal, since it will require me to keep the connection to the remote COMSOL server during the entire run, which is A. annoying and B. risky because these might take many hours to run and the connection can easily be severed, wrecking the job. I don't really have a workaround for this, so my current plan is to continue logging on to the remote server and running individual COMSOL files. As I'll discuss, this is probably necessary for now anyway since I ran in to some problems running these models, and might want to tweak the models on my local GUI before running them over a large (frequency) parameter space on sandbox1.
aLIGO parameters, but decrease F_load
Based on the plots generated from the first run, I estimated that the TM displacement was comparable to the TM thickness; just as a first pass I figured 'small displacement' could be easily defined as 'displacements smaller than the thickness of the fiber' or maybe 'smaller than the thickness of the surface layer'. If the displacement is directly proportional to the force, this meant I wanted to scale the previous F_load by 1e5, which I did. I wanted to just see what kinds of displacements I would get, so I asked to run the model only at 300Hz. After more than 45 min of running, the server threw an error that it couldn't find my STEP file for the horns, which were in a local directory but not in the remote (server's) directory.
I thought this was an odd error for two reasons. First, I had originally not even had the STEP in my local directory (it only needed to load once and then could be used many times for a given model in the GUI), and I was getting an error within a couple minutes of starting the job that it could not find the file. Adding the file to my local directory seemed to solve this problem, as the model was running for much longer. The error I got after it finally crashed was not the same error as before, but was still an error in loading the file, which makes me a bit confused about where it is actually looking for this file or when exactly the loading process happens. To solve this, I have copied the STEP to the remote directory and will run again with the import pointing to that remote file, which I suspect will solve the problem.
aLIGO with the mass scaled by 0.1, also decrease F_load by 1e5*0.1
Since the first run seemed to mostly work up to the outstanding questions mentioned above, I decided to also run the model in batch mode directly from the sandbox1 command line, just as Gautam had. That is, I ssh on to sandbox1 and run the command
$ /usr/local/comsol53/multiphysics/bin/glnxa64/comsol batch inputfile [inputname].mph outputfile [outputname].mph
For this run, I scaled the mass of the model by 0.1, so we now have a 4kg TM. Most lengths of the model then scale by 0.1^(1/3), but those related to the radius of the fibers (and the bend point of the fibers) scale by 0.1^(1/2) because we want to maintain a constaint stress in the fiber. The main length of the fiber remained the same so the modes would be in about the same place. I document these more thoroughly in the matlab script, which I should upload to the git. I scale the force by 1e5 as above (same back of the envelope calculation), with an additional factor of 0.1 to account for the lower mass.
The model ran in to some convergence issues after a few frequencies. I could solve this by changing the relative tolerance, but I think I will most likely instead pursue Gabriele's suggestion above and try to refine the mesh to improve the convergence. The solver nodes of the COMSOL models are still a bit mysterious to meI don't really know much about when different methods or measures of convergence are appropriate. Probably playing with these could both improve the accuracy and efficiency of the model... but it seems like a hefty undertaking. Might be worht the long term payoff though.
The first few frequencies did converge, but as I was extracting the resuls ancha, the COMSOL license server, went down so I'll have to extract them later.
NOTE TO SELF: WHAT'S WITH THE INLINE IMAGE QUALITY? DON't USE JPG !!!! 
125

Mon Jan 22 21:13:25 2018 
aaron  Mechanics  PonderSqueeze  tips from Shoaib  I talked with Shoaib about some changes I could make to the FEA model to improve convergence and reduce memory usage. Summary:
 use a hex mesh rather than tetrahedral
 Use more structured meshes. In particular, I can make an angled swept the mesh in the tapered portions rather than using a free mesh in these regions, defining the mesh only on one boundary
 Use a nonconformal mesh, so adjacent domains do not need to have matching meshes. This could allow transitions to coarser meshes as the fiber becomes thicker or contacts the horn
 try using curved elements so the curvature isn't just approximated by the number of regular mesh elements
 Might try scaling the model uniformly by some factor (100, 1000), which could avoid machine precision issues
I'll get these implemented this week and see if the computation goes through. 
126

Mon Jan 29 23:02:13 2018 
aaron  Mechanics  PonderSqueeze  tips from Shoaib  I started implemented some of these changes:
 Started the mesh with a boundary free quad mesh on the interface between the upper tapers and the main part of the fiber. I used the following size setting
 Maximum element size is fiber_taper_length, which I felt was a good characteristic maximum because it wouldn't make sense to have radial elements larger than the length of the taper itself. This setting does not limit the mesh (elements are much smaller than the maximum)
 Minimum element size: skin_depth/10, where I use the thickness of the surface layer as the characteristic length, and allow a mesh with elements substantially smaller than this characteristic length.
 Maximum element growth rate I set at 8, which from some trial and error isn't strongly affecting the mesh as long as it isn't brought much lower than this (so the mesh is allowed to get coarse relatively quickly, though other settings limit this)
 Curvature factor is 0.7, and actually does affect the mesh strongly; for now value seems to give a 'pretty' (symmetric to the eye) mesh.
 The resolution of narrow regions is set to 5, where this value again strongly affects the mesh and this value was chosen because it looks 'pretty'
 Note that in contrast to the previous way I had defined the mesh, I now am defining the surface and bulk of the mesh in the same step, which seems to help these two regions match up. Shoaib had mentioned using a nonconformal mesh, which I may end up implementing in a later update.
 Instead of doing a free mesh in each noncylindrical region, I continued with a swept mesh along the length of the fiber (still sweeping in sections, though I'm not sure this is necessary)
 I also slightly increased the number of mesh elements along the length of the fibers (ie, used more elements in each section of the mesh sweeping).
 In this model, I still have a free tet mesh on the horn and test mass, since those seemed a bit more involved to move over to a hex mesh (there is no free hex mesh in COMSOL that I could see right away, and I'm not sure that the horn can be meshed with a swept mesh due to its irregularity. I'll look in to this further, but I had no immediate solution for making a hex mesh in the horn) and anyway they are much coarser than the fibers so I don't expect them to cause much trouble.
These changes were fine with a 10um mesh, but ended up with much too fine a mesh (~80000 boundary mesh elements in the initial mesh between the taper and main fiber) due to meshing the surface and boundary in a single stage in the first step. I separate them out to get a coarser mesh in the bulk, and also make the resolution a bit coarser. 
129

Sun Aug 26 16:42:54 2018 
rana  Mechanics  Analysis  Test Mass Thermal Noise: Consistency Checks  For the Voyager test masses, we have been considering a barrel coating to increase the IR emissivity to increase the radiative cooling power. We also seek to estimate the added Brownian thermal noise that arises from this.
Dmitry Kopstov (from MSU) made a baseline model for this which we have been modifying. The latest is in the CryogenicLIGO git repo in the FEA directory as ETM_NASAbarrel.mph.
Comparison with Analytic Estimates:
Convergence Tests:

131

Thu Feb 14 12:38:51 2019 
Ching Pin  Mechanics   comsol modelling  So I did a simple comsol model of laser heating of a silicon disk, with only radiation, to see the temperature variation at steady state, which could be the limiting factor for high Q at 123 K, due to the thermalelastic effect.
The model just uses a simple 2 inch disc, at 0.028 cm thick, with the flats not incorporated in yet.
I had to search for silicon thermal conductivity and heat capacity at low temperatures, settling with k= 800 W/(m K) and C_p= 300 j/(kg K) from refering to papers. Will check LIGO documents for more accurate versions.
I put an arbitary boundary condition of constant temperature of 123 K on a spot .2 mm in radius, to simulate a beam.
Other arbitary values include 77 K for ambient and a surface emissivity of 0.5.
The laser is off center, because that it where the laser will enter the current setup.
We can see that the power required is .02 W, which seems reasonable.
The model is consistent with the analytic model I made with the laser beam at the centre of the disc. See last two figures.
I'm still trying to get the time dependence to work, as it is just giving me nonsense right now.
Some thoughts: beam radius affects the temperature variation quite significantly, with a fat beam (1 mm radius) having half the temperature variation as a beam of .2 mm radius
I think the halo is just a trick of the eye, but I could be wrong.
Things to do:
Find the time scale of the system, as we want to modulate the laser to adjust the temperature, which will then be run though the mode ringer to measure Q to find the zero crossing
Change the heat source to be an actual laser
Add in the solid mechanics part
Add in the sapphire lens underneath 
Attachment 1: Screenshot_from_20190214_124711.png


Attachment 2: Screenshot_from_20190214_150740.png


Attachment 3: Screenshot_from_20190214_154141.png


Attachment 4: graph.pdf


132

Fri Feb 15 21:05:31 2019 
Ching Pin  Mechanics   comsol modelling 
So I got the time dependence to work, but I'm not sure what went wrong in the first time anyways. I'll trying to get a sense of how long it takes for the temerature to semiequlibrate, and coming to grips with comsol as a whole. There seems to be some inaccuracies when the timing increases, so I'm having to figure out how to increase accuracy without drastically increasing compute time. On the bright side, I switched the model to heating via deposited beam, for a more accurate model.
Attached is comsol's handling of a deposited beam modulated sinusoidally with a frequency of 0.1 Hz and screwing up badly. Y axis is the average surface tempurature across the whole disc.

Attachment 1: Screenshot_from_20190215_214001.png


133

Tue Feb 19 19:52:53 2019 
Ching Pin  Mechanics   comsol modelling  The time step response to heating via laser (22.5 mW) is given in the attached picture, for 2 starting temperatures, 122.5 K and 122.8 K. We see that it takes fairly long to equlibrate, with a time constant of about 500 s, and is consistent across both temps. The y axis is average temperature across the surface of the disc, and the x axis is time. I believe that the heat distribution profile would be very similar with time, simply because of how much faster conduction is compared to radiation 
Attachment 1: Screenshot_from_20190219_195210.png


134

Fri Mar 1 19:33:40 2019 
Ching Pin  Mechanics   comsol modelling  I've changed the heating to be from two heat sources, to better model the situation with a heater and a laser. The heater deposits 22 mW, with the laser deposting .5 mW. The overall temperature distribution is smaller then before, as expected, but doesn't really change much. The heater is simulated with a deposited beam with a gassian beam profile with a standard deviation (s.d.) of 8 mm. The laser to the size has a .3 mm s.d. for contrast. I learned that while the deposited beam power doesn't care for emisitivity, it cares about the area the beam is incident with, so for example, if you increase the s.d. too much, you get less power deposited then what you enter.
I've also got modes to appear using solid mechanics, and I'm trying to see if there is a good way to get the frequency dependence with temperature nicely simulated. Changing the parameters of the simulation does give me my different frequencies, but I trying to find a way to do that over the time evolution of the model. I also got to check if the frequency shift is in line with what we are measuring. 
Attachment 1: Screenshot_from_20190301_193532.png


Attachment 2: Screenshot_from_20190301_193602.png


Attachment 3: Screenshot_from_20190301_194625.png


135

Mon Mar 4 17:22:07 2019 
Ching Pin  Mechanics   comsol modelling  I've updated the material properties to vary with temperature, mainly in the range of 90140 K. Using the parametric sweep function to vary the input power of the heater, we get the eigenfreqencies' dependence on temperature to show up. The fractional dependence of 1.3e5 /K around 123 K matches with what Aaron calculated in this elog entry, which is always a good sign that nothing is going horribly wrong. I've also added the flats to the silicon disc, for better accuracy. See the screenshots below showing the frequency shift with temperature. 
Attachment 1: Screenshot_from_20190304_191328.png


Attachment 2: Screenshot_from_20190304_191441.png


Attachment 3: Screenshot_from_20190304_191537.png


Attachment 4: Screenshot_from_20190304_191556.png


136

Wed Mar 6 09:51:18 2019 
Ching Pin  Mechanics   comsol modelling  So I tried adding the sapphire lens to the comsol model, and I am having teething issues. I can't seem to get the solver to converge, but I'm working on it. 
Attachment 1: Screenshot_from_20190306_093124.png


137

Thu Mar 7 10:10:37 2019 
Ching Pin  Mechanics   comsol modelling  There are no issues with the thermal side of the modeling, the issue seems to be with the structural mechanics side. I'm not sure what I'm doing wrong though, but it just isn't converging. In any case, seeing that this is my last day here, I'll just point out that the version without the lens is saved in cvs/cds/caltech/users/cp/current working model.mph, while the model with the lens is saved in the same folder under the file name testing with lens.mph, using optimus. There is also a small file edition of current working model, with a file name that is self evident. I'll leave it to aaron to upload that to git.
In any case, let me just put down some documentation and thoughts on this model: The physical parameters on the model are generally what we do know of silicon at these temperatures, with the exception of emissitivity, which was randomly given a parameter of .5. The model is currently absorbing 22 mW from the heater and .5 mW from the laser, which implies that the heater should be able to have 45 mW incident on the disc, which would in turn suggest that you would want it to at least dissapate 100 mW to account for the lack of direction from radiation. Because comsol's deposit beam power function does not care for emissitivity, it must be modified in tandem with it. 
Attachment 1: Screenshot_from_20190307_101047.png


3

Tue Apr 2 19:24:07 2013 
Emory Brown  Optics   First steps producing a Finite Element Model to find the internal Brownian noise of a LIGO test mass 
Motivation:
Reduction of Brownian thermal noise in future gravitational wave detectors is of significant interest. It has been suggested that changing the shape of the mirrors used may reduce the Brownian thermal noise. Before we can study how alterations in mirror shape effect Brownian thermal noise, we must be able to generate a finite element model which can compute the Brownian thermal noise in current mirror substrate so that the model may be tested against other calculations of that value using the fluctuationdissipation theorem.
Methods:
I began by constructing a model of a mirror design in COMSOL with the following parameters:
radius (R): 172.5 mm
height: 102 mm
Gaussian beam size (rBeam): 4 cm
TotalForce: 1 N (this is a placeholder force)
loge: log base 10 of e or about 0.4343
PressureCenter: TotalForce*(1/pi)*(rBeam^2)*(1e^(R^2/rBeam^2))^1
A cylinder was constructed using the radius and height above with its flat faces pointing in the z direction, and was given the properties of fused silica which COMSOL had built in.
Our first nondefault boundary condition was applying a force representing the laser beam onto one of the flat faces of the cylinder, surface 3, in the model. The force applied to this face was a Gaussian function PressureCenter*e^(sys2.r^2/rBeam^2) which was specified in COMSOL as a boundary load.
In order to keep the cylinder fixed in space, representing its being connected to wires supporting it and our ignoring violin modes, we need to apply another boundary condition. The simplest possibility would be to force surface 1 in the model, the face of the cylinder opposite the applied force, from moving, but page 5 of Liu and Thorne's paper demonstrates that we should instead apply a force of equal magnitude and opposite direction to that applied to surface 3, but distributed over the volume of the cylinder. To this end, we integrated the Gaussian force applied over the face of the mirror and determined that the total force applied was TotalForce=pi*PressureCenter(1e^(R^2/rBeam^2))rBeam^2/loge. It is more convenient to specify the total force than the pressure at the center of the mirror's face, so solving for PressureCenter we obtain PressureCenter=TotalForce*(loge/pi)*(rBeam^2)*(1e^(R^2/rBeam^2))^1. This opposing force was entered into COMSOL as a body load using the value TotalForce applied over the volume of the cylinder.
COMSOL was instructed to solve for a steady state given the above configuration and returned an error message that “The relative residual (19) is greater than the relative tolerance.” I increased the number of elements in the mesh and the analysis returned lower relative residuals (9.1) using a “finer” mesh. The computer being used did not have enough memory to use a finer mesh structure than that, but the lower relative residual indicates that using a finer mesh may solve the convergence problem.
Ideas for future work:
The simplest possibility is to perform the simulation again using a finer mesh on a computer with more memory and see if we can obtain a solution.
The function used to assign forces at the discrete points in the mesh is continuous, but the points are discrete. Using a discrete Gaussian function to determine the force at each point on the face may be worth trying.
Alternatively, the process could be handled by a Matlab script setup to first run a COMSOL simulation to determine the integrated force on the face of the cylinder, then normalize the central force based on that data such that the desired force is applied to the cylinder face. The script would then solve for the steady state solution.
We could also consider replacing the TotalForce applied over the cylinder's volume by a different boundary condition. Replacing it by a weak spring force on the back face of the mirror has been suggested. I think this is less likely to give good results than the above suggestions, but it may still be worth testing and seeing how its solution compares to values obtained from the fluctuationdissipation theorem.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.

5

Mon Jun 3 21:01:46 2013 
Emory Brown  Optics  General  Analytic Calculation of Thermal Noise due to Brownian Motion using Levin and Thorne's Methods  Motivation:
Reduction of Brownian thermal noise in future gravitational wave detectors is of significant interest. It has been suggested that changing the shape of the mirrors used may reduce the Brownian thermal noise. Before we can study how alterations in mirror shape effect Brownian thermal noise, we must be able to calculate the thermal noise analytically in order to compare it to our finite element model.
Methods:
I made Mathematica notebooks to perform calculation of the thermal noise. The first notebook implemented Levin's method directly.
Sx[f_] := (4*Kb*T/f)*(1  \[Sigma]^2)*1.87322*\[Phi]/(\[Pi]^3*E0*r0)
To test this against Levin's paper, the same values were used as in the paper, such that
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below. Additionally, the value for 100 Hz was explicitly computed and agreed with Levin's value of 8.7*10^40 m^2/Hz computed from his equation 15.
The notebook was modified to perform Thorne's method of computing the thermal noise for a finite test mass. This calculation was performed using equations 28, 35a, 54, 56, and 57. Equations 56 and 57 require the approximation made in eqn 37 which assumes a relatively low mass mirror. According to Liu and Thorne, this is a "rather good approximation for realistic parameter values." Performing the calculations again using this method gave Sx[100 Hz]=7.80081*10^40 m^2/Hz. So, for these parameters, the finite test mass provides a correction factor of about 10%. Another loglog plot of Sqrt[Sx[f]] against f was made using this method.
We can now modify the parameters above to match the values in our finite element model to verify the results of our finite element model.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.

Attachment 1: ComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.png


Attachment 2: DifferencePlot.png


Attachment 3: ThorneComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.nb

(* Contenttype: application/vnd.wolfram.mathematica *)
(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)
(* CreatedBy='Mathematica 9.0' *)
(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
... 2754 more lines ...

6

Tue Jun 4 17:25:11 2013 
Emory Brown  Optics  General  Comparison of First Results from COMSOL Model to the Analytic Solution  Motivation: Before we can attempt to modify the mirror designs to reduce the thermal noise due to Brownian motion in them, we must verify that our model works and that its results match with the analytically calculated ones.
Methods: I had previously constructed a model of a cylindrical mirror in COMSOL: http://nodus.ligo.caltech.edu:8080/COMSOL/3. I updated the values used in the COMSOL model to match those from Levin's paper so that they would match the values I had already computed analytically. These values are listed below:
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
Using the simplest of the boundary conditions I had attempted to implement, fixing the mirror face opposite the applied force, I ran a stationary solver on the model. After the solver had completed, I ran a volume integration of the strain energy in the mirror and obtained U_{max}=1.52887*10^10 J. I also ran a surface integral of the force applied to the surface to confirm that the total force, F_{0}, was the 1N that the COMSOL model had applied to the mirror's face, which it was.
Levin's equations 10 and 12 were combined to give S_{x}(f)=[(K_{b}*T) / (Pi*f)] * [U_{max} / F_{0}^2] * Phi
With the applied force of 1N and the value of U_{max}=1.52887*10^10 J, S_{x}(100 Hz)=2.0157*10^40 m^2 / Hz which agrees to within a factor of 4 with the results of the calculation based on Liu and Thorne's paper which gave a value of 7.80081*10^40 m^2 / Hz. A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below.
Future work:
The obvious next step to take in the project is to attempt to get the better boundary conditions to work, in particular the "gravitational" body load suggested by Liu and Thorne. We noted while working today that a much simpler case where we applied a force of 2N to one surface of a cylinder and an opposing force of 2N in the opposite direction did not converge. It may be worth working with this case and attempting to get it to converge in order to inform how we can make the more complicated case converge. If we can get that case to converge and it agrees with the analytic results, then we will be ready to start varying the relative sizes of the two mirror faces and determining the effect on thermal noise due to Brownian motion.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.

Attachment 1: ComsolComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.png


7

Wed Jun 5 13:31:41 2013 
Matt A  Optics  General  Analytic Calculation of Thermal Noise due to Brownian Motion using Levin and Thorne's Methods  Nice Work Emory,
Can you make a plot with both of the lines on the same plot, with a legend, units, and everything? I'd like to see these compared sideby side.
Quote: 
Motivation:
Reduction of Brownian thermal noise in future gravitational wave detectors is of significant interest. It has been suggested that changing the shape of the mirrors used may reduce the Brownian thermal noise. Before we can study how alterations in mirror shape effect Brownian thermal noise, we must be able to calculate the thermal noise analytically in order to compare it to our finite element model.
Methods:
I made Mathematica notebooks to perform calculation of the thermal noise. The first notebook implemented Levin's method directly.
Sx[f_] := (4*Kb*T/f)*(1  \[Sigma]^2)*1.87322*\[Phi]/(\[Pi]^3*E0*r0)
To test this against Levin's paper, the same values were used as in the paper, such that
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below. Additionally, the value for 100 Hz was explicitly computed and agreed with Levin's value of 8.7*10^40 m^2/Hz computed from his equation 15.
The notebook was modified to perform Thorne's method of computing the thermal noise for a finite test mass. This calculation was performed using equations 28, 35a, 54, 56, and 57. Equations 56 and 57 require the approximation made in eqn 37 which assumes a relatively low mass mirror. According to Liu and Thorne, this is a "rather good approximation for realistic parameter values." Performing the calculations again using this method gave Sx[100 Hz]=7.80081*10^40 m^2/Hz. So, for these parameters, the finite test mass provides a correction factor of about 10%. Another loglog plot of Sqrt[Sx[f]] against f was made using this method.
We can now modify the parameters above to match the values in our finite element model to verify the results of our finite element model.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.


8

Wed Jun 5 13:49:13 2013 
Deep Chatterjee  Optics  General  Difference in the Levin and Liu & Thorne results for thermal noise  The analytical expressions for thermal noise, as calculated by Liu & Thorne and Levin, was plotted as a function of frequency in a log  log plot using Matlab.
The value of the parameters were used from Levin's paper on thermal noise.
phi = 1*10^{7} (loss angle)
r0 = 1.56*10^{2 }(beam radius)
T = 300 K (temperature)
E0 = 7.18*10^{10 } (Young's modulus)
sigma = 0.16 (Poisson ratio)
The value of S_{x}(100) that given in Levin's paper (8.7*10^{40} m^{2} /Hz) while the LiuThorne value is 9.1*10^{40} m^{2} /Hz.

Attachment 1: Levin_throne_comparison.pdf


9

Wed Jun 5 16:56:53 2013 
Matt A.  Optics  General  Difference in the Levin and Liu & Thorne results for thermal noise  Good Work Deep,
Can you include the equations that you used to calculated these expressions?
Quote: 
The analytical expressions for thermal noise, as calculated by Liu & Thorne and Levin, was plotted as a function of frequency in a log  log plot using Matlab.
The value of the parameters were used from Levin's paper on thermal noise.
phi = 1*10^{7} (loss angle)
r0 = 1.56*10^{2 }(beam radius)
T = 300 K (temperature)
E0 = 7.18*10^{10 } (Young's modulus)
sigma = 0.16 (Poisson ratio)
The value of S_{x}(100) that given in Levin's paper (8.7*10^{40} m^{2} /Hz) while the LiuThorne value is 9.1*10^{40} m^{2} /Hz.

The Liu Thorne expression being
S = 4*kb*T/(pi*f) * (1  sigma^2) / (2*sqrt(pi)*E*sqrt(2)*r0) * phi
That of Levin goes as
S = 4*kb*T/f * (1  sigma^2) / (pi^3*E*r0) *I*phi
sigma : Poisson ratio
phi : loss angle
r0 : Beam spot radius
E : Young's modulus
I : This is a sum the value of which is approx. 1.873.22 [Eqn.(A6) of Levin  Internal Thermal Noise] 
10

Wed Jun 5 16:59:51 2013 
Matt A.  Optics  General  Matlab code for calculating the zeros of Bessel functions from GWINC  For calculating the Liu and Thorn U_0 + DU, you need to sum over the first N zeros of the firstorder Bessel function. Unfortunately, Matlab doesn't seem to come with a function to do this.
Rather than reinvent the wheel, we can just use the function used by GWINC. 
Attachment 1: besselzero.m

function x=besselzero(n,k,kind)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% besselzero.m
%
% Find first k positive zeros of the Bessel function J(n,x) or Y(n,x)
% using Halley's method.
%
% Written by: Greg von Winckel  01/25/05
... 79 more lines ...

11

Wed Jun 5 20:39:48 2013 
Deep Chatterjee  Optics  General  Conventional Thermal noise (Sec V) from Liu & Thorne  >The linearized PSD plots are created for the case of thermal noise in finite test mass (Sec. V) of Liu and Thorne.
S = 8*kb*T*phi*(U0 + delta_U)/omega
>The maximum energy due to stress is considered by an infinite sum here. A comparison has been made regarding the convergence of the sum.
>The two plots correspond to the cases of considering 10 and 100 terms in the sum respectively.
>The plot shows that the difference is not much and hence convergence is fast.
>The relative difference is plotted w.r.t S_100, the PSD considering 100 terms in the sum.
>The relative difference goes like abs(S_100  S_10)/S_100, where 10 and 100 represent the number of terms considered in the sum.
>The algorithm used to evaluate the sum involving Bessel functions was the one by GWINC. (http://nodus.ligo.caltech.edu:8080/COMSOL/10). 
Attachment 3: Thorne_thermal_noise.pdf


Attachment 4: Percentage_difference.pdf


13

Thu Jun 6 12:40:07 2013 
Matt A.  Optics  General  Conventional Thermal noise (Sec V) from Liu & Thorne  Good work Deep,
Can you write more on what this is and why you're doing it? We want our elog entries to be easy to understand for anybody who reads it.
Quote: 
The power spectrum of thermal noise in the case of Finite test masses differ slightly from that of the Infinite Test Mass. The expression for the PSD in the infinite test mass case is not a closed form solution but an infinite sum. In this post, a comparison has been made to check how fast does the sum (and as a result, the PSD) converge. The numerical values used for the calculation have been taken from the paper by Levin and that by Liu and Thorne.
>The linearized PSD plots are created for the case of thermal noise in finite test mass (Sec. V) of Liu and Thorne.
S = 8*kb*T*phi*(U0 + delta_U)/omega
>The maximum energy due to stress is considered by an infinite sum here. A comparison has been made regarding the convergence of the sum.
>The two plots correspond to the cases of considering 10 and 100 terms in the sum respectively.
>The plot shows that the difference is not much and hence convergence is fast.
>The relative difference is plotted w.r.t S_100, the PSD considering 100 terms in the sum.
>The relative difference goes like abs(S_100  S_10)/S_100, where 10 and 100 represent the number of terms considered in the sum.
>The algorithm used to evaluate the sum involving Bessel functions was the one by GWINC. (http://nodus.ligo.caltech.edu:8080/COMSOL/10).

I had made some minor errors in the expressions previously while typing the expressions. I have made the corrections and uploaded the new plots in place of the older ones. 
14

Thu Jun 6 17:00:11 2013 
Deep Chatterjee  Optics  General  Comparison of ITM and FTM geometries (FIG. 3 Liu & Thorne)  In section V of their paper on thermoelastic and homogenous thermal noise, Liu and Thorne have corrected the result of power spectral density of a Finite Test Mass from the result of BHV. Although the result of the infinte test mass has a closed form solution, that for the finite test mass is not closed and has to be approximated numerically. From the infinite series (Eq. (56), Eq.(57)) 100 terms were taken to approximate the sum. Considering that the convergence is fast (my last ELOG entry has a plot which shows little difference between considering 10 and 100 terms), 100 terms are a fair approximation.
The spectral density of thermal noise in the finite cavity is slightly lower than the corresponding infinite case. To highlight the difference between the two, they have plotted a quantity C_FTM which is a ratio of the linearized spectral densities of the Finite Test Mass and the Infinite Test Mass.
In this post, the same quantity has been numerically computed and plotted.
The above plot maybe compared to the one given on Fig. 3 of Liu & Thorne. Here r0 is the beam spot radius and CFTM is the quantity mentioned previously. A snapshot of the figure from Liu & Thorne is shown below
The match between the figures seem fair.
Since C_FTM is a ratio, its frequency independent (both S_ITM and S_FTM have 1/f dependence on frequency which cancels on taking the ratio). This is verified in the above plot where plots are made for 10, 100 and 1000 Hz and they fairly coincide. 
Attachment 2: C_FTM.pdf


Attachment 3: Fig3_Liu&Thorn.bmp

19

Wed Jun 12 13:52:44 2013 
Deep Chatterjee  Optics  General  Difference in COMSOL and Analytic solution for TR noise  In this post, I am reporting of the relative difference between the analytic and COMSOL results.
The file by the name 'Main_thermo_refractive.m' found in the SVN was slightly modified to generate the difference plot.
The quantity plotted is the absolute value of the difference between the two results divided by the result given by COMSOL.
The plot is in a semilogx plot. It can be seen the maximum deviation goes to a maximum of 6%.
Also the difference is higher at the higher frequencies.
The plot shows two cases  the case for which 500 terms in the series were summed to get the analytic solution and second
is the case where 1000 terms were summed. The plot does not change much due to the number of terms being summed
showing that the convergence is fast.

21

Mon Jun 17 14:59:20 2013 
Emory Brown  Optics  General  Improved Meshing in a COMSOL Model  In order to make computations more efficient and possibly allow the set of boundary conditions based on Liu and Thorne's suggestion of a gravitational force preventing bulk motion of the mirror (described in http://nodus.ligo.caltech.edu:8080/COMSOL/3) we want to be able to more optimally mesh our structures. In particular, we would like to have a finer mesh on the face of the mirror and especially near the center of the mirror. Doing so will allow us to significantly reduce the total number of elements contained in the mesh while keeping a large number in the regions which require them.
The best way to improve the meshing that I have found so far is to introduce a new geometry and mesh around it. In order to test this method, I constructed an additional cylinder in my model centered at the same location as the mirror. I gave the cylinder the same height as the mirror so that they can be changed by modifying a single variable and set the cylinder's radius equal to the Gaussian beam size of the laser. I then constructed a userdefined mesh composed of two free tetrahedral operations. The first of these I restricted to the smaller cylinder. By selecting a distribution using a fixed number of elements, this region can be meshed fairly uniformly and densely. The second meshing region can be specified again using the inner cylinder's top face as well as the edges running along its side. For this region, it seems optimal to use a predefined distribution type and select a geometric sequence. Using this method provides a much finer mesh in the areas where we want one, without wasting computational power performing more calculations in uninteresting regions. A screenshot demonstrating the results of this method is shown below.
I ran the previous analysis using the simplest boundary condition of a fixed back using this new mesh. The newly constructed mesh contained 7038 elements instead of the previous 7243 and obtained a Umax=1.57985*10^10 J instead of 1.52887*10^10 J. This is about a 3% difference in results. In the next few days, I will run the simulation again on a computer with more RAM using a finer version of the original mesh in order to confirm that the results of the new mesh agree better with it than with the current mesh with a similar number of elements. As a more immediate test, I constructed a mesh of 1940 elements using the new method and obtained a result of Umax=1.56085*10^10 J which is closer to the value obtained using the original meshing technique, though far enough away to indicate that they may not agree, which encourages me to run the original mesh design with more elements.
The next thing to consider is further improving the mesh by increasing the element density near the mirror's front face. From what I have seen it may be more difficult to implement both of these improvements together, in which case we can perform testing to determine which of the two methods provides better computational efficiency. 
Attachment 1: COMSOLMeshing.png


22

Tue Jun 18 15:26:50 2013 
Emory Brown  Optics  General  Attempting to Implement Liu and Thorne's Boundary Conditions  Liu and Thorne demonstrated in their paper that the optimal set of boundary conditions for a problem like the one I am working on is the balance the force applied to one of the mirror faces by an opposing force of equal magnitude distributed throughout the body of the mirror. When we first attempted to implement this boundary condition, we obtained an error from COMSOL indicating that a solution could not be generated as "The relative error is greater than the relative tolerance." I attempted to find a solution to this problem so that we can run our simulations using the correct boundary conditions.
I tried a number of settings changes and tweaks suggested on the COMSOL forums and in official documentation on the error (https://community.cmc.ca/docs/DOC1453), but without success. Most of the suggestions were not relevant to our model, and confusingly when the number of elements in the mesh was increased from 13628 to 24570 the relative residual increased contrary to our initial attempts in http://nodus.ligo.caltech.edu:8080/COMSOL/3. It seems like the next action to take is to attempt to generate a simpler model and attempt reproduce and fix the same error.
I also decided to run both an eigenfrequency and frequency domain solver on our model using improved meshing with 24570 elements. I first ran these solvers using the rigid back constraint and both of them returned Umax=1.49931*10^10 which gives S_{x}(100 Hz)=1.97673*10^40 m^2/Hz which agrees quite well with the solution obtained from our previous simulations which gave S_{x}(100 Hz)=2.0157*10^40 m^2 / Hz. This still significantly differs from the analytically obtained value of 7.80081*10^40 m^2 / Hz using Liu and Thorne's calculation.
I then changed the boundary condition to use Liu and Thorne's gravitational body load boundary condition and both of these methods returned Umax=1.17164*10^9 J which corresponds to S_{x}(100 Hz)=1.54472*10^39 m^2/HZ. This only differs from the analytic result by a factor of 1.98, which is a significant improvement over the previous simulations, but still a much more significant difference than I would hope for given that they now use the same boundary conditions and it seems unlikely that using a stationary state solver will give a different result when both the eigenfrequency and frequency domain solvers agree to 6 digits.
I will spend some time attempting to figure out why the stationary solution is returning the relative residual error. Afterwords, I will start to learn how to use Matlab in conjunction with COMSOL in order to allow us to vary the mirror's size as we want in a systematic manner. Also, Yanbei Chen suggested that we may want to add in a mirror coating later in the project and see how its Brownian noise responds to changes in the mirror substrate shape.
edit: Taking the advice of the error documentation, I also ran a time dependent solver which returned the same result as the eigenfrequency and frequency domain solvers.

27

Thu Jun 20 16:29:38 2013 
Deep Chatterjee  Optics  General  Differences on using a coarser mesh in calculating TR noise  In the code for evaluating TR noise, the parameter for automatic meshing (which can take values from 1 to 9  1 being the finest and 9 the coarsest mesh) was changed from 1 to 5 and 9.
The computation time surprisingly doesn't differ much. The maximum relative error also stays close to 6% as in the case of the finest mesh. The relative difference for two cases for a mesh
number of 5 and 9 are shown in the attachment.
The conclusion being that it is better to stick to the mesh number 1 which is the finest mesh.

29

Sun Jun 23 14:29:04 2013 
Koji  Optics  General  Differences on using a coarser mesh in calculating TR noise  What about the dependence on the size of the time slice?
It is unclear what the unit fo the frequency. Hz? ot logHz (which should be expressed as 10^1, 10^2, ...)?
Can you combine the result of the caluculation and the error together?
i.e. Combine the plot similar to http://nodus.ligo.caltech.edu:8080/COMSOL/4 and your plot with the horizontal axies lined up. 
30

Mon Jun 24 16:02:43 2013 
Emory Brown  Optics  General  Rebuilding a minimal version of our COMSOL model  I rebuilt our COMSOL model in order to see if I could find any errors in it and to clean out unnessesary features and parameters. The model returned the same results as before, and I did not find anything wrong with it while rebuilding it. Running an Eigenfrequency solver gives the first real mode as a sheer mode, indicating that a normal meshing may be insufficient, but finer meshes result in an out of memory during LU factorization error despite having just rebooted the computer used and having several gigs of unused RAM. I have also submitted a ticket to COMSOL support asking for assistance with the relative residual error, but the initial response I recieved was unhelpful. Hopefully I get a more helpful response soon. I have started working through a LiveLink for Matlab tutorial (http://www.rz.unikarlsruhe.de/~hf118/tutorial_comsol_matlab_livelink.pdf) I found from a third party and reading example Matlab code, since I am not familiar with Matlab. I will continue reading through this documentation for the next day or so, then hopefully will recieve useful help from the COMSOL support team. If so, I can start working on using Matlab Livelink to parameterize our model and vary the radius of the back face of the mirror to replicate Steve Penn's results. Otherwise, I may try doing so with the increase of relative tolerance which caused the stationary solution to converge last week. 
31

Mon Jun 24 18:00:26 2013 
Deep Chatterjee  Optics  General  The problem of TE Noise calculation for a beam passing through a cylindrical substrate  In this post, I mention of my immediate plan for the week and also describe the problem I am looking at in brief.
In the calculation of TR noise for the beam passing through a cylindrical substrate, Heinert has considered thermal fluctuations
by means of an sinusoidally oscillating heat source, scaled according to the intensity profile of the beam, present along the length of the cylinder.
The heat equation is then solved to evaluate the temperature field. The gradient of temperature is used to calculate the work dissipated.
The work dissipated then is used to find the spectral density according to the Fluctuation Dissipation Theorem.
On the other had, the case of TE noise is considered in the work by Liu and Thorne where the beam reflecting from one of the faces is
considered. In their work, they have used an oscillating pressure, scaled according to the beam, at the face. The equation of stress balance
is solved to get the strain field, the strain results in heating and as a result, a temperature gradient. The work dissipated is calculated using
the gradient of temperature and the spectral density is calculated using the FDT.
Now, in case of Heinert's work, it is to be noted that the heating term is the one that contains the TR coefficient, beta = dn/dT. This is the where
the TR noise that we are looking into comes in the picture. Liu and Thorne's work on the other had, never has an explicit heat source. In their case,
it is the strain that generates the heat implicitly. They have related the expansion with the temperature perturbation in their paper from an expression
in Landau Lifshitz's, Theory of Elasticity. The important point to note is that the physical parameters that characterize the TE noise like coefficient of
linear expansion, alpha, or Poisson ratio, sigma, come into the picture through this relation. Another point to note is that the expression
for the work dissipated ( W_diss ) uses the gradient of temperature( It is the same formula that Heinert has used). This expression is derivable from the heat
equation. Thus, one could have also done the exercise by injecting the right heat source and solving the heat equation instead, since its ultimately
fluctuations in temperature that cause these noises TR or TE.
Our problem is to evaluate TE noise for a beam that passes through the cylindrical substrate instead of reflection off the face. It is suspected that using an
oscillating pressure on the surface will not be the correct approach since the beam is going through the material and not just reflecting from the surface. We
want to solve it by means of a heat injection, as done in Heinert, calculating the gradient of temperature, the work dissipated and then the spectral density. It is realized that the
heat source should be oscillating but the correct coefficients is what is undetermined i.e. we realize the heat source is q_dot = [  ] * cos(omega * t) . In case
of Heinert it is at this point that the 'beta' comes in. However, in the TE case this is not yet determined. The literature doesn't deal with the case of beam goin through
a material. The equations in Heinert must be looked at more deeply to realize how the 'beta' comes in and then drawing an analogy, we may be able to figure out the
right heat source for the TE noise case.
Any comments on references, the approach that should be taken, or any thoughts on the problem is most welcome. 
32

Tue Jun 25 15:16:59 2013 
Arnaldo Rodriguez  Optics  General  Setting Up Looped Simulation for PID Controller  After setting up a COMSOL model that includes the heat flux from the laser and the ring heater, I've made the model solve over time manually by performing the solution process over a loop in MATLAB.
This allows for the future insertion of the PID controller object in the solving process, and the dynamic manipulation of the applied heat loads.
The following is an automatically generated plot of defocus effects as a function of time at 1 beam radius (54 mm), included in the program, with only the ring heater turned on in a tophat emission configuration and the total power being 5 watts.
The linear projection values needed to calculate the defocus effect are extracted directly from COMSOL, with no output files required through the use of the mphinterp command.
The behavior appears to be physical and is of the correct order of magnitude.
I've also attached the code that produced it for verification. (It requires COMSOL+Matlab to run.)

Attachment 2: ITMPlusRing.m

% function out = model
%
% ITMPlusRing.m
%
% Model exported on Jun 24 2013, 16:24 by COMSOL 4.3.0.151.
deltat = 600;
totalt = 3600*24;
t = [0:deltat:totalt];
omega = 54/1000;
... 275 more lines ...

33

Tue Jun 25 15:18:41 2013 
Emory Brown  Optics  General  A new boundary condition to attempt to avoid the relative residual error  I contacted COMSOL support about our difficulties with the relative residual error and was told "Typically for solid mechanics this error is because you have not constrained the body enough. So any solution + a rigid body transformation is also a solution. To remove this nonuniqueness, you need a solutions modulo the rigid body transformations. So try and constrain the body somehow." This is what the gravitationally balanced body load was supposed to do, but using the stationary solver has been nonconvergent and the eigenfrequency solver has generated odd modes, with the first real modes being sheer modes. These outputs indicate a problem with our model. Matt noticed yesterday that the meshing in our model was slightly nonsymmetrical, but we initially dismissed it as not being significant since it was a minor difference which we did not think could account for the large errors we are facing. At further consideration though, if some asymmetry arising from any source, even a numerical or rounding error, were present, that could cause a slight rotation of the object which would cause the object to experience a torque and minor sheer force. The stationary solution would not converge in this case, but the eigenfrequecy solution might, and if it did it would make sense to see sheer modes for some of its low frequency eigenmodes, as we have. One possible solution to this problem is to change the way we mesh the material to ensure a symmetrical distribution of nodes in the xy plane, probably by extruding lower dimensional meshed systems into our model. I am unsure if we would be able to implement this solution once we start to change the size of the radii of the object's faces. An alternate solution is to find another set of boundary conditions which should be equivalent to the gravitational body load constraint, but which are stable relative to minor perturbations of the system's conditions. I think that I have found another set of boundary conditions which should work and not be too difficult to implement in COMSOL.
The sides of the object should not move in the direction orthogonal both to their displacement from the center of the circular plane of the object even with them in the z direction or in the z direction, so the edge displaced from the center in the y direction should be fixed in the x direction.
The object's center of mass should not move which because of the previous condition can be reduced to not moving in the z direction.
I think that these boundary conditions should either be compatible with, or replace Liu and Thorne's boundary conditions in our model. I am going to spend some time attempting to implement these boundary conditions, see if they converge, and see if adding Liu and Thorne's gravitationally balanced force with them makes any difference, either in results (which it shouldn't) or the amount of time required to run. I am not sure how long this will take to implement, but I don't expect it to take more than a day or two. 
34

Wed Jun 26 13:52:06 2013 
Arnaldo Rodriguez  Optics  General  Verifying Relative Error in Defocus for Regular and ManualLoop Simulations  To verify the validity of the solutions produced by the manual simulation, I've calculated the relative error between the results from the manual code and the results produced by COMSOL normally.
The plot for the relative error in the defocus at r = 0 and r = 54 mm is shown below, in the case where only the ring heater is turned on at a total power of 5 W.
The following indicates that the maximum error is less than 0.01% (in percent error format).

35

Wed Jun 26 15:23:55 2013 
Arnaldo Rodriguez  Optics  General  Solving Time per Loop in Manual Dynamic Simulation  I've attempted to determine the solving time per loop as a function of the simulation time, in an attempt to identify any trends in the solving time for a constantly dynamic load.
The following is a plot of the solving time per loop as a function of the simulation time for a load which is constantly dynamic (sinusoidal in time, in this particular run).
The mesh size is normal (default in COMSOL) with heat fluxes from both the beam and the ring heater (as in the real case).
It is difficult to identify any particular behavior or trend due to the large amount of "noise" other than a trivial general increase after ~4 s. Mesh quality does not appear to influence solving time per loop significantly.
Work must be done to reduce the total solving time for the simulation, which in this case amounted to 18 and a half minutes (1.1109e+03 seconds).

36

Wed Jun 26 16:46:32 2013 
Emory Brown  Optics  General  A convergent stationary solution  By adding additional constraints as described in the previous eLog to the model I was able to get a convergent solution using the stationary solver. I cleared meshes, solutions, and history and uploaded the file in its current form. Using the stationary solver, we get Umax=1.51646*10^10 J instead of the 1.52887*10^10 J obtained using the rigid back constraint. This is less than a 1% difference, so it is still off from the analytic solution by a factor of ~4. I will look through the model and analytic calculations again and see if I can find anything which could contribute more than a few percent difference. 
Attachment 1: ConvergingModelMinimal62613.mph

37

Thu Jun 27 11:26:18 2013 
Emory Brown  Optics  General  Missing factor of 4 and a comparison of analytic and simulated Brownian noise  In our calculation of Sx(f) from Umax prior to this, we were missing a factor of 4. The correct formula from Liu and Thorne's formula 58 is Sx[f] := 4*Kb*T*Umax*Phi/Pi*f*F0^2). When we correct this formula, we get results which agree fairly well with the analytic results. The Mathematica script attached was used to perform both calculations of the thermal noise spectrum (using SxComsol and SxFTMThorne).
SxCOMSOL(100 Hz) = 7.99735*10^40 m^2/Hz
SxFTMThorne(100 Hz) = 7.80081*10^40 m^2/Hz
These results differ by about 2.5%. We also need to verify that the result converges for increasing mesh size.
Number of Elements

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

I attempted to run the simulation with 24360 elements in the mesh, but the computer I was running it on repeatedly crashed.
edit: Running on a more powerful computer I got the following additional values.
Number of Elements

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

24181

1.52169

61772

1.52281

Given these additional values, it appears that our simulation will converge to a value for increasing mesh size and that it agrees fairly well with the analytic result. 
38

Thu Jun 27 15:06:13 2013 
Arnaldo Rodriguez  Optics  General  PID Function in Manual Simulation  I have inserted a rudimentary PID function into the manual simulation code as a way to test whether or not the PID function is changing the defocus values in the desired manner.
I am currently determining the ratio of ring heater power to the steadystate defocus as a way of measuring the scale of the response.
This ought to give a good way of measuring the scale needed to convert the calculated actuator response into an actual load.
I've attached the rudimentary code below. (The actuator isn't feeding into the heater at the moment, but inserting the "actuation" variable into the load expression is all that is required.)

Attachment 1: ITMPlusRingPlusPID.m

% function out = model
%
% ITMPlusRing.m
%
% Model exported on Jun 24 2013, 16:24 by COMSOL 4.3.0.151.
u = zeros(1,2);
e = zeros(1,3);
mtime = 0;
deltat = 600;
... 292 more lines ...

Attachment 2: PID.m

function [u, e] = PID(s, u, e, Ki, Kp, Kd, deltat, s_target)
u(2) = u(1);
e(3) = e(2);
e(2) = e(1);
e(1) = s  s_target;
u(1) = u(2) + Ki*e(1)*deltat + Kp*(e(1)  e(2)) + Kd*(e(1)  2*e(2) + e(3))/deltat;
end

39

Thu Jun 27 17:12:33 2013 
Arnaldo Rodriguez  Optics  General  Defocusing in PControlled Manual Simulation  To test the effectiveness of the PID loop, I've only given a nonzero value to the K_p constant and varied it over a series of runs to verify the response by the controller.
The following is a plot of the baseline compared to the Pactive simulations for different proportionality gain constants. I've yet to determine whether all the trends in behavior are physical.
Each different case has an evident droop (steadystate error), and seem to be destabilizing as K_p increases. However, most cases reach steadystate quicker than the baseline.
It should be noted that the ring heater is not correcting any effects from the beam, as it is turned off. The ring heater is correcting itself; or rather, its own initial heat pulse of 5 W * delta_t (600s) = 3000 J before the controller becomes active. 
40

Fri Jun 28 11:54:10 2013 
Emory Brown  Optics  General  Maintaining a constant volume while modifying the radii of the test mass faces  When changing the dimensions of the mirror, we want to keep a constant volume. To do so, we must modify the height of the test mass. The area of the test mass if R is the radius of the front mirror face and r the radius of the back mirror face is (1/3)*pi*h*(R^2+R*r+r^2)=area, so using the values currently in our model, R=r=0.17 m and h=0.2 m, area=0.0181584055 m^3=pi*0.2*0.17^2 m^3, so to ensure that the volume of the test mass remains constant, h=0.01734/(R^2+R*r+r^2). Note that this does give 0.2 m when R=r=0.17 m. 
42

Mon Jul 1 13:01:58 2013 
Arnaldo Rodriguez  Optics  General  Brief Parameter Study on Pactive Control Loop  The following is a plot of the defocus as a function of time with only the proportional gain being nonzero. Only three cases were run before the program crashed.
The dotted line indicates the setpoint.
One can clearly see the loop becomes unstable and oscillates out of control at a proportional gain of 240, and can see the characteristic steadystate offset of the Pcontroller from the set point develop in the other cases. 
43

Mon Jul 1 17:31:40 2013 
Arnaldo Rodriguez  Optics  General  Parameter Study on Sampling Time for Pactive Loop  The following is a plot of the defocus as a function of time for a fixed proportionality gain of 150, and different listed sampling times.
The discrepancy between the behaviors is enough to warrant a proper mathematical explanation; the max. overshoot, for example, is inversely proportional to the sampling rate.

44

Tue Jul 2 13:32:30 2013 
Emory Brown  Optics  General  Using Matlab to run our COMSOL simulations  We want to be able to test a number of values for the ratio of the radii of the frustum which will be much easier once we can get Matlab to automate running our COMSOL model. I have been able to output the COMSOL model to a .m file and was able to figure out a way to export the Umax value we want to Matlab (see the attached .m file).
The model ran in Matlab after debugging it for a short while and returned the same result as the COMSOL model had by itself. I setup a script to loop over some values of RatioRadii and feed them to the model. The first value it fed to the model worked and returned 1.2099e10 J for RatioRadii=0.1, but for RatioRadii=0.2 the model had a
meshing error. I am going to have to see if there is some way we can modify the mesh to fix this error. When I briefly looked into it today, I found that for several RatioRadii which I tried, the object could be meshed for certain finenesses, but not others. I am going to spend tomorrow working on the first progress report, but after that is done I will
work on fixing this meshing error. 
Attachment 1: model_7_1_13.m

function [Umax] = model(ratioRadii)
% ratioRadii must be a scalar
% model_7_1_13.m
%
% Model exported on Jul 1 2013, 17:08 by COMSOL 4.3.2.189.
import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
... 220 more lines ...

Attachment 2: looper.m

function [RatioRValues,UmaxValues] = looper
% Iterates over a range of RatioRValues
RatioRValues=[0.1:0.1:10]
UmaxValues=[]
for r=RatioRValues,
UmaxValues = [UmaxValues model_7_1_13(r)];
UmaxValues
end
... 4 more lines ...

45

Wed Jul 3 00:04:46 2013 
Koji  Optics  General  Using Matlab to run our COMSOL simulations  Please use zipped files instead of uploading extremely long text files on the elog.
Or use wiki or svn instead of elog for uploading them.
And appropriate word wrapping so that we can read the entries. Please. 
46

Thu Jul 4 15:23:47 2013 
Deep Chatterjee  Optics  General  Attempt to compute TE noise following Liu and Thorne  Liu and Thorne, in their paper have analytically solved for the TE noise spectra considering the test mass to be a half infinite space.
An attempt was made to get similar results using a COMSOL model. In this post, the rough profile of the PSD plot is given and is compared
to Liu and Thorne's result.
Steps followed by Liu and Thorne:
> Oscillating pressure, scaled according to beam profile, is applied one of the surfaces.
> The stress balance equation is solved borrowing algorithm from Landau Lifshitz.
> The stress and the temperature perturbation is related using an equation given in Chapter 1 of LL
> The W_diss is calculated from the gradient of the temperature.
> PSD is calculated using FDT due to Levin
> The resulting plot looks like the following  ***Only the frequency dependence is plotted here leaving out the prefactors which will
appear just as an additive constant in the figure. The aim was to see how the plot of the simulation looks like.
Steps followed in the COMSOL model:
> 2D axisymmetric model was used since the applied was applied radially symmetrically.
> The dimensions of the test mass(cylinder) was defined larger as compared to the beam spot radius.
> The physics added was Solid Stress and Heat Transfer in Solids.
> Thermal expansion was enabled in the Solid Stress physics
> Pressure work was enabled in Heat Transfer physics which takes into consideration heat due to stresses i.e. adds this heat source to the heat eqn
> Heat insulation wsa applied to all surfaces.
> All surfaces apart from the one where pressure was applied was held fixed to mimic infinity
> A custom mesh was put.
> The model was exported as '.m'
> It was called for frequency range 10^(0 : 0.25 : 4)
> Following Liu and Thorne, the temperature gradient was integrated over the volume using mphint2() function in Matlab
> The time average was taken over the second half of the cycle to avoid transient which are assumed to die off after a couple of initial cycles.
> All the prefactors were ignored while plotting
> The plot shown below is integral{grad_T}^2} / omega^2 which is proportional to PSD
> We just wanted to see if at all we got anything looking similar to the analytic result of Liu and Thorne (shown above)
The above plot is when the radius is 100 times the beam spot and height of the cylinder is 200 times the beam spot
Other values of radius and height were also tried. The one below shows similar plot as above but for radius 10 times
the beam spot and height 20 times the same
While the one below is for radius = 1000 * beam spot and height = 2000 * beam spot
Here, it is seen that the straight line with a negative slope is seen to some extent.
However, the prefactors, when put in, will show the actual PSD.
Any comments or suggestions related to the model construction in COMSOL or otherwise is welcome. 
Attachment 1: Liu&Thorne_type.bmp

Attachment 5: 2D_almost_infinte_TE.mph

47

Fri Jul 5 14:06:24 2013 
Deep Chatterjee  Optics  General  Attempt to compute TE noise following Liu and Thorne 
Quote: 
Liu and Thorne, in their paper have analytically solved for the TE noise spectra considering the test mass to be a half infinite space.
An attempt was made to get similar results using a COMSOL model. In this post, the rough profile of the PSD plot is given and is compared
to Liu and Thorne's result.
Steps followed by Liu and Thorne:
> Oscillating pressure, scaled according to beam profile, is applied one of the surfaces.
> The stress balance equation is solved borrowing algorithm from Landau Lifshitz.
> The stress and the temperature perturbation is related using an equation given in Chapter 1 of LL
> The W_diss is calculated from the gradient of the temperature.
> PSD is calculated using FDT due to Levin
> The resulting plot looks like the following  ***Only the frequency dependence is plotted here leaving out the prefactors which will
appear just as an additive constant in the figure. The aim was to see how the plot of the simulation looks like.
Steps followed in the COMSOL model:
> 2D axisymmetric model was used since the applied was applied radially symmetrically.
> The dimensions of the test mass(cylinder) was defined larger as compared to the beam spot radius.
> The physics added was Solid Stress and Heat Transfer in Solids.
> Thermal expansion was enabled in the Solid Stress physics
> Pressure work was enabled in Heat Transfer physics which takes into consideration heat due to stresses i.e. adds this heat source to the heat eqn
> Heat insulation wsa applied to all surfaces.
> All surfaces apart from the one where pressure was applied was held fixed to mimic infinity
> A custom mesh was put.
> The model was exported as '.m'
> It was called for frequency range 10^(0 : 0.25 : 4)
> Following Liu and Thorne, the temperature gradient was integrated over the volume using mphint2() function in Matlab
> The time average was taken over the second half of the cycle to avoid transient which are assumed to die off after a couple of initial cycles.
> All the prefactors were ignored while plotting
> The plot shown below is integral{grad_T}^2} / omega^2 which is proportional to PSD
> We just wanted to see if at all we got anything looking similar to the analytic result of Liu and Thorne (shown above)
The above plot is when the radius is 100 times the beam spot and height of the cylinder is 200 times the beam spot
Other values of radius and height were also tried. The one below shows similar plot as above but for radius 10 times
the beam spot and height 20 times the same
While the one below is for radius = 1000 * beam spot and height = 2000 * beam spot
Here, it is seen that the straight line with a negative slope is seen to some extent.
However, the prefactors, when put in, will show the actual PSD.
Any comments or suggestions related to the model construction in COMSOL or otherwise is welcome.

Here is a plot showing the profiles of Liu and Thorne and the COMSOL simulation. Since the constants were not taken into account, one expects the plot to be shifted by a constant equal to the logarithm of the prefactor that has been ignored.
The green 'o' is just a 1/omega^2 plot in the loglog scale which is how Liu and Thorne's result go (this is not linearized PSD, but is doesn't matter as long as we want to see the profiles). The other plot in blue is the COMSOL simulation.
The plot below is a difference between them. One expects it to be a constant but there are a few jitters. Once again, it is emphasized that the numerical values do not give any information related to the actual difference which will only be calculable
once the correct physical constants are put in the formula. I will post another entry once that is done

Attachment 2: L&T_COMSOL_difference.bmp

48

Mon Jul 8 15:47:23 2013 
Emory Brown  Optics  General  Meshing error when varying geometry  I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long. 
49

Mon Jul 8 18:39:31 2013 
Matt A.  Optics  General  Meshing error when varying geometry 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

Great Emory,
Can you put these in a plot so it's easier to read? 
50

Mon Jul 8 19:04:35 2013 
Matt A.  Optics  General  Meshing error when varying geometry 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

When we talked today, I thought that these frequencies were a little low, so I ran a really simple COMSOL model of a cylinder with radius 0.17 m and height 0.2 m (attached).
I solved for the first 10 modes (the first 6 being trivial), and the frequencies I got are listed below:
Mode # Freq [Hz] Type
7 5950 Butterfly
8 5950 Butterfly
9 8106 Drumhead
10 8226 Shear
These are more in the frequency range that I expected.
**What material are you using? I used fused silica.
**I used 'finer' meshing. Perhaps your more complicated meshing is screwing with the solutions, although I'd expect it to make it stiffer, not softer.
**Are you applying force during the eigenvalue calculation? That might make it softer.
**Perhaps some of your constraints are moving the modes around. My model has zero constraints: it's just a cylinder floating in space.
Lets talk about this tomorrow.

Attachment 1: SimpleCylinder.mph

51

Tue Jul 9 14:45:04 2013 
Emory Brown  Optics  General  Meshing error when varying geometry  I tried removing the constraints and applied force and was able to get close to the same frequencies for the butterfly modes as the simple model. By also removing the interior cylinder from our model I was able to restore the butterfly modes as well. I can see why removing the constraints and forces would make a difference, but I don't understand how the interior region which was only used to aid in making the meshing in the central region more dense would have such a strong effect. With that removed as well, we get the following results using a finer mesh (plot attached):
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz) Any lower modes (Hz)
0.1 8777.8 7253.3 3059.7
0.3 5636.9 7264.7 4202.3
0.5 6068.5 7308.0 5980.9
0.75 7341.5 7183.0 no
1 8110.2 5951.3 no
1.1 7302.8 5308.8 no
1.3 5744.3 4087.9 no
This data is more in line with what we would expect, and indicates that the optimal region is probably between ratioR=0.6 and ratioR=0.8. I will run the looping script we designed to map out how Umax changes when varying ratioR, particularly in this interesting region.
This line was used to generate the plot from data.
figure(1); plot(ratio,drum); hold on; plot(ratio,butterfly,'rs'); plot(ratio(1:3),others,'g'); hleg1=legend('Drumhead','Butterfly','Other'); title('Eigenfrequencies of mode types for different ratioR values'); xlabel('ratioR'); ylabel('lowest eigenfrequency'); shg;
Quote: 
Quote: 
I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge. The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully. I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned. This region of sizes which throw errors varies for different values of ratioR. By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.
I also noted our COMSOL model had an error: ratioR was applied to the wrong face. This was corrected as was the .m file.
I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:
ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)
0.1 3624.8 3802.4
0.3 4375.9 3785.4
0.5 4366.5 3650.0
.75 3709.0 3041.85
1 3000 2186.3
1.1 2736.2 1893.0
1.3 2247.6 1426.4
The mirror's front face is 0.17 [m]. Its Gaussian beam size is 0.0156 [m]. Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

When we talked today, I thought that these frequencies were a little low, so I ran a really simple COMSOL model of a cylinder with radius 0.17 m and height 0.2 m (attached).
I solved for the first 10 modes (the first 6 being trivial), and the frequencies I got are listed below:
Mode # Freq [Hz] Type
7 5950 Butterfly
8 5950 Butterfly
9 8106 Drumhead
10 8226 Shear
These are more in the frequency range that I expected.
**What material are you using? I used fused silica.
**I used 'finer' meshing. Perhaps your more complicated meshing is screwing with the solutions, although I'd expect it to make it stiffer, not softer.
**Are you applying force during the eigenvalue calculation? That might make it softer.
**Perhaps some of your constraints are moving the modes around. My model has zero constraints: it's just a cylinder floating in space.
Lets talk about this tomorrow.


Attachment 1: Eigenfrequency.png


52

Wed Jul 10 14:15:20 2013 
Emory Brown  Optics  General  Comparison of Umax and the frequency of the lowest eigenmode against ratioR  I was able to run the looper program designed last week over the interval ratioR=[0.5, 0.9] with a step size of 0.01 and obtained the plot below with a minimal value of Umax=1.4723e10 J obtained for ratioR=0.71, which is about 3.2% less than the value of Umax=1.5202e10 J when ratioR=1. This is not very significant, but the shifting upward of the frequencies of the lowest eigenmodes as described in the previous post is. Additionally, these curves both have their optimal points in the same region, with a ratioR value of about 0.7.
I also ran a higher resolution test to see how the frequency of the lowest real eigenmode responded to changes in ratioR over this region. The plot generated is also attached with a maximal (most favorable) value of 7210.0 Hz at ratioR=0.74.
Data from both of these runs is included in the attached tarball stored as txt files delineated by tabs.
The next few things to try are substituting silicon for fused silica as the material in the model and seeing how the curves change and seeing the response from reversing the orientation of the optimal test mass (ratioR~0.72). I expect that the orientation change will result in much higher noise and a lower minimal real eigenfrequency. 
Attachment 1: UmaxValuesVsRatioR.png


Attachment 2: FrequencyLowestEigenmodes.png


Attachment 3: data.tar.gz

