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.

> 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

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 Abernathy-desktop /dev/tty (v4.3) (ancha/1718 1016), start Mon 2/4 19:52

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 fluctuation-dissipation 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)*(1-e^(-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 non-default 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(1-e^(-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)*(1-e^(-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 fluctuation-dissipation theorem.

Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitational-wave test masses. Physical Review D, 62(12), 122002.

Thermo-refractive noise in a finite cylindrical/infinite test mass with Levin's approach

Location of the codes: 40m SVN repository
comsol/thermo-refractive/

This code realizes Levin's calculation on thermo-refractive 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/S0375-9601(03)00473-0

- The code applies gaussian-shaped 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".

Motivation:

 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:

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₀,  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₀^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 log-log 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.

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 side-by side.

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¹⁰   (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² /Hz) while the Liu-Thorne value is 9.1*10^-40 m² /Hz.

Good Work Deep,

Can you include the equations that you used to calculated these expressions?

 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]

 For calculating the Liu and Thorn U_0 + DU, you need to sum over the first N zeros of the first-order Bessel function. Unfortunately, Matlab doesn't seem to come with a function to do this.

Rather than re-invent the wheel, we can just use the function used by GWINC.

>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).

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.

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.

 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.

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 C-FTM 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.

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.T-T_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 53-61]. 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 end-point 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 121-122

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');

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 quatity '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.

You'll notice that the range command that is given to COMSOL is in single quotes, that means that all Matlab is doing is feeding to COMSOL a string, just as you would send it numerical values as a string.

It is COMSOL that evaluates this range command, so it can be different from Matlab's range command.

 **The point where the feature called 'Transient' is spoken about - By selecting 'Time Dependent' under the 'Study' in COMSOL desktop, the line that gets added to the Matlab script.

So it was the 'Time Dependent' feature that was selected while creating the model

The file 'thermo_refractive_COMSOL_1D_axisymmetric.m' found in the SVN repository https://nodus.ligo.caltech.edu:30889/svn/trunk/comsol/thermo-refractive/ performs the data extraction from the COMSOL simulation

and computes quantities as given in Heinert's paper.

This is the where the dissipation is calculated from the temperature gradient and the use of FDT is made to evaluate the linearized PSD.

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The first few lines of the code extracts the variable values out of the structure 'param' and stores it in new variables for the ease of coding.

The important portion comes in with the iterative structure (for loop) in line 29.

A point wise summary of the steps done would be as follows -

* The finite elements in space and time are defined in the arrays r0 and t0 according to the time and radial slice steps dt and dR until R, the radius and the t_end, the time until which simulation runs

* A matrix called dat is created which is used to store the values of the temperature gradient as returned by COMSOL.

It stores the values as a function of r along the columns and as a function of t along the rows.

Thus, moving vertically down the matrix in a column would give the values of the temperature gradient for a fixed r as a function of t

* Two other row vectors datrabs and datrphs having lengths of the array r0 are used to store the values of the square of the Fourier coefficients and the phase angle.

* The extraction of data from COMSOL is done in the lines 71-74

    for i2=1:length(t0)
       tmp = mphinterp(model,'ht.gradTr','coord',r0,'T',t0(i2));
       dat(:,i2) = tmp;
    end

* After the above step, the matrix dat is filled with the values of the temperature gradient as mentioned previously.

* Next to evaluate to the Fourier coefficients, the second half of the simulation time is used. Probably to let the transients die down as mentioned in the ELOG post by the author.

* To find out the Fourier components, we multiply the desired function by sin(omega*t) or cos(omega*t) and integrate over the period of the function. There is also a prefactor of 1/L where L is the length of the period

over which the function is defined.

* To do the above exercise, the sine and cosine of omega*t is separately evaluated and stored in skernel and ckernel. Note that the time interval as mentioned above is the second half of the simulation time.

* Now for all the radial slices, the Fourier coefficients are extracted using the statements in for loop in line 83

    for i2=1:length(r0)
       tseries=detrend(dat(i2,(length(t0)+1)/2:length(t0)));
       tmp1=mean(tseries.*skernel);
       tmp2=mean(tseries.*ckernel);

       datrabs(i2)=2*sqrt(tmp1^2+tmp2^2);
       datrphs(i2)=atan2(tmp2,tmp1)/pi*180;
    end

Using the definition of mean of a function ' f ' in the continuous case as f_mean = integral{ f(t)dt} / integral{ dt } , over the values of t

tmp1 has used ' f ' as tseries.*skernel which means that one multiplies the time signal with sin(omega*t) and integrates. tmp2 has used cosine instead

Thus, tmp1 and tmp2 are the Fourier coefficients of tseries at the frequency fmod defined right close to the first for loop in line 29.

* The coefficients are squared and added in the line 88

       datrabs(i2)=2*sqrt(tmp1^2+tmp2^2);

This quantity how much of the frequency ' fmod' is present in the time signal.

datrphs stores the same for the phase shift phi. Note however, all of this is done as a function of r. The plotting is done in the following lines. Note that while running the code, as a result, the plots are seen to
get updated each time the loop [line 29] iterates.

* Following this, is the part where W_diss, the dissipation is calculated using formula [The expression used does not match Eq. 15 in Heinert's. Any comments regarding the formula used for W_diss are welcome].

The heinert's eq(15) goes as - W_diss = pi*H*kappa/T0 * integral{ grad_T^2 * rdr }

The quantity datrabs is squared and the integral mentioned above is calculated w.r.t. r using the trapz algorithm from Matlab.

The value for each frequency i.e. each 'fmod' is converted to a string and displayed on screen.

* The reason why the grad_T was changed to its frequency domain before integration is suspected to be something related to the Parseval's Theorem. However, more details or correction on the same is welcome.

In this post, I am reporting of the relative difference between the analytic and COMSOL results.

The file by the name 'Main_thermo_refractive.m' found in the SVN was slightly modified to generate the difference plot.

The quantity plotted is the absolute value of the difference between the two results divided by the result given by COMSOL.

The plot is in a semilogx plot. It can be seen the maximum deviation goes to a maximum of 6%.

Also the difference is higher at the higher frequencies.

The plot shows two cases - the case for which 500 terms in the series were summed to get the analytic solution and second
is the case where 1000  terms were summed. The plot does not change much due to the number of terms being summed
showing that the convergence is fast.

This week I will update my previous eLog entry with a nicely labeled plot showing both lines on the same plot, then I intend to implement improved meshing on the mirror in the COMSOL model I have constructed, increasing the number of elements on the mirror's face and the central axis of the mirror such that there are more elements in the regions where the applied force is greatest.  This should result in more accurate and faster to compute simulations, allowing us to increase the number of elements in the mesh and possibly reduce the residual and allow us to use the better gravitationally balanced boundary condition and obtain an answer which converges.  Afterwords, I intend to look at sample code showing how to use Matlab in conjunction with COMSOL since we will need to do so for my project.  If there is time remaining later in the week, I will attempt to replicate some of the results Steve Penn generated in order to verify my model.

In order to make computations more efficient and possibly allow the set of boundary conditions based on Liu and Thorne's suggestion of a gravitational force preventing bulk motion of the mirror (described in http://nodus.ligo.caltech.edu:8080/COMSOL/3) we want to be able to more optimally mesh our structures.  In particular, we would like to have a finer mesh on the face of the mirror and especially near the center of the mirror.  Doing so will allow us to significantly reduce the total number of elements contained in the mesh while keeping a large number in the regions which require them.

The best way to improve the meshing that I have found so far is to introduce a new geometry and mesh around it.  In order to test this method, I constructed an additional cylinder in my model centered at the same location as the mirror.  I gave the cylinder the same height as the mirror so that they can be changed by modifying a single variable and set the cylinder's radius equal to the Gaussian beam size of the laser.  I then constructed a user-defined mesh composed of two free tetrahedral operations.  The first of these I restricted to the smaller cylinder.  By selecting a distribution using a fixed number of elements, this region can be meshed fairly uniformly and densely.  The second meshing region can be specified again using the inner cylinder's top face as well as the edges running along its side.  For this region, it seems optimal to use a predefined distribution type and select a geometric sequence.  Using this method provides a much finer mesh in the areas where we want one, without wasting computational power performing more calculations in uninteresting regions.  A screenshot demonstrating the results of this method is shown below.

I ran the previous analysis using the simplest boundary condition of a fixed back using this new mesh.  The newly constructed mesh contained 7038 elements instead of the previous 7243 and obtained a Umax=1.57985*10^-10 J instead of 1.52887*10^-10 J.  This is about a 3% difference in results.  In the next few days, I will run the simulation again on a computer with more RAM using a finer version of the original mesh in order to confirm that the results of the new mesh agree better with it than with the current mesh with a similar number of elements.  As a more immediate test, I constructed a mesh of 1940 elements using the new method and obtained a result of Umax=1.56085*10^-10 J which is closer to the value obtained using the original meshing technique, though far enough away to indicate that they may not agree, which encourages me to run the original mesh design with more elements.

The next thing to consider is further improving the mesh by increasing the element density near the mirror's front face.  From what I have seen it may be more difficult to implement both of these improvements together, in which case we can perform testing to determine which of the two methods provides better computational efficiency.

Liu and Thorne demonstrated in their paper that the optimal set of boundary conditions for a problem like the one I am working on is the balance the force applied to one of the mirror faces by an opposing force of equal magnitude distributed throughout the body of the mirror.  When we first attempted to implement this boundary condition, we obtained an error from COMSOL indicating that a solution could not be generated as "The relative error is greater than the relative tolerance."  I attempted to find a solution to this problem so that we can run our simulations using the correct boundary conditions.

I tried a number of settings changes and tweaks suggested on the COMSOL forums and in official documentation on the error (https://community.cmc.ca/docs/DOC-1453), but without success.  Most of the suggestions were not relevant to our model, and confusingly when the number of elements in the mesh was increased from 13628 to 24570 the relative residual increased contrary to our initial attempts in http://nodus.ligo.caltech.edu:8080/COMSOL/3.  It seems like the next action to take is to attempt to generate a simpler model and attempt reproduce and fix the same error.

I also decided to run both an eigenfrequency and frequency domain solver on our model using improved meshing with 24570 elements.  I first ran these solvers using the rigid back constraint and both of them returned Umax=1.49931*10^-10 which gives S_x(100 Hz)=1.97673*10^-40 m^2/Hz which agrees quite well with the solution obtained from our previous simulations which gave S_x(100 Hz)=2.0157*10^-40 m^2 / Hz.  This still significantly differs from the analytically obtained value of 7.80081*10^-40 m^2 / Hz using Liu and Thorne's calculation.

I then changed the boundary condition to use Liu and Thorne's gravitational body load boundary condition and both of these methods returned Umax=1.17164*10^-9 J which corresponds to S_x(100 Hz)=1.54472*10^-39 m^2/HZ.  This only differs from the analytic result by a factor of 1.98, which is a significant improvement over the previous simulations, but still a much more significant difference than I would hope for given that they now use the same boundary conditions and it seems unlikely that using a stationary state solver will give a different result when both the eigenfrequency and frequency domain solvers agree to 6 digits.

edit: Taking the advice of the error documentation, I also ran a time dependent solver which returned the same result as the eigenfrequency and frequency domain solvers.

 We have been encountering an error in COMSOL for a while of the form "Failed to find a solution.; The relative residual (###) is greater than the relative tolerance.; Returned solution is not converged.; -Feature: Stationary Solver 1 (sol1/s1); -Error: Failed to find a solution."  The error has occurred when attempting to find a stationary solution in models with boundary loads and no fixed constraints (preventing an edge from moving).  We wanted to determine what the error was caused by to allow us to use the stationary solver, or at least to confirm that the error was not indicative of a problem in our model.  To this end, I designed a few very simple COMSOL models in which I was able to reproduce the behaviour and attempted to determine the root of the issue.

I first constructed a somewhat similar model to ours using a cylinder of fused silica with all of the default values and a normal meshing.  I applied a boundary load of 1N on one of the faces and ran a stationary solver.  As expected, the solver failed to converge since it had no boundary condition preventing it from accelerating continuously.  Applying a force of 1N on the opposing face resulted in the same error as above, which replicates the previous error since the system is failing to converge in a case where it should.  I decided also to make an even simpler 2-D model of a square.  Applying 1N forces to opposing sides on the square again returned the error above.  

Both of these models were able to be evaluated using at least an eigenfrequency solver as noted on the primary model in the previous eLog.  I looked on the COMSOL forums and read through some more of their documentation and saw the suggestion in response to this error to use a time-dependent solver and simply view times after the system will have settled to a stationary state (#2 https://community.cmc.ca/docs/DOC-1453).  I attempted this on the test models and both of the time dependent solutions converged without error to their expected solutions (compression between the faces on which forces were being applied).  This may be a sub-optimal computational method though as even in the simple cylinder case with 6133 elements and a simple force profile, the solution took several minutes to run.  For the cylindrical model, I evaluated the strain energy  using both an eigenfrequency and time dependent solver and obtained the same result using both of the solvers.  The eigenfrequency solver evaluates much more quickly than the time dependent solver, and in the primary model as I noted in my previous eLog, the strain energy  obtained using the eigenfrequency and frequency domain solvers agreed, so it seems that the best manner in which to proceed is to use the eigenfrequency solver to compute the strain energy.

The source of the error is still unknown, but given these tests, it seems very unlikely to be indicative of a problem in our model.  We still have a very significant disagreement between the simulated results and the calculated values.  I am going to spend the next day or so looking through both the COMSOL model and the analytic calculation and checking them for errors which could cause this discrepancy.  I will then start reading the documentation on Livelink for Matlab and try to implement it.

Emory:

I don't understand how you're getting the strain energy from the eigenvalue solver. It is my understanding that the eigenvalue solver will only give you the strain energy at a particular eigenfrequency. We're interested in the strain energy from the beam deformation at frequencies below the first eigenmode. 

While discussing TE noise with Yanbei Chen, it was realized that the expression for the work dissipated W_diss was derivable from the inhomogenous Heat equation with a source term.
The exercise was tried out to some success and can be found in the attachment. The attachment describes the steps roughly.
The important point to note is the fact that while Liu and Thorne considered stresses developed in the material by means of the heat balance equation, they have ultimately resorted to the
expression for W_diss = 1/T * integral{ kappa * grad(T)^2 rdr } to calculate the dissipation. They have used an expression whch relates the expansion, Theta, to the gradient of temperature.

It is suspected that the fundamental approach is to consider a source of heat and evaluating the dissipation. However, Liu and Thorne considered applying pressure probably because of the
physical scenario of the fluctuations in the mirror face.

Discussion over the same would be helpful.

I used the same integrator that we had setup for use with the stationary solver.  I had expected it to either fail or return a very different result, but assumed that the fact that it was returning the same result as when it was used after the previous solvers indicated that it was able to use that solver to find the value.  After seeing this comment I went back and checked and the solver was setup to use the last stationary solution regardless which solver last ran.  I will run the previous tests which relied on using this integrator in conjunction with non-stationary solvers again and see if they actually agree and how the results they give compare to the analytic solution.  I will spend some time today trying this.  Unfortunately, both the time dependent and frequency domain solvers take a long time to converge, so running these tests may take most of the day, but I can start reading through some of the Livelink for Matlab documentation while they run.

edit:

I ran the time-dependent solver for the rigid back boundary condition and computed the strain energy for the final two time-steps.  Both of them gave a value of 1.58425*10^-10 J which corresponds to Sx(100 Hz)=2.08871*10^-40 which is about 3.6% greater than the value obtained from the stationary solver.  I don't understand why these values would differ at all.  When I tried to run the time dependent solver with the spring back and gravitationally opposing force boundary conditions an error message was returned: "Failed to find consistent initial values., Last step is not converged.; -Feature: Time-Dependent Solver 1 (sol4/t1).; -Error: Failed to find consistent initial values."

I'm going to spend some more time looking at the COMSOL model attempting to find anything which could be causing this error and reading any relevant documentation I can find.

edit 2:

I spent most of the day attempting to find the source of either of these errors.  Possible solutions I found and tried included increasing or decreasing how fine the mesh is, increasing the acceptable tolerances, and increasing the time interval in the time dependent solver.  I attempted all of these and none of them worked.  Surprisingly, when I increased the acceptable tolerance level for the stationary solver, it returned a greater relative residual which does not make sense to me.

I also took the simple 3d cylinder constructed yesterday and was able to replicate the errors with it.  When I increased the time interval on that case for the time dependent solver, it converged and gave a result.  I was able to get a Umax value for the final time step, which should be equivalent to the one which we would expect a stationary solver to return.  Further increasing the number of time steps the primary model computed did not cause it to converge.  I will run this again with many more time steps and see if it converges.  Even if it does, this doesn't seem like a good way to do our computations as it takes a long time to complete a solution, but seeing whether it converges may give us information on what the problem with the stationary solution is, in particular if increasing the number of time steps does cause it to converge, then that would indicate that finding a way to make the simple case converge would probably work for our model as well.

In the code for evaluating TR noise, the parameter for automatic meshing (which can take values from 1 to 9 -  1 being the finest and 9 the coarsest mesh) was changed from 1 to 5 and 9.
The computation time surprisingly doesn't differ much. The maximum relative error also stays close to 6% as in the case of the finest mesh. The relative difference for two cases for a mesh
number of 5 and 9 are shown in the attachment.

The conclusion being that it is better to stick to the mesh number 1 which is the finest mesh.

  We have been seeing an error when attempting to use a stationary solver in conjunction with a set of boundary conditions which does not fix the face of the mirror opposite the applied force.  I have tried a number of settings changes and tweaks in order to attempt to get the stationary solver to converge on our model.  I have found that by switching to the PARDISO solver and increasing the relative tolerance to 500, the solution will converge (Umax=1.49866*10^-10 J).  This does require increasing the relative tolerance greatly, which is concerning.  After doing this, I also found that by increasing the relative tolerance to 700 and using the SPOOLES solver it also converged giving Umax=1.498653*10^-10 J.  The fact that these agree quite well indicates that the increase in relative tolerance may not have harmed our values.  If this were a more complicated system which we would expect to have behaviour which could cause our solvers to get stuck on a set of values I would be more concerned, but I think that this may be a workable fix in this case.  These values of Umax give Sx(100 Hz)=1.97586*10^-40 m^2/Hz which differs by about 5% from our previous value of Sx(100 Hz)=2.08291*10^-40 m^2/Hz.  Unfortunately, this seems to indicate that the difference in results between the COMSOL model and our direct computation is not due to a difference in boundary conditions.  I will spend some time looking for more useful documentation on the relative residual, relative tolerance, and LU factorization: out of memory (despite having more RAM availiable) error, then I will work through our COMSOL model again, possibly remaking it, and check it for any errors which could result in the disparity between our simulated and directly computed results.

What about the dependence on the size of the time slice?

It is unclear what the unit fo the frequency. Hz? ot logHz (which should be expressed as 10^1, 10^2, ...)?

Can you combine the result of the caluculation and the error together?
i.e. Combine the plot similar to http://nodus.ligo.caltech.edu:8080/COMSOL/4 and your plot with the horizontal axies lined up.

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

In this post, I mention of my immediate plan for the week and also describe the problem I am looking at in brief.

In the calculation of TR noise for the beam passing through a cylindrical substrate, Heinert has considered thermal fluctuations
by means of an sinusoidally oscillating heat source, scaled according to the intensity profile of the beam, present along the length of the cylinder.
The heat equation is then solved to evaluate the temperature field. The gradient of temperature is used to calculate the work dissipated.
The work dissipated then is used to find the spectral density according to the Fluctuation Dissipation Theorem.

On the other had, the case of TE noise is considered in the work by Liu and Thorne where the beam reflecting from one of the faces is
considered. In their work, they have used an oscillating pressure, scaled according to the beam, at the face. The equation of stress balance
is solved to get the strain field, the strain results in heating and as a result, a temperature gradient. The work dissipated is calculated using
the gradient of temperature and the spectral density is calculated using the FDT.

Now, in case of Heinert's work, it is to be noted that the heating term is the one that contains the TR coefficient, beta = dn/dT. This is the where
the TR noise that we are looking into comes in the picture. Liu and Thorne's work on the other had, never has an explicit heat source. In their case,
it is the strain that generates the heat implicitly. They have related the expansion with the temperature perturbation in their paper from an expression
in Landau Lifshitz's, Theory of Elasticity. The important point to note is that the physical parameters that characterize the TE noise like coefficient of
linear expansion, alpha, or Poisson ratio, sigma, come into the picture through this relation. Another point to note is that the expression
for the work dissipated ( W_diss ) uses the gradient of temperature( It is the same formula that Heinert has used). This expression is derivable from the heat
equation. Thus, one could have also done the exercise by injecting the right heat source and solving the heat equation instead, since its ultimately
fluctuations in temperature that cause these noises TR or TE.

Our problem is to evaluate TE noise for a beam that passes through the cylindrical substrate instead of reflection off the face. It is suspected that using an
oscillating pressure on the surface will not be the correct approach since the beam is going through the material and not just reflecting from the surface. We
want to solve it by means of a heat injection, as done in Heinert, calculating the gradient of temperature, the work dissipated and then the spectral density. It is realized that the
heat source should be oscillating but the correct coefficients is what is undetermined i.e. we realize the heat source is q_dot = [- - -] * cos(omega * t) . In case
of Heinert it is at this point that the 'beta' comes in. However, in the TE case this is not yet determined. The literature doesn't deal with the case of beam goin through
a material. The equations in Heinert must be looked at more deeply to realize how the 'beta' comes in and then drawing an analogy, we may be able to figure out the
right heat source for the TE noise case.

Any comments on references,  the approach that should be taken, or any thoughts on the problem is most welcome.

After setting up a COMSOL model that includes the heat flux from the laser and the ring heater, I've made the model solve over time manually by performing the solution process over a loop in MATLAB.

This allows for the future insertion of the PID controller object in the solving process, and the dynamic manipulation of the applied heat loads.

The following is an automatically generated plot of defocus effects as a function of time at 1 beam radius (54 mm), included in the program, with only the ring heater turned on in a top-hat emission configuration and the total power being 5 watts.

The linear projection values needed to calculate the defocus effect are extracted directly from COMSOL, with no output files required through the use of the mphinterp command.

The behavior appears to be physical and is of the correct order of magnitude.

I've also attached the code that produced it for verification. (It requires COMSOL+Matlab to run.)

 I contacted COMSOL support about our difficulties with the relative residual error and was told "Typically for solid mechanics this error is because you have not constrained the body enough. So any solution + a rigid body transformation is also a solution. To remove this non-uniqueness, you need a solutions modulo the rigid body transformations. So try and constrain the body somehow."  This is what the gravitationally balanced body load was supposed to do, but using the stationary solver has been non-convergent and the eigenfrequency solver has generated odd modes, with the first real modes being sheer modes.  These outputs indicate a problem with our model.  Matt noticed yesterday that the meshing in our model was slightly non-symmetrical, but we initially dismissed it as not being significant since it was a minor difference which we did not think could account for the large errors we are facing.  At further consideration though, if some asymmetry arising from any source, even a numerical or rounding error, were present, that could cause a slight rotation of the object which would cause the object to experience a torque and minor sheer force.  The stationary solution would not converge in this case, but the eigenfrequecy solution might, and if it did it would make sense to see sheer modes for some of its low frequency eigenmodes, as we have.  One possible solution to this problem is to change the way we mesh the material to ensure a symmetrical distribution of nodes in the x-y plane, probably by extruding lower dimensional meshed systems into our model.  I am unsure if we would be able to implement this solution once we start to change the size of the radii of the object's faces.  An alternate solution is to find another set of boundary conditions which should be equivalent to the gravitational body load constraint, but which are stable relative to minor perturbations of the system's conditions.  I think that I have found another set of boundary conditions which should work and not be too difficult to implement in COMSOL.

To verify the validity of the solutions produced by the manual simulation, I've calculated the relative error between the results from the manual code and the results produced by COMSOL normally.

The plot for the relative error in the defocus at r = 0 and r = 54 mm is shown below, in the case where only the ring heater is turned on at a total power of 5 W.

The following indicates that the maximum error is less than 0.01% (in percent error format).

I've attempted to determine the solving time per loop as a function of the simulation time, in an attempt to identify any trends in the solving time for a constantly dynamic load.

The following is a plot of the solving time per loop as a function of the simulation time for a load which is constantly dynamic (sinusoidal in time, in this particular run).

The mesh size is normal (default in COMSOL) with heat fluxes from both the beam and the ring heater (as in the real case).

It is difficult to identify any particular behavior or trend due to the large amount of "noise" other than a trivial general increase after ~4 s. Mesh quality does not appear to influence solving time per loop significantly.

Work must be done to reduce the total solving time for the simulation, which in this case amounted to 18 and a half minutes (1.1109e+03 seconds).

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

 In our calculation of Sx(f) from Umax prior to this, we were missing a factor of 4.  The correct formula from Liu and Thorne's formula 58 is Sx[f] := 4*Kb*T*Umax*Phi/Pi*f*F0^2).  When we correct this formula, we get results which agree fairly well with the analytic results.  The Mathematica script attached was used to perform both calculations of the thermal noise spectrum (using SxComsol and SxFTMThorne).

SxCOMSOL(100 Hz) = 7.99735*10^-40 m^2/Hz

SxFTMThorne(100 Hz) = 7.80081*10^-40 m^2/Hz

These results differ by about 2.5%.  We also need to verify that the result converges for increasing mesh size.

I attempted to run the simulation with 24360 elements in the mesh, but the computer I was running it on repeatedly crashed.

edit: Running on a more powerful computer I got the following additional values.

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

 I have inserted a rudimentary PID function into the manual simulation code as a way to test whether or not the PID function is changing the defocus values in the desired manner.

I am currently determining the ratio of ring heater power to the steady-state defocus as a way of measuring the scale of the response.

This ought to give a good way of measuring the scale needed to convert the calculated actuator response into an actual load.

I've attached the rudimentary code below. (The actuator isn't feeding into the heater at the moment, but inserting the "actuation" variable into the load expression is all that is required.)

 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.

 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.

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.

 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.

 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.

 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.

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.

Liu and Thorne, in their paper have analytically solved for the TE noise spectra considering the test mass to be a half infinite space.

An attempt was made to get similar results using a COMSOL model. In this post, the rough profile of the PSD plot is given and is compared
to Liu and Thorne's result.

Steps followed by Liu and Thorne:

> Oscillating pressure, scaled according to beam profile, is applied one of the surfaces.

> The stress balance equation is solved borrowing algorithm from Landau Lifshitz.

> The stress and the temperature perturbation is related using an equation given in Chapter 1 of LL

> The W_diss is calculated from the gradient of the temperature.

> PSD is calculated using FDT due to Levin

> The resulting plot looks like the following - ***Only the frequency dependence is plotted here leaving  out the prefactors which will
appear just as an additive constant in the figure. The aim was to see how the plot of the simulation looks like.

Steps followed in the COMSOL model:

> 2D axisymmetric model was used since the applied was applied radially symmetrically.

> The dimensions of the test mass(cylinder) was defined larger as compared to the beam spot radius.

> The physics added was Solid Stress and Heat Transfer in Solids.

> Thermal expansion was enabled in the Solid Stress physics

> Pressure work was enabled in Heat Transfer physics which takes into consideration heat due to stresses i.e. adds this heat source to the heat eqn

> Heat insulation wsa applied to all surfaces.

> All surfaces apart from the one where pressure was applied was held fixed to mimic infinity

> A custom mesh was put.

> The model was exported as '.m'

> It was called for frequency range 10^(0 : 0.25 : 4)

> Following Liu and Thorne, the temperature gradient was integrated over the volume using mphint2() function in Matlab

> The time average was taken over the second half of the cycle to avoid transient which are assumed to die off after a couple of initial cycles.

> All the prefactors were ignored while plotting

> The plot shown below is integral{grad_T}^2} / omega^2   which is proportional to PSD

> We just wanted to see if at all we got anything looking similar to the analytic result of Liu and Thorne (shown above)

The above plot is when the radius is 100 times the beam spot and height of the cylinder is 200 times the beam spot
Other values of radius and height were also tried. The one below shows similar plot as above but for radius 10 times
the beam spot and height 20 times the same

While the one below is for radius = 1000 * beam spot and height = 2000 * beam spot

Here, it is seen that the straight line with a negative slope is seen to some extent.

However, the prefactors, when put in, will show the actual PSD.

Any comments or suggestions related to the model construction in COMSOL or otherwise is welcome.

 Here is a plot showing the profiles of Liu and Thorne and the COMSOL simulation. Since the constants were not taken into account, one expects the plot to be shifted by a constant equal to the logarithm of the prefactor that has been ignored.

The green 'o' is just a 1/omega^2 plot in the loglog scale which is how Liu and Thorne's result go (this is not linearized PSD, but is doesn't matter as long as we want to see the profiles). The other plot in blue is the COMSOL simulation.

The plot below is a difference between them. One expects it to be a constant but there are a few jitters. Once again, it is emphasized that the numerical values do not give any information related to the actual difference which will only be calculable
once the correct physical constants are put in the formula. I will post another entry once that is done

 I attempted to run the Matlab code built last week, but received an error that the mesh failed to generate on the inner domain/failed to respect edge element on geometry edge.  The error occurred on the second computation with ratioR=0.2, but for ratioR=0.1 the simulation completed successfully.  I reproduced the error in COMSOL and noted that for meshes with size at least as course as "course" or extremely fine, a mesh could be constructed, but for sizes in between an error was returned.  This region of sizes which throw errors varies for different values of ratioR.  By specifying an edge selection of the central line running down the frustum and another along the top edge of the inner cylinder and meshing points on those edges first the rest of the meshing is forced to use those nodes as a seed.

I also noted our COMSOL model had an error: ratioR was applied to the wrong face.  This was corrected as was the .m file.

I added an eigenfrequency study and ran it for several values of ratioR, obtaining the following:

ratioR Lowest Drumhead mode (Hz) Lowest Butterfly mode (Hz)

0.1 3624.8 3802.4

0.3 4375.9 3785.4

0.5 4366.5 3650.0

.75 3709.0 3041.85

1 3000 2186.3

1.1 2736.2 1893.0

1.3 2247.6 1426.4

The mirror's front face is 0.17 [m].  Its Gaussian beam size is 0.0156 [m].  Its height is given by 0.01734[m^3]/(R^2*(1+ratioR+ratioR^2)) such that when the two faces have the same radii it is 0.2 [m] long.

 Great Emory,

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

 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.