Here, I describe some test cases to see if COMSOL's solutions are agreeing with some simple analytical solutions.  Right now, I have two plots showing COMSOL's solution and my analytical solution on separate plots.  I will be plotting there difference to see if they really match up.

 It does appear that the simplified model is only relevant for the simulations.  To quote Heinert: "An efficient computation is only possible for the simple model, as the advanced model would require an element of size more than 10⁶ ."  I have run Koji's code that replicates Heinert's figure 3.  I have attached the resulting temperature distribution and noise amplitude curve.  In the noise amplitude curve, the red line is the analytical result, while the dots are from COMSOL.

The next step is to convert this code to an efficient complex time-independent solution.  As stated before, my main concern here is whether COMSOL actually solves the right equation in the stationary case.

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 re-edited your comments to make them more easily read.

I used http://www.sciweavers.org/free-online-latex-equation-editor to convert your latex to png files for easy reading.

At this point, my are goals are to:

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 zeroth-order 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."

At this point, my are goals are to 1) convert the time-dependent 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 time-dependent 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 steady-state 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 zeroth-order 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.

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.

 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: thick-open 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.

 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.

 The analysis is making nice eigenmode and stress mode models, but the displacement experiment needs work. Should be fixed by monday.

The simpler models of the optical mounts are finished, they will be run through the comsol analysis software soon.

see pictures below:

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.

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:

I then started Comsol/Matlab using a slight modification of the commands given in the Comsol 4.3 help file "Running A COMSOL M-file In Batch Mode Without Display":

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.

For the material parameters of Sapphire at 300K, the TE and TR Noise profiles, though not very close, lie close to within an order of magnitude. Sapphire has a positive coefficient of linear expansion. We just inverted the sign of this quantity
and the ran the codes that puts heat and pressure simultaneously to the test mass. The total noise looks to be lower than the TR noise which is greater.

Mirror radius 0.25[m]

Mirror height 0.46[m]

Beam Radius 0.09[m]

If we have physical parameters which make the two Noise sources come closer to each other and then flip the sign of alpha, we may be able to see some noise reduction to a greater extent.

After the simulations have been found to match to a fair extent with the analytic results by Heinert and Liu and Thorne, the attempt is check out the parameters for which the TE and TR noise are
close to each other. This was done with the analytic results. The frequency range 10 - 1000 Hz was looked into. In between this range, the quantity that was minimized was the absolite value
of the logarithm of the ratio between the TR and TE noise. The fminsearch function was used to minimize the mentioned quantity. The parameters that were changed were - conductivity, thermo
optic coefficient and coefficient of linear expansion. The reason for choosing these three were -

> TR noise is independent of coefficient of linear expansion

> TE noise is independent of thermo optic coefficient

> The power law dependence on conductivity is different for the TE and TR cases as can be seen from the analytic expressions

once the code returned the optimized parameters, these values were plugged in and the results were plotted.

**Note that the minimization was done for frequencies between 10 to 1000 Hz

An error was being encountered in the computation of the TR noise lately. It was observed that while running the simulations in the case of the materials which have a lower value of the thermal diffusivity (silicon / sapphire at room temperature), the simulated result were slightly off from the analytic result. On the other hand, if the simulation was run with a material of higher diffusivity(same materials at lower temperature), the match would be better. The reason being the transient solution not dying off significantly during the period of the simulation. Since a time average was being taken of the quantity integral{grad_T ^2},  the transient contributed to the integral. To get the correct value, the fourier coefficient of the time signal of integral{grad_T ^2},  was extracted at twice the frequency of the pressure oscillation. The reason being that the signal was squared. Extracting the response at this frequency after the integration is logical since the integration is over space while the response we extract is over time.

The same procedure was also applied to the TE noise calculation. However, this time we obtained similar result as the case where this procedure was not applied but a simple time averaging was performed. The tail of the plot, at high frequency, is still seen to deviate from the analytic result of Liu and Thorne as was the case previously. A plot is attached showing the spectrum for Fused silica at 290K. The conclusion being that the transients do not affect the TE noise calculation - the plot stayed the same even after filtering them out. This is probably because unlike the TR case which has a heat source present along the cylinder axis, the TE noise calculation involves applying pressure only on the face of the cylinder, and the transient do not contribute much to the volume integral.

The issue related to the difference between the analytic and simulated values has been resolved. The codes seems to give reasonable match
between the analytic and simulated case. There is, however, a difference between the formulas being used from the previous cases. Note that
the 1/2 in front of Eq.(15) of Heinert is a because the time average has already been considered. However, in the present codes, the volume
integral of grad_T is evaluated in COMSOL and exported as a function of time. It is then averaged in MATLAB. Thus the factor of 1/2 is to be
omitted in this case(see Liu and Thorne Eq.(5). The presence of this extra factor of 1/2 was giving error in the last upoaded plots. From the
relative difference plot, one can see the maximum difference between COMSOL and analytic results go upto 7% but for most of the graph
it is close to 1% which is a fair result.

The discrepancy related to the difference between the analytic and COMSOL results has been partially addressed. Attached is another
plot showing the comparison. The ratio this time between the COMSOL results and the analytic results is between 0.7 - 0.8. This difference
will be looked into. It is, however, observed that the difference is not a constant factor - it has to do with the model file.

 PDF please instead of EPS or BMP or JFIF or TARGA or GIF or ascii art.

In this post I simulate the procedure of calculating the TR noise for finite cavities as proposed by Heinert and check for a
match.

The technique of performing all necessary calculations in COMSOL and exporting the results was applied to the TR codes.
It was seen that the codes gives similar output as the technique of extraction of Fourier coefficients in place of time averaging
as has been done in the codes of Koji Arai. One can see the output as the present code runs to be similar to the previous ones
found in the SVN.

However, the results in the present case were off by a constant factor close to 100. This maybe due to some 'm' - 'cm' or similar difference between
analytic calculations and COMSOL values of parameters. Although, it has not been found yet, the correction is hopeful to be
found soon.

The codes give results similar to the analytic result for other values of the mirror radius and beam radii (apart from the constant
factor I have mentioned above). One may have a look at the trend of the graphs between analytic and simulated values in the plots
attached. These plots are for the case when the mirror radius = 25m while the beam radius = 9 cm i.e. the original radii were 0.25m
and 9cm respectively i.e. the ratio has been changed by a order of 2.

As mentioned before the reason for the constant factor difference will be looked into.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

 Here is another plot with a mesh size slightly finer than the default Extra fine mesh in COMSOL.

One may notice that the value for the final frequency i.e. 10000Hz is different from the previous plot. 
It maybe that the error for the higher frequencies is a result of the FEA. However, it may also be that
the appropriate boundary conditions required for an infinite model break down at high frequencies.

In this post I report of the results of TE noise simulated by COMSOL for the TE noise of Infinite test masses.

The aim was to follow the procedure by Liu and Thorne in their analytic calculations so that the same model could be used for the other
geometries.

The simulation is done in a different way than the TR simulations. It was observed that the output given by COMSOL by the use of commands
like mphinterp() or taking an export resulted in certain discrepancies between the results computed in COMSOL and that read by MATLAB.

Thus, the volume integration of the temperature gradient is performed in COMSOL itself and the results of the integration for each time
are sent to files. Matlab read these values and time averages them to get the result as in the paper (Sec. 2 of Liu and Thorne).

The errors expected are
> Fourier analysis is not done at all. This would have involved exporting data which, as mentioned before is giving errors

> The numerical errors by COMSOL are therefore not filtered off.

> The plot differs from the analytic solution for larger frequencies over 3000 Hz.

> It is to be noted from the paper by Liu and Thorne that the TE noise for the finite and infinite case are not very different. In
fact the correction factor goes as O(1). Thus, differences between finite and inifinte cases are unlikely to be prominent
in the log scale plots

 The codes are put as a zip file. Corrections made to the codes will be uploaded as a reply.

 The dependence of the Power Spectral Density on the parameters F0, which is the pressure amplitude, has been somewhat corrected for in the codes attached. The mistake lay in the fact that
extraction of the Fourier coefficients was being done for the ht.gradTmag, the Fourier coefficients of which are not the correct ones to be used in this case since ht.gradTmag is not the correct time signal.
In the present case the Fourier coefficients are extracted from the r and z components, squared and added to carry out the calculation. However, the results are still off from the L&T results by order of 5.
The codes and the plots are attached.

The red values are the L&T result while the blue sketch is the simulation

 I ran sweeps over the ratioR range from 0.5 to 0.9 for test masses with modified height, and thus mass, which were otherwise equivalent to our previous test mass.  The attached plot was generated and indicates that the optimal ratioR value for the test mass was shifted significantly, for the half mass test mass into lower frequency regions, and for the 10x mass test mass into higher frequency ranges.  I took the data collected from the half mass run and fit a polynomial of degree 2 to it in Matlab using the command p=polyfit(x,y,2) and obtained p=1.089*10^-10*x^2 - 1.024*10^-10*x +1.774*10^-10 which has a minimum value for ratioR=0.47.  This fits with the plot generated as it appears to be near its minumum value, but I should run the simulation more around this value to see if it fits the expected behaviour.  I also did this fitting for the 10x mass curve, but based on the greater variance in the curve over this region of ratioR values I don't expect the predicted minimum to be as accurate, so if we want to find it we will have to search that region of ratioR values. For the 10x mass test mass, p=4.459*10^-10*x^2 - 8.859*10^-10*x + 6.289*10^-10 which has a minimum of about 1.

I will run simulations in these regions overnight, and edit this post with the plots/results.

 Following the previous reply regarding the use of mphinterp(), the appropriate changes in code were made. The overhead time due to the call to COMSOL was responsible for the longer
simulation time. Right now, for relatively large no. of total slices (about 27000) just the extraction of data happens within a few seconds. However, for the larger number of slices, the codes 
will probably still take some time to complete. The codes are attached.

 I am replying to report of the modification made to the script from the last one. The problem of too many radial and z slices have been avoided by using a slicing based on a gaussian
distribution both along the r and z directions - i.e. the slicing is heavy in the region around the surface is applied and much thinner towards the edges of the cylinder. The slicing has
been separately handled by a function named giveSlices() which does the slicing and returns the values as two arrays corresponding to r and z.
The extent to which the fineness of slicing is controlled in the geometry is by controlling the parameter SD in the code which is the standard deviation of the gaussian according to which
the slicing is controlled as being fine throughout or fine in the centre and thin with increasing r.

This method reduces the number of slices by almost 2 orders for significant large values of the dimensions of the test mass while it is expected to give similar results since slicing the
test mass finely at the edge is not required as the gradient of temperature almost falls to zero.

It is suspected that the long time of simulation is attributed to calling mphinterp() command a large number of times. In my next modification, it will be tried to use the command just once
followed by the proper data extraction from the output given by the command.

 We had noticed that when we ran an eigenfrequency analysis, there was often a node (see attached screenshot) on the face opposite the applied force.  This indicated that the noise might be lowered by reversing the face the force was applied to, but it also seems like this would be equivalent to increasing the ratioR value and making the primary face smaller while retaining the mass of the test mass.  It seems worthwhile to see how the Umax values compare for a beam incident on either face of a test mass in the optimal region with Umax=0.71.  For the initial configuration, Umax=1.47228*10^-10 J and changing the side of the mirror the force is applied to we get Umax=1.47148*10^-10 J, which is about the same value.

 I performed the same analysis we previously did on a test mass made of silicon (using the built-in silicon within COMSOL) and obtained the attached plots.  The plots show more of a difference between the two solutions optimal points than for fused silica, but they are still quite close.  All data is included in the attached tarball.  The values for ratioR=1 or a cylindrical test mass place the frequency of the lowest real eigenmode at 8491 Hz and the Umax value of  6.5484*10^-11 J.

The present aim is the calculation of the TE noise like Liu and Thorne. Although, the results of Liu and Thorne are for the case of an infinite mirror, in our case, we try
to model a test mass with large dimensions in an attempt to get closer to the result due to Liu and Thorne

Following my previous post - The presents results given by the COMSOL simulations gave the same profile as Liu and Thorne's result which has a frequency dependence
of 1/omega^2 but was displaced from the analytic result by a constant factor which was dependent on the applied pressure on the face of the test mass which should not be
the case.

Following the steps followed Liu and Thorne, we had constructed the test mass as a cylinder, large compared to the beam spot size. An oscillating pressure was applied
to one of the faces and the from the temperature gradient generated in the process due to strains, one can calculate the work dissipated. The process involves integrating
the gradient of temperature over the entire geometry and taking the time avrage. However, this did not give the correct results, so it was decide to extract the Fourier coefficient
of the signal and perform the integration on it, as done in the case of TR niose by Koji Arai.

I mention the steps it was done -

> The data extracted was stored in a 3D array in matlab. using mphinterp - each dimension for r, z, t

> The number of radial slices and z slices is defined by the user previously

> For each r,z value there was a time signal.

> The fourier coefficient corresponding to the pressure oscillation frequency was extracted from the last three cycles
(In theory only one cycle should suffice)

> The above step was done so that random data( which is small in magnitude) generated by COMSOL is avoided.

> The fourier coefficients are stored in a 2D array corresponding to r and z values.

> The integration was performed using trapz twice once on each dimension to get the total volume integral

> The rest of the calculation is same as the previous script - plugging in the prefactors in the formulae and plotting

**The primary problem with this script is that the data extraction takes significant amount of time - For a mirror with 200 times the radius of the beam spot and 200 radial slices, it takes close to 6 minutes
to evaluate work dissipated for each frequency. The script solves for 16 frequency values. The results for small number of radial slices does not follow the straight line profile. For larger number of slices
the program takes a longer time. The results have not been checked yet.

A common trap in COMSOL:

When you run a time series simulation, COMSOL tries to be "clever" by skipping some time steps
and returns you an interpolated result in order to accelerate the calculation.
But for certain types of applications, such as time series response analysis by applying periodic force/heat,
this "clever" adaptive algorithm messes up our calculation.

You can apply forced time steps by giving a specific option:
- In Model Builder, open your study. Then open your Solver
- Click "Time-Dependent Solver" in Model Builder
- You get bunch of options in Time-Dependent Solver pane.
  Open "Time Stepping".
- Change "Method" from "Generalized alpha" to "BDF".
- Change "Steps taken by solver" to "Strict"

Then, you should be OK.

This can be done by

model.sol('sol1').feature('t1').set('tstepsbdf', 'strict');

in COMSOL with MATLAB. (Of course sol1 and t1 should be changed appropriately.)

 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.4723e-10 J obtained for ratioR=0.71, which is about 3.2% less than the value of Umax=1.5202e-10 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.

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;

 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.

 Great Emory,

Can you put these in a plot so it's easier to read?

 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.

 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

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.

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.

 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.2099e-10 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.

 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.

 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 steady-state offset of the P-controller from the set point develop in the other cases.

In this post I mention of the approach that would be taken to solve for the TE noise spectrum for a beam inside the substrate cavity.

Completing the same would be my plan for the next week.

In the COMSOL simulation, a gaussian pressure will be applied on both the flat surface of the cylinder. The strain would be calculated
using the Structural mechanics module and the temperature field would be calculated using the Heat Transfer module. The approach
is similar to that followed in Liu and Thorne. The Work dissipated would be calculated using the gradient of  T. The spectral density is
then calculated using the FDT.

In the attachment, I have written briefly about the material I have studied this week along with my thoughts on the problem of the TE noise
calculation.

 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.

 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 P-active 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 (steady-state error), and seem to be destabilizing as K_p increases. However, most cases reach steady-state 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.