ID 
Date 
Author 
Type 
Category 
Subject 
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.

83

Sun Jun 22 23:44:56 2014 
Sam Moore  Optics   Going through Heinert's 'TR noise of cylindrical test masses' paper  At this point, my are goals are to 1) convert the timedependent heat equation into stationary, complex form, 2) use the Levin approach to calculate the TR noise given this stationary PDE, and 3) verify the results in COMSOL
I have looked at the Heinert paper and converted the timedependent partial differential Heat equation to a stationary, complex one. I did this by assuming a temperature distribution of the form T(\vec{r},t) = T_0 ( \vec{r} ) e^{i \omega t }, where \omega is the frequency of oscillation of the input heat field: q( \vec{r} , t) = q( \vec{r} ) e^{i \omega t}. Hence, I am assuming the steadystate solution of the temperature field. From this ansatz, the PDE becomes
C_p T(\vec{r}) + \frac{ i \kappa}{ \omega } \nabla^2 T ( \vec{r} ) = q ( \vec{ r} ).
This complex conversion already seems to have been done in the paper, since the stationary solution of the form
T(r, z) = \sum_{n = 0}^{\infty} \sum_{m = 0}^{ \infty} T_{n m} J_0 (k_n r) \cos ( \ell_m z) is assumed.
k_n and \ell_m are obtained from the adiabatic boundary conditions. Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12. ( although I'm not quite sure how to take the second derivative of a zerothorder bessel function).
I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules. At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

84

Mon Jun 23 11:40:17 2014 
Sam Moore  Optics   Going through Heinert's 'TR noise of cylindrical test masses' paper 
Quote: 
At this point, my are goals are to 1) convert the timedependent heat equation into stationary, complex form, 2) use the Levin approach to calculate the TR noise given this stationary PDE, and 3) verify the results in COMSOL
I have looked at the Heinert paper and converted the timedependent partial differential Heat equation to a stationary, complex one. I did this by assuming a temperature distribution of the form T(\vec{r},t) = T_0 ( \vec{r} ) e^{i \omega t }, where \omega is the frequency of oscillation of the input heat field: q( \vec{r} , t) = q( \vec{r} ) e^{i \omega t}. Hence, I am assuming the steadystate solution of the temperature field. From this ansatz, the PDE becomes
C_p T(\vec{r}) + \frac{ i \kappa}{ \omega } \nabla^2 T ( \vec{r} ) = q ( \vec{ r} ).
This complex conversion already seems to have been done in the paper, since the stationary solution of the form
T(r, z) = \sum_{n = 0}^{\infty} \sum_{m = 0}^{ \infty} T_{n m} J_0 (k_n r) \cos ( \ell_m z) is assumed.
k_n and \ell_m are obtained from the adiabatic boundary conditions. Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12. ( although I'm not quite sure how to take the second derivative of a zerothorder bessel function).
I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules. At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

I carried out the calculations for Heinert's simple model in greater detail and was able to obtain equation 12 (with the exception of a plus sign instead of a minus sign in the denominator; I guess the Fourier transform method gives a minus sign. In the end, that doesn't matter, since the modulus square yields the same result in the denominator.) I obtained the same spectral density S_z (\omega), with the exception of a missing factor of R^2. I don't know where this R^2 went, or why it is that the spectral density is associated with the extrinsic variable z (instead of r, say). It's supposed to be a "homogeneous [noise] readout along the z direction." 
85

Mon Jun 23 16:10:17 2014 
Matt A.  Optics   Going through Heinert's 'TR noise of cylindrical test masses' paper 
Quote: 
At this point, my are goals are to 1) convert the timedependent heat equation into stationary, complex form, 2) use the Levin approach to calculate the TR noise given this stationary PDE, and 3) verify the results in COMSOL
I have looked at the Heinert paper and converted the timedependent partial differential Heat equation to a stationary, complex one. I did this by assuming a temperature distribution of the form T(\vec{r},t) = T_0 ( \vec{r} ) e^{i \omega t }, where \omega is the frequency of oscillation of the input heat field: q( \vec{r} , t) = q( \vec{r} ) e^{i \omega t}. Hence, I am assuming the steadystate solution of the temperature field. From this ansatz, the PDE becomes
C_p T(\vec{r}) + \frac{ i \kappa}{ \omega } \nabla^2 T ( \vec{r} ) = q ( \vec{ r} ).
This complex conversion already seems to have been done in the paper, since the stationary solution of the form
T(r, z) = \sum_{n = 0}^{\infty} \sum_{m = 0}^{ \infty} T_{n m} J_0 (k_n r) \cos ( \ell_m z) is assumed.
k_n and \ell_m are obtained from the adiabatic boundary conditions. Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12. ( although I'm not quite sure how to take the second derivative of a zerothorder bessel function).
I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules. At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

Hi Sam,
I know it seems a bit distracting, but it's sometimes nice to write up your elog entries in a more readable form. I've reedited your comments to make them more easily read.
I used http://www.sciweavers.org/freeonlinelatexequationeditor to convert your latex to png files for easy reading.
At this point, my are goals are to:
1) convert the timedependent heat equation into stationary, complex form,
2) use the Levin approach to calculate the TR noise given this stationary PDE, and
3) verify the results in COMSOL
I have looked at the Heinert paper and converted the timedependent partial differential Heat equation to a stationary, complex one.
I did this by assuming a temperature distribution of the form:
where ω is the frequency of oscillation of the input heat field:
Hence, I am assuming the steadystate solution of the temperature field.
From this ansatz, the PDE becomes:
This complex conversion already seems to have been done in the paper, since the stationary solution of the form:
is assumed. k_{n} and ℓ_{m} are obtained from the adiabatic boundary conditions.
Moreover, if I plug this form for T(r, z) into the stationary PDE, I get Heinert's eq. 12.
( although I'm not quite sure how to take the second derivative of a zerothorder bessel function).
I've started to try out the cylindrical geometry in COMSOL, using the structural mechanics and heat transfer modules.
At this point, I'm not quite sure how to specify the PDE that I'm interested in, nor am I sure about how to specify the boundary conditions.

91

Sat Jul 5 13:04:32 2014 
Sam Moore  Optics   Heinert Model TR Noise Verification  Agreement with Heinert's paper for cylindrical TR noise has now been achieved. Using the stationary state assumption to calculate the temperature profile, the computation time was reduced compared to the previous timedependent approach. Here are the plots showing the agreement. I have shown the plots for a 1D axisymmetric model, in addition to a full 3D model in COMSOL. Both give the same result.
What went wrong? In the 1D axisymmetric case, it turns out that COMSOL has the incorrect cylindrical coordinate Laplacian for the coefficient form PDE interface. I corrected for this by expanding the Laplacian with the product rule, giving an extra kappa/r ''convection coefficient'' in the PDE interface. Meanwhile, the 3D case worked from the beginning, since its Laplacian was in cartesian coordinates.
The next goal is to use a stationary state method to calculate TE noise for a test case.

Attachment 1: threeD_cylinderTRnoise.png


Attachment 2: oneDcylinderTRnoise_copy.png


92

Mon Jul 7 19:47:00 2014 
Koji  Optics   Heinert Model TR Noise Verification  How close are these FEA calculations with the analytical values?
Can you plot residual too? (Put analytical values, 1D, abs(1D  analytical), 3D, and abs(3D  analytical) all together.) 
94

Fri Jul 11 10:58:18 2014 
not Koji  Optics   Heinert Model TR Noise Verification 
Quote: 
How close are these FEA calculations with the analytical values?
Can you plot residual too? (Put analytical values, 1D, abs(1D  analytical), 3D, and abs(3D  analytical) all together.)

Here are the plots with their fractional residuals: abs(S_COMSOL  S_analytical)/S_analytical 
Attachment 1: heinert_analyticalTest_residual_threeD.eps


Attachment 2: heinert_analyticTest_residual_oneD.eps


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


77

Wed Nov 13 23:55:44 2013 
Chris Couste  Optics  Analysis  Optical Mount Vibrational Analysis  Project: Vibrational Analysis of Optical Mounts
Goal: Use COMSOL to run finite element analysis on simplified models of different types of optical mounts available to us, in order to find which shape/style/material reacts the least to external sound pollution. Once the few best candidates have been identified, develop test to experimentally determine optimal mount configuration and material.
Part 1: FEM
Process: Simultaneously build models in SolidWorks and design analysis in COMSOL
Solidworks: base/shaft models based on measurements of actual optical mounts, optic holder models matching mass and moment of inertia from downloaded models from their manufacturers. ~Progress: base/shaft done, optic holders almost done.
COMSOL: Design analysis that first does a general eigenfrequency analysis to find general vicinity of modes, then does a full modal frequency analysis. ~Progress: finished, ready for models to be imported. 
78

Mon Nov 18 15:56:59 2013 
Chris Couste  Optics  Analysis  Representative Models  The simpler models of the optical mounts are finished, they will be run through the comsol analysis software soon.
see pictures below:

Attachment 1: Custom_Closed.JPG


Attachment 2: custom_open.JPG


Attachment 3: Stock_Closed.JPG


Attachment 4: Stock_Open.JPG


79

Fri Nov 22 21:04:53 2013 
Chris Couste  Optics  Analysis  analysis  The analysis is making nice eigenmode and stress mode models, but the displacement experiment needs work. Should be fixed by monday. 
80

Fri Jan 24 17:26:38 2014 
Chris Cousté  Optics  Analysis  Mount Analysis Functional!  The comsol eigenmode analysis is complete, and the only thing left to do on this part of the project is to run the analysis on a range of different configurations of optical mounts as well as a range of materials. This compilation will be posted on this elog in the next few weeks, due to the fairly long runtime of the analysis software. 
81

Sun Feb 16 21:10:56 2014 
Chris Cousté  Optics  Analysis  Optical Mount data compilation 1: Aluminum  The time is finally here! this is the highest displacement of each mount in its lowest few eigenfrequencies, using 6061 Aluminum as a material. pictures will be added in a future log, because I'm going to make them into one file. Other materials will also be tested to see if there is variance in these findings, but only relevant data will be posted.
tl;dr findings: thickopen is the best combination because it has less displacement in eigenmodes and its eigenmodes are all very high/spread apart in frequency.
long:
(key) Eigenfrequency (Hz), highest displacement (mm), direction
Thin base, open mount:
315, 8.07, leaning back (laser will angle up)
316, 7.94, leaning to the side (laser won't change direction if the optic is flat, i.e. point of contact moves along surface of optic)
924, 10.50, twist to the side (laser will angle to the side)
1949, 8.89, twist about optic, slight lean forward (point of contact will move slightly, laser may angle down)
2144, 10.47, another leaning forward mode
thin base, open mount:
322, 13.29, leaning forward and translating to the side
328, 13.36, same as above
994, 17.66, twisting to the side (laser will angle to the side)
1943, 13.86, intense leaning forward
2100, 14.27, twist about optic
thick base, closed mount:
1329, 13.23, leaning back
1347, 13.32, leaning back and to the side
3577, 15.46, twisting to the side
thick base, open mount:
1326, 9.69, leaning to the side
1335, 9.77, leaning back
3504, 13.88, twisting to the side

All of these frequencies were checked for all four models, and the max displacements in those cases ranged from 100 picometers to 10 femtometers, so it's pretty reasonable to base decisions off of the displacements at the eigenmodes.
Based on this, it seems that for both thick and thin, open is the best, considering it displaces significantly less than closed (which is reasonable because there is less mass at the top to be thrown around). The choice between thick and thin depends on what frequencies we believe the mounts will be subjected to. If we are probably not going to have much sound at frequencies above 1000Hz, then thick is the better option (which is nice because it is the stock stalk instead of the custom made one). However, if there will be plenty of high noise, the thick is still probably the better option because it has fewer eigenmodes in the same range of frequencies.
here's a sample of two pictures of the relevant eigenmodes of the winner.

Attachment 1: aluminum2.png


Attachment 2: aluminum3.png


105

Tue Jun 28 15:50:45 2016 
Joy Westland  Mechanics  Analysis  A Simple Model of the Modal Analysis of a Cantilever Circular/Cylinder Cross Section in ANSYS tutorial  Here is a tutorial to implement a simple Modal Analysis of a Cantilever Cylinder
1. Open the ANSYS workbench
2. Drop and drag the "Modal" analysis system into the project schematic
3. Right click on "Engineering Data" and edit the material. There are predetermined properties found in the "Engineering Data Sources" (Right click on the description or press the books in the top right hand corner). For this tutorial, fused silica was created with these parameters:
 Young's Modulus E: 7.2E10 Nm^2 (Isotropic Elasticity)
 Density p: 2200 kgm^3 (Density)
 Specific Heat c: 770 J/kgK (Specific Heat)
 Isotropic Thermal Conductivity k: 1.38 W/mk (Isotropic Thermal Conductivity)
 Coefficient of Thermal Expansion: 3.9E7 (Isotropic Secant Coefficient of Thermal Expansion)
 Reference Temperature: 21 C (Isotropic Secant Coefficient of Thermal Expansion)
 Poisson's Ratio: 0.17 (Isotropic Elasticity)
In order to put in the parameters, just drag and drop from the toolbox. Also make sure the filter is turned off (that's the filter symbol next to the engineering data sources/books in the top right hand corner) this will allow a user to few all the parameters.
4. Return to the project, and open/edit the geometry. The cantilever that was used for this tutorial is a cylinder with a diameter of 2 mm with a length of 200 mm, shown in image 1.
 How to Create a geometry:
 Open up the geometry
 Decide which plane you want to start drawing in: XY Plane, ZX Plane, YZ Plane
 Click on "new sketch" in the 4th row from the top after "Sketch1" it's blue sketch
 Whatever plane you decide to choose, press the "Look at Face/Plane/Sketch" icon in the top right hand corner (The person looking at a face)
 This will allow the view to go to the view where you want to sketch on
 Click on "Sketching" which is next to the "Modeling" right above the "Details View"
 Here will be Several Options
 Draw: This is where you can decide what shape you want to use
 Modify: This is the extras such as fillets, chamfers and corners
 Dimensions: This is where you can select what type of dimensions you will be using
 Constraints: This is where you can constrain lines or connect lines
 When you want to connect lines together, click on each end and make them "Coincident" this way the shape will connect when you change dimensions
 Settings: This is where you can show a grid or snap grid
 Sketch Color:
 Teal: Underdefined: The sketch does not have enough dimensions to constrain and tell the program what exactly the shape is
 Blue: Defined: The sketch is fully defined and does not need any more dimensions
 Red: Overdefined: The sketch has repeating dimensions that causes errors
 For this tutorial, draw a circle. I made mine coincident to the middle by pressing the vertex in the middle and one of the axis and then repeat the process with the second axis to center the circle
 Dimensions: The diameter is 2 mm, under dimensions to display the number, click "Display" and then "Value" instead of "name"
 To change the dimensions, the "Details View" in the bottom left hand corner will say D1 and from there you can change the diameter
 From there in order to create the length of the cylinder, press "Extrude" at the top (4th row) towards the right. From there it will ask what geometry, click on your sketch and click "Apply". Make sure the "Operation" is "Add Material" to Extrude.
 Then choose a "Depth" in this case it's 200 mm. Then click "Generate" the lightning bolt next to the "New Sketch" symbol.
 Your cylinder should look like image 1.
5. Return to the project and open the Model tab. With that open, go to the subtab of "Geometry", named "Solid". Under Solid go to "Assignment" under "Material" to assign Fused Silica or whatever material you decide to choose.
6. Go to the "Mesh" tab. In this case the "Element Size" under the "Sizing" tab is 0.04 and the "Relevance Center" as "Fine". Right click on the mesh and press "Generate mesh".
7. Under the Modal tab, go to the Analysis Settings and change the "Max Modes to Find" for ANSYS to calculate. In this tutorial, the amount of modes that were used was 17.
8. Right Click the Modal and press on "Fixed Support". This will make the bar a cantilever bar, once the geometry is set to fix one face of the bar. Press on one of the faces and click on "apply" on the Geometry.
9. Right click on the solution and press "Solve". Let ANSYS run the modes through.
10. On the Solution tab, the "Tabular Data" is listed but the Total Deformation has not been listed. In order to do so, select all the frequency of the modes in the column and right click and press "Create Mode Shape Results".
11. Once loaded, the total deformation will have a lightning bolt next to each entry, right click and press "Solve" or "Evaluate All Results". This will make all the entries have a green check mark.
12. This will give the user the ability to animate each entry. This is down by clicking the play (sideways triangle button). ANSYS will run through the simulation for that mode that is selected, shown in image 2.
13. Using Mathematica (or another computational system) input the analytical solution for a Cantilever bar fixed to one end.
 w_{n}=(k_{n}L)^{2}Sqrt(EI/mL^{4})
 Where:
 w_{n} is the frequency measured in radians/second
 E = Young's Modulus
 I = moment of cross section
 I_{rectangle }= ba^{3}/12 (b and a are the sides)
 I_{circle} = (Pi/64)(d^{4}) (d is the diameter)
 m = mass per unit length
 L = length of the bar
 k_{n} relates to the amount of nodes
 k_{1}L = 1.875, k_{2}L = 4.694, k_{3}L = 7.855, k_{4}L = 10.996, k_{5}L = 14.137
 n greater than 5: k_{n}L = (2n 1)(Pi/2)
 Convert w_{n} to f measured in Hz
14. Compare only to the bar actually bending, not twisting or contracting. There are modes that are the same due to the symmetry of the bar. In image 3 the underlined frequencies compare to the analytical calculations (Mathematica) and the computational calculations (ANSYS). 
Attachment 1: Simple_Cantilever_200_mm_Long_Geometry.PNG


Attachment 2: Simple_Cantilever_200_mm_Long.PNG


Attachment 3: Comparison_between_ANSYS_and_Calculations_Fused_Silica_200_mm_L.PNG


106

Tue Jun 28 16:54:23 2016 
Nikhil Mathur  Mechanics  Analysis  Ansys 14.5 Indroductory Tutorial: Modal Frequency Convergence  This tutorial will go through how to show frequency convergence for a cylindrical cantilever using Ansys 14.5 and Mathematica. 
Attachment 1: ConvergenceTutorial.pdf.gz

107

Wed Jun 29 14:14:44 2016 
Joy Westland  Mechanics  Analysis  ANSYS Tutorials with Basic Meshing  Here are a series of tutorials for basic meshing principles from ANSYS Meshing Basics:
 ANSYS Meshing Fine Mesh Basic Tutorial 1: https://www.youtube.com/watch?v=sZIX3CJkWBE
 ANSYS Meshing Method Basic Tutorial 2: https://www.youtube.com/watch?v=ATSnvFUbEk4
 ANSYS Meshing Refinement Mesh Basic Tutorial 3: https://www.youtube.com/watch?v=QVfJiydCu2g
 ANSYS Meshing Sizing Mesh Basic Tutorial 4: https://www.youtube.com/watch?v=ODJ11YDC6ac
 ANSYS Meshing Spheres of Influence Basic Tutorial 5: https://www.youtube.com/watch?v=AQlGccK_X0
 Body of Influnce 6: https://www.youtube.com/watch?v=3E1p1w32jt0
 Section Plane 7: https://www.youtube.com/watch?v=hhZOJSESJEU
 Mapped and Match Control 8: https://www.youtube.com/watch?v=3pogywR8Vgw
 MultiZone and Inflation 9: https://www.youtube.com/watch?v=vxVBvANHuB0
 Sizing Soft/Hard 10: https://www.youtube.com/watch?v=D8UzdC5Gmhs

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:

143

Fri Jul 21 16:42:57 2023 
Rajashik  Optics  Analysis   Initially, ideal conditions were considered, without any deformations or thermal heating effects in the mirrors. This offered a baseline model from which subsequent tests could draw comparisons.
Following this, more complex scenarios were introduced into the model: mirrorshape defects caused by suspension and thermal lensing due to heat absorption. These heat absorption rates were set at 0.5 ppm in the ETM coating and 0.3 ppm in the ITM coating.
The lensing calculations for these scenarios were performed using the power observed in the nondeformed mirror case. This lensing needs to be made dynamic in future codes to get the modeling correct. Nonetheless, introducing these conditions resulted in the emergence of higherorder modes (HOMs). Unfortunately, there seems to be a loss of power in the cavity that I haven't been able to resolve yet. Notably, this power loss persisted across the range of mirror position scans showing that this is not due to just a change in resonance length.
Also, the Hello Vinet Function in Finesse 3 seems to yield a higher degree of lensing than typically reported in the literature. I need to check with someone about what assumption I am taking incorrectly here. It's plausible that the heightened lensing effect is contributing to the observed power loss in the cavity, possibly due to high curvatures.
Looking forward, my focus will be on determining the source of this lensing discrepancy first. This analysis will involve a more detailed examination of the variables contributing to thermal lensing and potential causes for the higherthanexpected lensing observed in the simulation. Ultimately, the aim is to create a dynamic deformation model that can more accurately represent these phenomena within the aLIGO setup before introducing ring heaters and controls. 
Attachment 1: Cavity_Scan_Before_Adding_Defects.png


Attachment 2: Transmitted_Power_Before_Adding_Defects.png


Attachment 3: Cavity_Scan_After_Adding_Defects.png


Attachment 4: Transmitted_Power_After_Adding_Defects.png


82

Thu Apr 17 16:52:09 2014 
Evan  General  Characterization  Interpreting logfiles and picking a solver  Here are two entries by Walter Frei on the Comsol blog that I've found useful.
Solving Linear Static Finite Element Models: Tells you how to interpret all those numbers that Comsol dumps into its logfiles.
Solutions to Linear Systems of Equations: Direct and Iterative Solvers: Explains a bit more about Comsol's solvers. Apparently, PARDISO is usually fast and SPOOLES is usually slow (the tradeoff is apparently that SPOOLES uses less memory). I've found this to be personally true with at least one model I've worked with. 
4

Thu May 2 14:00:36 2013 
Koji  General  Configuration  test mass TR with Levin's approach  Thermorefractive noise in a finite cylindrical/infinite test mass with Levin's approach
Location of the codes: 40m SVN repository
comsol/thermorefractive/
This code realizes Levin's calculation on thermorefractive noise
doi:10.1016/j.physleta.2007.11.007
and duplicates the result of D. Heinerts paper
DOI: 10.1103/PhysRevD.84.062001
Also the result is compared with Braginsky's result in 2004.
doi:10.1016/S03759601(03)004730
 The code applies gaussianshaped heat into a cylindrical mirror.
 The heating/cooling is sinusoidal and the dissipation (heat flow) is calculated in COMSOL.
 The time series result was analyzed in MATLAB to extract the single coefficient corresponds to the transfer function.
This way the effect of the initial transient was avoided.
 Unfortunately direct measurement of frequency response in COMSOL was not available as the heat flow is not modal.
If we make a fourier analysis of the partial differential equation and solve it in COMSOL using arbitrary PDE solver,
we may turn this time dependent analysis into static analysis.
All of the calculation was driven from MATLAB. So you have to launch "COMSOL with MATLAB". 
Attachment 1: thermo_refractive_1D_axisym_result.pdf


76

Mon Sep 16 10:07:32 2013 
Evan  General  Configuration  COMSOL 4.3 on the OS X command line  If you're running Matlab scripts that iterate over simulation parameters (à la Tara), you might find it useful to be able to run everything from the command line (instead of launching the "Comsol with Matlab" GUI application).
First, the comsol command wasn't in my path, so I symlinked it to someplace where bash could find it:
ln s /Applications/COMSOL43a/bin/comsol /opt/local/bin/comsol
I then started Comsol/Matlab using a slight modification of the commands given in the Comsol 4.3 help file "Running A COMSOL Mfile In Batch Mode Without Display":
comsol server > server_log.txt &
matlab nodesktop nosplash r "run /Applications/COMSOL43a/mli/mphstart; comsol_script; exit"
The first command starts the Comsol server in the background and cats its output into a log file. The second command launches Matlab, runs the initialization script which makes Matlab aware of Comsol, runs my Comsol/Matlab script (comsol_script.m), and then exits. 
99

Tue Sep 30 11:30:27 2014 
Nic, Dmass, Evan  General  Configuration  Gravity in Comsol  Here is a set of slides by Yoichi Aso on how to handle gravity in Comsol. 
100

Sat Sep 5 15:04:43 2015 
Dennis Coyne  Mechanics  Configuration  summary of FEA modal model to State Space model  At the 2014 commissioning workshop, I presented a summary of my efforts in converting finite element modal models into state space models:
https://dcc.ligo.org/LIGOG1400099
I also provided the attached written summary for a report on the workshop that Aidan Brooks and Gabriele Vajente are preparing. 
Attachment 1: FEM2SS_summary.docx

101

Sat Sep 5 15:17:41 2015 
rana  General  Configuration  FEA logs merged  I moved the only entry from the 'ENG_FEA' log into the COMSOL log and then renamed that logbook as 'FEA' since we don't need two FEA logs.
Also renamed 'AdhikariLab' log as ATF. 
102

Fri Dec 4 18:32:02 2015 
rana  General  Configuration  COMSOL: remote server w/ matlab  This summarizes how to get the remote comsol server to run. COMSOL 5.1.0.234 is now on tegmeni thanks to Larry.
On the server:
rana@tegmeni~> /usr/local/comsol51/multiphysics/bin/comsol server login force
This starts up a comsol server instance, listening on port 2036. 'login force' will ask you to supply a username and password which you make up. You will have to enter these later from the client side.
On my laptop, from the MATLAB prompt:
>> addpath('/Applications/COMSOL51/Multiphysics/mli/')
>> mphstart('tegmeni.ligo.caltech.edu', 2036,'uname','pword')
That's it! Now you just run the matlab script which runs the COMSOL file.
I'm attaching a tarball of the .mph file (written by Dmitry Kopstov from MSU) and the associated matlab scripts which change the parameters, as well as looping over test mass thickness to produce a plot of Brownian noise PSD v. 
Attachment 1: BrownianSweep.tgz

110

Mon Jul 24 15:35:34 2017 
Mariia  General  Configuration  Running Comsol to Matlab  WIth Gautam's help, I have created a user directory in 40 meter Lab and copied Rana's documents (MATLAB coating files) from flash card into it. After that, from this elog by Rana : COMSOL: remote server w/ matlab from Fri Dec 4 18:32:02 2015, ran the matlab document BarrelCoating which resulted in the following error
Messages:
Could not obtain license for#LiveLink for MATLAB
Could not obtain license for LiveLink for MATLAB.
License error: 4.
Licensed number of users already reached.
Has anybody seen this before?

111

Mon Jul 24 15:54:26 2017 
Koji  General  Configuration  Running Comsol to Matlab  The number of licenses already used by whom / still remains can be confirmed by running the following command on a comsolinstalled linux machine
$ cd /usr/local/comsol51/multiphysics/license/glnxa64
$ ./lmstat c ../license.dat a

112

Wed Jul 26 20:14:46 2017 
rana  General  Configuration  Running Comsol to Matlab  I've just tried this out on my desktop machine using COMSOL 5.1 and its still working. Which COMSOL is installed on optimus at the 40m ? 
113

Fri Jul 28 15:48:58 2017 
Mariia  General  Configuration  Comsol batch for windows  Using the written path from elog by ericq: Computer Scripts/Programs, Comsol can be run from the directory on the distant computer: /cvs/cds/caltech/apps/linux64/comsol51/bin/glnxa64/comsol batch inputfile Model1.mph outputfile Model_solv.mph. To transfer files from Linux to Windows : the command pscp.
The method:
1. To download PuTTY and make a coonection with a distant computer.
2. Linux terminal: ssh optimus % allows to go into the whole system
3. Windows terminal: pscp "the path of the file" controls@nodus.ligo.caltech.edu:/users/.../ % allows to transfer mphfile from Windows to Linux
4. Linux terminal: /cvs/cds/caltech/apps/linux64/comsol51/bin/glnxa64/comsol batch inputfile Model1.mph outputfile Model_solv.mph % calculation on the distant computer, without outputfile solution would be stored in the input file
5. Windows terminal: pscp controls@nodus.ligo.caltech.edu:/users/.../Model_solv.mph Documents\ % copy file from Linux to Windows, Documents is the name of the folder in Windows

139

Tue Aug 11 11:16:29 2020 
aaron  General  Configuration  COMSOL with Matlab without display  When running comsol with matlab interface on sandbox1, it is usually most convenient to ssh with screen forwarding (eg 'CY') and launch COMSOL with matlab by following the instructions in the livelink manual.
Sometimes, it is necessary to run COMSOL without any display available. In that case, the Instructions in the manual are a little unclear. Here are the detailed steps that let me run my script '/home/amarkowi/metamaterials/run_spiral.m' with no screen forwarding.
1. ssh onto sandbox1 by entering the following at your laptop command prompt
password: [enter sandbox1 credentials]
2. start a tmux session for starting the COMSOL server
aaron.markowitz@sandbox1:~$ tmux
3. start a comsol server from the tmux session
aaron.markowitz@sandbox1:~$ comsol54 mphserver login force port 2020
Username: whoever
Password: whatever
Confirm password: whatever
4. detach the tmux session by pressing ‘ctrlb’ followed by ‘d’
5. (optionally, you can start a new tmux session for your matlab work by running tmux again at your main sandbox prompt)
6. Start matlab by running
>> matlab nodesktop mlnosplash
6. Add the comsol directory to the matlab path by running at the matlab prompt
>> addpath(‘/localhome/comsol54/multiphysics/mli/‘)
6. Start the matlab with comsol interface by running the following at the matlab prompt
7. change into the directory containing the script, and run it
>> cd /home/amarkowi/metamaterials
>> run_spiral

140

Tue Aug 11 16:35:07 2020 
aaron  General  Configuration  COMSOL: remote server w/ matlab  To run COMSOL on sandbox1 with no graphical Interface, here are the steps that worked for me (Tue Aug 11 16:35:51 2020) from a mac on the Caltech VPN.
1. ssh onto sandbox1 with screen forwarding (Y). Make sure you have a compatible version of XQuartz or a substitute. C specifies data compression, which may be useful on slow connections
password: [enter sandbox1 credentials]
2. Launch Matlab with COMSOL 5.4, specifying no graphical interface, by running
aaron.markowitz comsol54 mphserver matlab nodesktop mlnosplash script_name’
3. If there is still a splash screen from the COMSOL server, you will have to specify nosplash by adding the following line to your .bash_profile (in your home directory)
export COMSOL_MATLAB_INIT=’matlab nosplash’
4. You can run whatever comsol script you need. Make sure that in your script you import the comsol functions by calling the following
import com.comsol.model.*
import com.comsol.model.util.*

1

Thu Dec 13 18:47:18 2012 
rana  General  General  COMSOL: who's using up all the licenses?  When you want to find out who's using up all the licenses, you can run lmstat a do find out.
Specifically, with Mac COMSOL, I do:
> cd /Applications/COMSOL43a/license
> maci64/lmstat a c license.dat
Users of COMSOL: (Total of 5 licenses issued; Total of 1 license in use)
"COMSOL" v4.3, vendor: LMCOMSOL
floating license
dmassey c22042.local (v4.3) (ancha/1718 1074), start Thu 12/13 18:26
etc.

2

Tue Feb 5 13:15:13 2013 
Dmass  General  General  COMSOL: who's using up all the licenses? 
Quote: 
When you want to find out who's using up all the licenses, you can run lmstat a do find out.
Specifically, with Mac COMSOL, I do:
> cd /Applications/COMSOL43a/license
> maci64/lmstat a c license.dat
Users of COMSOL: (Total of 5 licenses issued; Total of 1 license in use)
"COMSOL" v4.3, vendor: LMCOMSOL
floating license
dmassey c22042.local (v4.3) (ancha/1718 1074), start Thu 12/13 18:26
etc.

Whenever I try to open a model (.mph file) which I built in COMSOL and has never been touched by any CAD software, it uses the CADIMPORT license, and we only have 2 of these. I can open a different mph file without the CADIMPORT being used.
I don't yet understand what makes something you build in COMSOL and save in COMSOL try to use the CADIMPORT module.
CADIMPORT licenses currently in use: (2/2)
ligo m5 m5 (v4.3) (ancha/1718 202), start Wed 1/30 11:39
mattabe Abernathydesktop /dev/tty (v4.3) (ancha/1718 1016), start Mon 2/4 19:52 
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


12

Wed Jun 5 20:54:59 2013 
Deep Chatterjee  General  General  Bessel Function roots  During the process of evaluating the PSD from Sec. V of Liu and Thorn, I chanced to write a simpler root finding algorithm applying bisection to find the roots of J1(x).
The difference between this and the algorithm by
Greg von Winckel goes as
The difference is of the order 10^{7} and can be reduced by reducing the tolerance. However, It should however be noted that bisection is a crude algorithm for rough usage and differences become pronounced for larger n.

Attachment 2: bessel_zeros.m

function X = bessel_zeros(k, n)
%BESSEL_ZEROS : calculates the first k zeros of the bessel function Jn
% k : No. of ROOTS to be evaluated
% n : Order of the Bessel function Jn
% X : stores the roots in serial order i.e. X(1) gives the first root,
% X(2) gives the second and so on
X = zeros(1,k);% empty array of lenght k
count = 1;%this acts as a counter to k
dx = 0.1;%step size within which, by assumption, no roots exist
x = 0;
... 22 more lines ...

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

15

Fri Jun 7 17:41:41 2013 
Deep Chatterjee  General  General  Disscussion of the code of a 'comsol with matlab' model file  In this post the matlab code to build a model using COMSOL with matlab is analyzed from Koji Arai's codes on the calculation of TR noise in the 1D axisymmetric case. The post primarily describes the various commands used to interface with COMSOL.
It should be noted that the matlab script analyzed here is not the master script that will perform the simulation. The "main" file is called the "Main_thermo_refracive.m" found in the SVN repository TR noise_comsol
>The main script defines a structure called 'param' which stores all the parameters including the material properties and those related to running the simulation in comsol. This structure 'param' is hence passed on as arguments to all the other functions (in the SVN rep)
that compute the analytical solution or the COMSOL simulations.
>The part of the code ( the second if...end block) in "Main_thermo_refracive.m" that calls the COMSOL simulations, calls the function "thermo_refractive_COMSOL_1D_axisymmetric".
>In this function, the parameters are extracted from the structure 'param' and stored in separate variables. The time step in the simulation, the end time, mirror radius are defined and the function
that generates the model "thermo_refractive_COMSOL_1D_axisymmetric_model" is called which prepares the model in COMSOL.
> it should be noted that it is much of much ease to the user to define a simple model in COMSOL and export it in Matlab and then analyzed the same.
I had created a simple such file which a simple 3D metal bar is created and added some physics and study to see how the matlab file.
I have attached the same which also describes the relevant packages and subroutines in the comments. However, I had not erased the COMSOL history before the export hence it has significant amount of superfluous code. One may however have a look at the same

>The first thing done by the model builiding function is that it extracts a few parameters and stores them in separate variables in the first lines
f_mod = param.COMSOL.simulation_modulation_freq;
dt = param.COMSOL.simulation_time_step;
t_end = param.COMSOL.simulation_end_time;
>Since COMSOL accepts string arguments, the time steps are defined as a string
trange = ['range(0,' num2str(dt) ',' num2str(t_end) ')'];
**There is however a query regarding the keyword range which in Matlab returns the difference between the maximum and minimum of the list passed to it. range() passed to COMSOL would probably create an array fro 0  t_end in steps of dt.
The above line converts the quantity 'dt' and 't_end' to a string and concatenates it along with the rest of the string to make a valid string that would be understood by COMSOL.
>The COMSOL class is called next by the following two import statements (after a series of displays on the screen showing the status of the simulation).
>Managing models( like creation and destruction) is handled by the modelUtil class in COMSOL with Matlab and hence to create a model modelUtil.create() is called. "model.remove()" is used to destroy one such object. Line 29 says 
model = ModelUtil.create('Model');
>The above statement created an object called 'model' in Matlab and and a model in COMSOL by the name of 'Model'. The various attributes of this object is defined by the next few statements. Like,
[line 31]  model.name('thermo_refractive_COMSOL_1D_axisymmetric_model.mph');
The above statement assigns the filename of this COMSOL model created. Note that, since in Matlab, the Matlab object 'model' is used to invoke the functions i.e. 'model.name()' and not 'Model.name()'.
> The 'model.param' contains all functions related to setting and describing the parameters in the model.
*Note that the 'param' following keyword 'model' connected by the '.' has nothing to do with the 'param' structure used in the script "Main_thermo_refracive.m" which have been named anything else
> Now the model.param.set(<P>,<expr>) is used to give the parameter P an expression expr both being string.
The model.param.descr(<P>,<des>) has the same format but gives the description of the parameter P as des (something that is understandable in common terms).
This can be understood in the statements in [line 46  51]
model.param.set('beam_shape', '2/(pi*beam_size^2)*exp(2*r^2/beam_size^2)');
model.param.descr('beam_shape', '');
model.param.set('beam_intensity', 'beam_power*beam_shape');
model.param.descr('beam_intensity', '');
model.param.set('dT', 'mod1.TT_amb');
model.param.descr('dT', '');
An empty quote in the description implies that no description is given. The above statements define the parameters mentioned alongside them according to Heinerts paper.
> [line 5361]. the string 'var1' is used to tag all the global variables that are created. variables are expressions created out of the parameters.
Just like the above case, 'set()' is used to create variables and give them an expression, here as one can see 4 variables are created T_amb, beam_size, beam_power. and f_mod.
model.variable.create('var1');
model.variable('var1').set('T_amb', [num2str(param.material.temperature) '[K]']);
model.variable('var1').descr('T_amb', '');
model.variable('var1').set('beam_size', [num2str(param.beam.radius) '[m]']);
model.variable('var1').descr('beam_size', '');
model.variable('var1').set('beam_power', [num2str(param.COMSOL.beam_power) '[W]']);
model.variable('var1').descr('beam_power', '');
model.variable('var1').set('fmod', num2str(f_mod));
model.variable('var1').descr('fmod', '');
> The next section deals with the geometry. In this case the model is a 1D axisymmetric model. Hence one has the the digit 1 which specifies the dimension and axisymmetric is set to 'true'
model.geom.create('geom1', 1);
model.geom('geom1').axisymmetric(true);
> The run() function is to build the geometry
model.geom('geom1').run;
> The 'Interval' feature is present in the case of 1D models. The following line creates an interval feature called 'i1'
model.geom('geom1').feature.create('i1', 'Interval');
model.geom('geom1').feature('i1').set('p2', num2str(param.mirror.radius));
The above line is used to set the value of the right endpoint as the radius of the mirror. The string 'p2' is used to denote the right end point..
So this step creates a cylinder of radius equal to the radius of the mirror. In case of an infinite mirror an extra interval is added from the the radius to twice the radius.
It would help if more description on "interval" is left as comments and on why the extra interval was added and how does it make the calculation different.
> The interval is run using the 'run()' command.
> The material definition is added in the following section. Probably the data seems to have been taken from an external file rather than the COMSOL material library. Details of the file would be of help from the author.
The following lines cannot be understood. Comsol help on"propertyGroup" gives no results found
model.material('mat1').propertyGroup('def').set('heatcapacity', [num2str(param.material.specific_heat_per_volume) '[J/(kg*K)]']);
model.material('mat1').propertyGroup('def').set('density', [num2str(param.material.density) '[kg/m^3]']);
model.material('mat1').propertyGroup('def').set('thermalconductivity', [num2str(param.material.thermal_conductivity) '[W/(m*K)]']);
> The meshing is done on auto.
> As far as the Physics is concerned i.e. what kind of results we expect out of this geometry, the 'HeatTransfer' is selected.
A 1D heat source is created called 'hs1'
model.physics('ht').feature.create('hs1', 'HeatSource', 1);
and the Heat Source is time varying and given by beam_intensity*sin(2*fmod*pi*t). This can be seen in lines 121122
model.physics('ht').feature('hs1').selection.set([1]);
model.physics('ht').feature('hs1').set('Q', 1, 'beam_intensity*sin(2*fmod*pi*t)');
> If the mirror is infinite, something called 'InfiniteElements' is applied to the outer interval. More details on "infiniteElements" and why was it used would be helpful.
> A study involving the feature "Transient" was used. However, it is required to know what study was implemented from the GUI since the term "Transient" was not found under Study Steps in the GUI
> The time step tweaking of COMSOL is stopped by the following line 166
model.sol('sol1').feature('t1').set('tstepsbdf', 'strict');

Attachment 1: Iron_Bar_test1_model.m

function out = model
%
% Iron_Bar_test1_model.m
%
% Model exported on Jun 10 2013, 10:13 by COMSOL 4.3.0.151.
% This is an exported file from COMSOL in the '.m' format. I am adding
% comments to the lines.
% As one goes along the steps one followed in COMSOL and the corrsponding
% Matlab file, the purpose of the various subroutines become clear. The
... 430 more lines ...

