40m QIL Cryo_Lab CTN SUS_Lab TCS_Lab OMC_Lab CRIME_Lab FEA ENG_Labs OptContFac Mariner WBEEShop
  COMSOL elog, Page 1 of 3  Not logged in ELOG logo
ID Date Author Typeup Category Subject
  1   Thu Dec 13 18:47:18 2012 ranaGeneralGeneralCOMSOL: who's using up all the licenses?

 When you want to find out who's using up all the licenses, you can run lmstat -a do find out.

Specifically, with Mac COMSOL, I do:

> cd /Applications/COMSOL43a/license
> maci64/lmstat -a -c license.dat

Users of COMSOL:  (Total of 5 licenses issued;  Total of 1 license in use)

  "COMSOL" v4.3, vendor: LMCOMSOL
  floating license

    dmassey c22042.local (v4.3) (ancha/1718 1074), start Thu 12/13 18:26

etc.

  2   Tue Feb 5 13:15:13 2013 DmassGeneralGeneralCOMSOL: who's using up all the licenses?

Quote:

 When you want to find out who's using up all the licenses, you can run lmstat -a do find out.

Specifically, with Mac COMSOL, I do:

> cd /Applications/COMSOL43a/license
> maci64/lmstat -a -c license.dat

 

Users of COMSOL:  (Total of 5 licenses issued;  Total of 1 license in use)

  "COMSOL" v4.3, vendor: LMCOMSOL
  floating license

    dmassey c22042.local (v4.3) (ancha/1718 1074), start Thu 12/13 18:26

etc.

 

Whenever I try to open a model (.mph file) which I built in COMSOL and has never been touched by any CAD software, it uses the CADIMPORT license, and we only have 2 of these. I can open a different mph file without the CADIMPORT being used.

I don't yet understand what makes something you build in COMSOL and save in COMSOL try to use the CADIMPORT module.

CADIMPORT licenses currently in use: (2/2)

    ligo m5 m5 (v4.3) (ancha/1718 202), start Wed 1/30 11:39
    mattabe Abernathy-desktop /dev/tty (v4.3) (ancha/1718 1016), start Mon 2/4 19:52

  4   Thu May 2 14:00:36 2013 KojiGeneralConfigurationtest mass TR with Levin's approach

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

Attachment 1: thermo_refractive_1D_axisym_result.pdf
thermo_refractive_1D_axisym_result.pdf
  12   Wed Jun 5 20:54:59 2013 Deep ChatterjeeGeneralGeneralBessel Function roots

During the process of evaluating the PSD from Sec. V of Liu and Thorn, I chanced to write a simpler root finding algorithm applying bisection to find the roots of J1(x).

The difference between this and the algorithm by

Greg von Winckel goes as
Difference_between_algorithms.png

The difference is of the order 10-7 and can be reduced by reducing the tolerance. However, It should however be noted that bisection is a crude algorithm for rough usage and differences become pronounced for larger n.

 

Attachment 2: bessel_zeros.m
function X = bessel_zeros(k, n)
%BESSEL_ZEROS : calculates the first k zeros of the bessel function Jn
%   k : No. of ROOTS to be evaluated
%   n : Order of the Bessel function Jn
%   X : stores the roots in serial order i.e. X(1) gives the first root,
%   X(2) gives the second and so on
X = zeros(1,k);% empty array of lenght k
count = 1;%this acts as a counter to k
dx = 0.1;%step size within which, by assumption, no roots exist
x = 0;
... 22 more lines ...
  15   Fri Jun 7 17:41:41 2013 Deep ChatterjeeGeneralGeneralDisscussion of the code of a 'comsol with matlab' model file

In this post the matlab code to build a model using COMSOL with matlab is analyzed from Koji Arai's codes on the calculation of TR noise in the 1D axisymmetric case. The post primarily describes the various commands used to interface with COMSOL.

It should be noted that the matlab script analyzed here is not the master script that will perform the simulation. The "main" file is called the "Main_thermo_refracive.m" found in the SVN repository TR noise_comsol

>The main script defines a structure called 'param' which stores all the parameters including the material properties and those related to running the simulation in comsol. This structure 'param' is hence passed on as arguments to all the other functions (in the SVN rep)

that compute the analytical solution or the COMSOL simulations.

>The part of the code ( the second if...end block) in "Main_thermo_refracive.m" that calls the COMSOL simulations, calls the function "thermo_refractive_COMSOL_1D_axisymmetric".

>In this function, the parameters are extracted from the structure 'param' and stored in separate variables. The time step in the simulation, the end time, mirror radius are defined and the function

that generates the model "thermo_refractive_COMSOL_1D_axisymmetric_model" is called which prepares the model in COMSOL.

> it should be noted that it is much of much ease to the user to define a simple model in COMSOL and export it in Matlab and then analyzed the same.

I had created a simple such file which a simple 3D metal bar is created and added some physics and study to see how the matlab file.

I have attached the same which also describes the relevant packages and subroutines in the comments. However, I had not erased the COMSOL history before the export hence it has significant amount of superfluous code. One may however have a look at the  same

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

>The first thing done by the model builiding function is that it extracts a few parameters and stores them in separate variables in the first lines

f_mod = param.COMSOL.simulation_modulation_freq;
dt    = param.COMSOL.simulation_time_step;
t_end = param.COMSOL.simulation_end_time;

 

>Since COMSOL accepts string arguments, the time steps are defined as a string

trange = ['range(0,' num2str(dt) ',' num2str(t_end) ')'];

**There is however a query regarding the keyword range which in Matlab returns the difference between the maximum and minimum of the list passed to it. range() passed to COMSOL would probably create an array fro 0 - t_end in steps of dt.

The above line converts the quantity 'dt' and 't_end' to a string and concatenates it along with the rest of the string to make a valid string that would be understood by COMSOL.

>The COMSOL class is called next by the following two import statements (after a series of displays on the screen showing the status of the simulation).

>Managing models( like creation and destruction) is handled by the modelUtil class in COMSOL with Matlab and hence to create a model modelUtil.create() is called.  "model.remove()" is used to destroy one such object. Line 29 says -

model = ModelUtil.create('Model');

>The above statement created an object called 'model' in Matlab and and a model in COMSOL by the name of 'Model'. The various attributes of this object is defined by the next few statements. Like,

[line 31] - model.name('thermo_refractive_COMSOL_1D_axisymmetric_model.mph');

The above statement assigns the filename of this COMSOL model created. Note that, since in Matlab, the Matlab object 'model' is used to invoke the functions i.e. 'model.name()' and not 'Model.name()'.

> The 'model.param' contains all functions related to setting and describing the parameters in the model.

*Note that the 'param' following keyword 'model' connected by the '.' has nothing to do with the 'param' structure used in the script "Main_thermo_refracive.m" which have been named anything else

> Now the model.param.set(<P>,<expr>) is used to give the parameter P an expression expr both being string.

The model.param.descr(<P>,<des>) has the same format but gives the description of the parameter P as des (something that is understandable in common terms).

This can be understood in the statements in [line 46 - 51]

model.param.set('beam_shape', '2/(pi*beam_size^2)*exp(-2*r^2/beam_size^2)');
model.param.descr('beam_shape', '');
model.param.set('beam_intensity', 'beam_power*beam_shape');
model.param.descr('beam_intensity', '');
model.param.set('dT', 'mod1.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');

 

Attachment 1: Iron_Bar_test1_model.m
function out = model
%
% Iron_Bar_test1_model.m
%
% Model exported on Jun 10 2013, 10:13 by COMSOL 4.3.0.151.

% This is an exported file from COMSOL in the '.m' format. I am adding
% comments to the lines.
% As one goes along the steps one followed in COMSOL and the corrsponding
% Matlab file, the purpose of the various subroutines become clear. The
... 430 more lines ...
  16   Mon Jun 10 18:37:09 2013 Matt A.GeneralGeneralResponse to question in: Disscussion of the code of a 'comsol with matlab' model file (partially complete).

 

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.

 

  17   Tue Jun 11 10:21:05 2013 Deep ChatterjeeGeneralGeneralRe: Disscussion of the code of a 'comsol with matlab' model file

Quote:

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

model.study('std1').feature.create('time', 'Transient');   [line 134]

> The time step tweaking of COMSOL is stopped by the following line 166

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

 

 **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

  18   Tue Jun 11 17:04:33 2013 Deep ChatterjeeGeneralGeneralAnalyzing the file 'thermo_refractive_COMSOL_1D_axisymmetric' by Koji Arai

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.

  20   Mon Jun 17 11:06:17 2013 Emory BrownGeneralGeneralPlans for the week of 6-17-13 through 6-22-13

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.

  23   Wed Jun 19 14:59:13 2013 Emory BrownGeneralGeneralThe Relative Residual Convergence Error

 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.

  24   Thu Jun 20 10:26:38 2013 Matt A.GeneralGeneralThe Relative Residual Convergence Error

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. 

 

Quote:

 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.

 

  25   Thu Jun 20 12:41:01 2013 Deep ChatterjeeGeneralGeneralThe expression for the "work dissipated" in the TE and TR noise calculations

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.

Attachment 1: dissipated_work.pdf
dissipated_work.pdf dissipated_work.pdf
  26   Thu Jun 20 12:51:44 2013 Emory BrownGeneralGeneralThe Relative Residual Convergence Error

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.

 

Quote:

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. 

 

Quote:

 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.

 

 

  28   Sun Jun 23 14:00:05 2013 Emory BrownGeneralGeneralManipulating the Relative Tolerance

  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.

 
After doing this, I was able to find more information on the relative tolerance in the COMSOL release.book (page ~30 https://www.tuegrid-doc.uni-tuebingen.de/dokuwiki/lib/exe/fetch.php?id=hpc-uni%3Asoftware-docs%3Acomsol%3Acomsol&cache=cache&media=hpc-uni:software-docs:comsol:release.pdf).  It appears that the relative tolerance and the Factor in error estimate values jointly determine the maximum allowed difference between successive estimates when using an iterative solver.  I decided to try to get a better idea for this behaviour, so I increased the Relative tolerance to 10000 and surprisingly obtained the same result as before.  I think I am going to recreate our COMSOL model without any unnessesary things implemented and only include a stationary solver, then run this test again as it seems like this should have had some affect on the output.
  41   Fri Jun 28 14:03:07 2013 Deep ChatterjeeGeneralGeneralApproach for TE noise for beam in cavity

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.

Attachment 1: LL_derivations.pdf
LL_derivations.pdf LL_derivations.pdf LL_derivations.pdf LL_derivations.pdf
  53   Wed Jul 10 17:52:18 2013 KojiGeneralGeneralHow to force time steps for a time dependent simulation

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

  54   Thu Jul 11 19:16:12 2013 Deep ChatterjeeGeneralGeneralTE noise problem - Matlab script

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.

Attachment 1: TE_model.m
% TE calculation using 2D axisymmetric model. This script calls the model
% file, extracts grad_T and does the calculations for the different
% frequency
% All parameter are defined in this script and passed to the COMSOL model
% via arguments
% Parameters, Quantitites in [m]
param.beam.radius = 0.05;
param.mirror.radius = 100 * param.beam.radius; 
%Define the radius of the mirror as 100 times that of the beam radius
param.mirror.height = 200 * param.beam.radius;
... 94 more lines ...
Attachment 2: almost_infinite_2D_axisym_TE_4.m
function out = almost_infinite_2D_axisym_TE_4(param)
% Putting all the parameters in the variables
r0 = param.beam.radius;
R = param.mirror.radius;
H = param.mirror.height;
dt = param.COMSOL.simulation_time_step;
fmod = param.COMSOL.fmod;
t_end = param.COMSOL.simulation_end_time;
p0 = param.constant.F0/(pi*r0^2);% The amplitude of Osc. pressure
T_amb = param.constant.T_amb;
... 121 more lines ...
  57   Fri Jul 12 19:43:32 2013 Deep ChatterjeeGeneralGeneralTE noise problem - Matlab script

Quote:

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.

 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.

 

Attachment 1: codes_TE_calc.zip
  58   Sun Jul 14 15:13:46 2013 Deep ChatterjeeGeneralGeneralTE noise problem - Matlab script

Quote:

Quote:

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.

 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.

 

 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.

Attachment 1: codes_TE_calc_(2).zip
  60   Mon Jul 15 16:05:29 2013 Deep ChatterjeeGeneralGeneralTE noise problem - Matlab script

Quote:

Quote:

Quote:

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.

 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.

 

 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.

 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.

Jul_15.bmp

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

Attachment 1: codes_TE_calc.zip
  76   Mon Sep 16 10:07:32 2013 EvanGeneralConfigurationCOMSOL 4.3 on the OS X command line

If you're running Matlab scripts that iterate over simulation parameters (à la Tara), you might find it useful to be able to run everything from the command line (instead of launching the "Comsol with Matlab" GUI application).

First, the comsol command wasn't in my path, so I symlinked it to someplace where bash could find it:

ln -s /Applications/COMSOL43a/bin/comsol /opt/local/bin/comsol

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

comsol server > server_log.txt &

matlab -nodesktop -nosplash -r "run /Applications/COMSOL43a/mli/mphstart; comsol_script; exit"

The first command starts the Comsol server in the background and cats its output into a log file. The second command launches Matlab, runs the initialization script which makes Matlab aware of Comsol, runs my Comsol/Matlab script (comsol_script.m), and then exits.

  82   Thu Apr 17 16:52:09 2014 EvanGeneralCharacterizationInterpreting logfiles and picking a solver

Here are two entries by Walter Frei on the Comsol blog that I've found useful.

Solving Linear Static Finite Element Models: Tells you how to interpret all those numbers that Comsol dumps into its logfiles.

Solutions to Linear Systems of Equations: Direct and Iterative Solvers: Explains a bit more about Comsol's solvers. Apparently, PARDISO is usually fast and SPOOLES is usually slow (the tradeoff is apparently that SPOOLES uses less memory). I've found this to be personally true with at least one model I've worked with.

  99   Tue Sep 30 11:30:27 2014 Nic, Dmass, EvanGeneralConfigurationGravity in Comsol

Here is a set of slides by Yoichi Aso on how to handle gravity in Comsol.

  101   Sat Sep 5 15:17:41 2015 ranaGeneralConfigurationFEA logs merged

I moved the only entry from the 'ENG_FEA' log into the COMSOL log and then renamed that logbook as 'FEA' since we don't need two FEA logs.

Also renamed 'AdhikariLab' log as ATF.

  102   Fri Dec 4 18:32:02 2015 ranaGeneralConfigurationCOMSOL: remote server w/ matlab

This summarizes how to get the remote comsol server to run. COMSOL 5.1.0.234 is now on tegmeni thanks to Larry.

On the server:

rana@tegmeni|~> /usr/local/comsol51/multiphysics/bin/comsol server -login force

This starts up a comsol server instance, listening on port 2036. '-login force' will ask you to supply a username and password which you make up. You will have to enter these later from the client side.

On my laptop, from the MATLAB prompt:

>> addpath('/Applications/COMSOL51/Multiphysics/mli/')

>> mphstart('tegmeni.ligo.caltech.edu', 2036,'uname','pword')

That's it! Now you just run the matlab script which runs the COMSOL file.

I'm attaching a tarball of the .mph file (written by Dmitry Kopstov from MSU) and the associated matlab scripts which change the parameters, as well as looping over test mass thickness to produce a plot of Brownian noise PSD v. 

Attachment 1: BrownianSweep.tgz
  109   Sun May 7 18:22:35 2017 ranaGeneralVoyagerVoyager ITM: Radiative cooling with cold shield and cold CP

I took Aidan's COMSOL model for the ITM from a couple years ago and updated it with some more details:

  1. Through radiative cooling only, the ITM is cooled to 103 K. Taking it to 123 K will be accomplished by adding a ring heater to the ITM.
  2. Assume 3 W of heating from main laser beam onto ITM HR face.
  3. Emissivity of ITM barrel is 0.95. Emissivity of HR* and AR faces is 0.5.
  4. The CP and the Inner Shield are kept fixed at 80 K. This is to simulated the effect of having conductive cooling with cold straps. This needs to be checked in more detail by actually modeling thermal straps.
  5. Emissivity of the CP is 0.9.
  6. The total length of the inner shield is 5 m. The CP is at z = 0 m and the ITM is at z = 2.25 m. We should check what the result would be if the shield is ~1m shorter or longer.

In the attached image, I have made one quadrant of the tubular cryo shield transparent just for clarity - the actual modelled tube is 3 cm thick, made of aluminum, has an emissivity of 0.95 on the inside and 0.03 on the outside (to simulate what we would get from polished aluminum or gold coating).

This files is in our GitLab: https://git.ligo.org/rana-adhikari/CryogenicLIGO/blob/master/FEA/ITM-ColdShield-CP.mph

*I am suspicious of just using a single emissivity number for the AR and HR coatings. Since we are concerned with wavelengths which are long w.r.t. the coating thickness, it may be that the HR and AR coatings have a complicated wavelength dependence in the 5-50 micron band.

Attachment 1: ITM-ColdShield-CP.png
ITM-ColdShield-CP.png
  110   Mon Jul 24 15:35:34 2017 MariiaGeneralConfigurationRunning Comsol to Matlab

WIth Gautam's help, I have created a user directory in 40 meter Lab and copied Rana's documents (MATLAB coating files) from flash card into it. After that, from this elog by Rana : COMSOL: remote server w/ matlab from Fri Dec 4 18:32:02 2015,  ran the matlab document BarrelCoating which resulted in the following error 

Messages:
Could not obtain license for#LiveLink for MATLAB
 
Could not obtain license for LiveLink for MATLAB.
 
License error: -4.
Licensed number of users already reached.
 
Has anybody seen this before?


 

  111   Mon Jul 24 15:54:26 2017 KojiGeneralConfigurationRunning Comsol to Matlab

The number of licenses already used by whom / still remains can be confirmed by running the following command on a comsol-installed linux machine

$ cd /usr/local/comsol51/multiphysics/license/glnxa64
$ ./lmstat -c ../license.dat -a

 

  112   Wed Jul 26 20:14:46 2017 ranaGeneralConfigurationRunning Comsol to Matlab

I've just tried this out on my desktop machine using COMSOL 5.1 and its still working. Which COMSOL is installed on optimus at the 40m ?

  113   Fri Jul 28 15:48:58 2017 MariiaGeneralConfigurationComsol batch for windows

Using the written path from elog by ericq: Computer Scripts/Programs, Comsol can be run from the directory on the distant computer: /cvs/cds/caltech/apps/linux64/comsol51/bin/glnxa64/comsol batch -inputfile Model1.mph -outputfile Model_solv.mph. To transfer files from Linux to Windows : the command pscp.

The method:

1. To download PuTTY and make a coonection with a distant computer.

2. Linux terminal: ssh optimus % allows to go into the whole system

3. Windows terminal: pscp "the path of the file" controls@nodus.ligo.caltech.edu:/users/.../ % allows to transfer mph-file from Windows to Linux

4. Linux terminal: /cvs/cds/caltech/apps/linux64/comsol51/bin/glnxa64/comsol batch -inputfile Model1.mph -outputfile Model_solv.mph % calculation on the distant computer, without outputfile solution would be stored in the input file

5. Windows terminal: pscp controls@nodus.ligo.caltech.edu:/users/.../Model_solv.mph Documents\ % copy file from Linux to Windows, Documents is the name of the folder in Windows 

 

 

  114   Mon Jul 31 22:18:57 2017 ranaGeneralGeneralusing more than 12 cores in matlab

Since 2014, the limit of 12 workers using the matlab parallel computing toolbox has been lifted. Today, I was able to get this to work. There's a trick.

Usually, when you start up matlab and run a parallel thing like 'parfor', it just uses a default profile 'local' which limits you to 12 workers. You can try to ask for more by doing 'parpool(40)' for 40 workers, but it will tell you that NumWorkers = 12 and you're out of luck. So instead:

myCluster = parcluster('local')

myCluster.NumWorkers = 40;

saveProfile(myCluster);

parpool('myCluster', 40)

It seems that it needs the max # of workers and the requested number of workers to be 40 to use 40, otherwise you'll just get 12 (as of matlab 2016a).

Attachment 1: Screen_Shot_2017-07-31_at_10.11.35_PM.png
Screen_Shot_2017-07-31_at_10.11.35_PM.png
  127   Sat Mar 17 15:27:48 2018 ranaGeneralGeneralfile size >> small

When saving your COMSOL files do these two things to make the files much smaller (good for saving in version control and sharing):

  1. File -> Compact History
  2. Preferences -> Files -> Optimize for File Size (not speed)
  128   Mon Aug 20 15:44:56 2018 ranaGeneralGeneralfile size >> small

Also,

  1. click 'Clear Mesh' under the mesh menu
  2. 'Clear Solutions' under the Study menu

In this way the file sizes will be ~100 kB instead of 10's of MB.

Quote:

When saving your COMSOL files do these two things to make the files much smaller (good for saving in version control and sharing):

  1. File -> Compact History
  2. Preferences -> Files -> Optimize for File Size (not speed)

 

  130   Sun Aug 26 19:21:27 2018 ranaGeneralVoyagerVoyager ITM: Radiative cooling with cold shield and cold CP

this is a time dependent model of the previous steady-state one

  • Cold Shield and CP held at a constant 60 K
  • 3 W heat input to the ITM from the main laser beam
  • radiative cooling to the shield
  • ITM barrel emissivity = 0.9
  • ITM HR/AR emissivity = 0.5/0.5

So the cooldown time w/o a heat switch is ~4 days. Since this is less than the usual pumpdown time required to open the gate valves on the beamtubes, perhaps no heatswitch or in-vac cryogens are required.

Attachment 1: ITM-Cooldown.pdf
ITM-Cooldown.pdf
Attachment 2: CoolDown.webm
  138   Tue May 12 14:16:28 2020 KojiGeneralGeneralFEA tutorial resources

cf. Forwarded email from Stephen

1) Tuesday Demo - Basics of FEA Meshing G2000696
2) CIT SYS User Guides, How to Use the FEA User Group T2000295
3) CIT SYS User Guides, How to Use the ANSYS Learning Hub T2000236

Fabrice's SAMS piezo actuator second prototype E1900383

  139   Tue Aug 11 11:16:29 2020 aaronGeneralConfigurationCOMSOL with Matlab without display

When running comsol with matlab interface on sandbox1, it is usually most convenient to ssh with screen forwarding (eg '-CY') and launch COMSOL with matlab by following the instructions in the livelink manual.

Sometimes, it is necessary to run COMSOL without any display available. In that case, the Instructions in the manual are a little unclear. Here are the detailed steps that let me run my script '/home/amarkowi/metamaterials/run_spiral.m' with no screen forwarding.

1. ssh onto sandbox1 by entering the following at your laptop command prompt
Aaron’s-latpop $ ssh aaron.markowitz@sandbox1.ligo.caltech.edu
                           password: [enter sandbox1 credentials]
 
2. start a tmux session for starting the COMSOL server
aaron.markowitz@sandbox1:~$ tmux
 
3. start a comsol server from the tmux session
aaron.markowitz@sandbox1:~$ comsol54 mphserver -login force -port 2020
                                                    Username: whoever
                                                    Password: whatever
                                                    Confirm password: whatever
 
4. detach the tmux session by pressing ‘ctrl-b’ followed by ‘d’
 
5. (optionally, you can start a new tmux session for your matlab work by running tmux again at your main sandbox prompt)
 
6. Start matlab by running
>> matlab -nodesktop -mlnosplash
 
6. Add the comsol directory to the matlab path by running at the matlab prompt
>> addpath(‘/localhome/comsol54/multiphysics/mli/‘)
 
6. Start the matlab with comsol interface by running the following at the matlab prompt
>>   mphstart(’sandbox1.ligo.caltech.edu’, 2020, ‘whoever’, ‘whatever’)
 
7. change into the directory containing the script, and run it
>> cd /home/amarkowi/metamaterials
>> run_spiral
  140   Tue Aug 11 16:35:07 2020 aaronGeneralConfigurationCOMSOL: remote server w/ matlab

To run COMSOL on sandbox1 with no graphical Interface, here are the steps that worked for me (Tue Aug 11 16:35:51 2020) from a mac on the Caltech VPN.

1. ssh onto sandbox1 with screen forwarding (-Y). Make sure you have a compatible version of XQuartz or a substitute. -C specifies data compression, which may be useful on slow connections
Aaron’s-laptop $ ssh -CY aaron.markowitz@sandbox1.ligo.caltech.edu
                           password: [enter sandbox1 credentials]
 
2. Launch Matlab with COMSOL 5.4, specifying no graphical interface, by running
aaron.markowitz comsol54 mphserver matlab -nodesktop -mlnosplash script_name’

3. If there is still a splash screen from the COMSOL server, you will have to specify nosplash by adding the following line to your .bash_profile (in your home directory)

export COMSOL_MATLAB_INIT=’matlab -nosplash’

 
4. You can run whatever comsol script you need. Make sure that in your script you import the comsol functions by calling the following
import com.comsol.model.*
import com.comsol.model.util.*
  100   Sat Sep 5 15:04:43 2015 Dennis CoyneMechanicsConfigurationsummary of FEA modal model to State Space model

At the 2014 commissioning workshop, I presented a summary of my efforts in converting finite element modal models into state space models:

https://dcc.ligo.org/LIGO-G1400099

I also provided the attached written summary for a report on the workshop that Aidan Brooks and Gabriele Vajente are preparing.

Attachment 1: FEM2SS_summary.docx
  105   Tue Jun 28 15:50:45 2016 Joy WestlandMechanicsAnalysisA Simple Model of the Modal Analysis of a Cantilever Circular/Cylinder Cross Section in ANSYS tutorial

Here is a tutorial to implement a simple Modal Analysis of a Cantilever Cylinder

1. Open the ANSYS workbench

2. Drop and drag the "Modal" analysis system into the project schematic

3. Right click on "Engineering Data" and edit the material. There are predetermined properties found in the "Engineering Data Sources" (Right click on the description or press the books in the top right hand corner). For this tutorial, fused silica was created with these parameters:

  • Young's Modulus E: 7.2E10 Nm^-2 (Isotropic Elasticity)
  • Density p: 2200 kgm^-3 (Density)
  • Specific Heat c: 770 J/kgK (Specific Heat)
  • Isotropic Thermal Conductivity k: 1.38 W/mk (Isotropic Thermal Conductivity)
  • Coefficient of Thermal Expansion: 3.9E-7 (Isotropic Secant Coefficient of Thermal Expansion)
  • Reference Temperature: 21 C (Isotropic Secant Coefficient of Thermal Expansion)
  • Poisson's Ratio: 0.17 (Isotropic Elasticity)

In order to put in the parameters, just drag and drop from the toolbox. Also make sure the filter is turned off (that's the filter symbol next to the engineering data sources/books in the top right hand corner) this will allow a user to few all the parameters.

4. Return to the project, and open/edit the geometry. The cantilever that was used for this tutorial is a cylinder with a diameter of 2 mm with a length of 200 mm, shown in image 1.

  • How to Create a geometry:
    • Open up the geometry
    • Decide which plane you want to start drawing in: XY Plane, ZX Plane, YZ Plane
    • Click on "new sketch" in the 4th row from the top after "Sketch1" it's blue sketch
    • Whatever plane you decide to choose, press the "Look at Face/Plane/Sketch" icon in the top right hand corner (The person looking at a face)
    • This will allow the view to go to the view where you want to sketch on
    •  Click on "Sketching" which is next to the "Modeling" right above the "Details View"
    • Here will be Several Options
      • Draw: This is where you can decide what shape you want to use
      • Modify: This is the extras such as fillets, chamfers and corners
      • Dimensions: This is where you can select what type of dimensions you will be using
      • Constraints: This is where you can constrain lines or connect lines
        • When you want to connect lines together, click on each end and make them "Coincident" this way the shape will connect when you change dimensions
      • Settings: This is where you can show a grid or snap grid
    • Sketch Color:
      • Teal: Underdefined: The sketch does not have enough dimensions to constrain and tell the program what exactly the shape is
      • Blue: Defined: The sketch is fully defined and does not need any more dimensions
      • Red: Overdefined: The sketch has repeating dimensions that causes errors
    • For this tutorial, draw a circle. I made mine coincident to the middle by pressing the vertex in the middle and one of the axis and then repeat the process with the second axis to center the circle
    • Dimensions: The diameter is 2 mm, under dimensions to display the number, click "Display" and then "Value" instead of "name"
    • To change the dimensions, the "Details View" in the bottom left hand corner will say D1 and from there you can change the diameter
    • From there in order to create the length of the cylinder, press "Extrude" at the top (4th row) towards the right. From there it will ask what geometry, click on your sketch and click "Apply". Make sure the "Operation" is "Add Material" to Extrude.
    • Then choose a "Depth" in this case it's 200 mm. Then click "Generate" the lightning bolt next to the "New Sketch" symbol.
    • Your cylinder should look like image 1.

5. Return to the project and open the Model tab. With that open, go to the subtab of "Geometry", named "Solid". Under Solid go to "Assignment" under "Material" to assign Fused Silica or whatever material you decide to choose.

6. Go to the "Mesh" tab. In this case the "Element Size" under the "Sizing" tab is 0.04 and the "Relevance Center" as "Fine". Right click on the mesh and press "Generate mesh".

7. Under the Modal tab, go to the Analysis Settings and change the "Max Modes to Find" for ANSYS to calculate. In this tutorial, the amount of modes that were used was 17.

8. Right Click the Modal and press on "Fixed Support". This will make the bar a cantilever bar, once the geometry is set to fix one face of the bar. Press on one of the faces and click on "apply" on the Geometry.

9. Right click on the solution and press "Solve". Let ANSYS run the modes through.

10. On the Solution tab, the "Tabular Data" is listed but the Total Deformation has not been listed. In order to do so, select all the frequency of the modes in the column and right click and press "Create Mode Shape Results".

11. Once loaded, the total deformation will have a lightning bolt next to each entry, right click and press "Solve" or "Evaluate All Results". This will make all the entries have a green check mark.

12. This will give the user the ability to animate each entry. This is down by clicking the play (sideways triangle button). ANSYS will run through the simulation for that mode that is selected, shown in image 2.

13. Using Mathematica (or another computational system) input the analytical solution for a Cantilever bar fixed to one end.

  • wn=(knL)2Sqrt(EI/mL4)
    • Where:
    • wn is the frequency measured in radians/second
    • E = Young's Modulus
    • I = moment of cross section
      • Irectangle = ba3/12 (b and a are the sides)
      • Icircle = (Pi/64)(d4) (d is the diameter)
    • m = mass per unit length
    • L = length of the bar
    • kn relates to the amount of nodes
      • k1L = 1.875, k2L = 4.694, k3L = 7.855, k4L = 10.996, k5L = 14.137
      • n greater than 5: knL = (2n -1)(Pi/2)
  • Convert wn to f measured in Hz

14. Compare only to the bar actually bending, not twisting or contracting. There are modes that are the same due to the symmetry of the bar. In image 3 the underlined frequencies compare to the analytical calculations (Mathematica) and the computational calculations (ANSYS).

Attachment 1: Simple_Cantilever_200_mm_Long_Geometry.PNG
Simple_Cantilever_200_mm_Long_Geometry.PNG
Attachment 2: Simple_Cantilever_200_mm_Long.PNG
Simple_Cantilever_200_mm_Long.PNG
Attachment 3: Comparison_between_ANSYS_and_Calculations_Fused_Silica_200_mm_L.PNG
Comparison_between_ANSYS_and_Calculations_Fused_Silica_200_mm_L.PNG
  106   Tue Jun 28 16:54:23 2016 Nikhil MathurMechanicsAnalysisAnsys 14.5 Indroductory Tutorial: Modal Frequency Convergence

This tutorial will go through how to show frequency convergence for a cylindrical cantilever using Ansys 14.5 and Mathematica.

Attachment 1: ConvergenceTutorial.pdf.gz
  107   Wed Jun 29 14:14:44 2016 Joy WestlandMechanicsAnalysisANSYS Tutorials with Basic Meshing

Here are a series of tutorials for basic meshing principles from ANSYS Meshing Basics:

  1. ANSYS Meshing Fine Mesh Basic Tutorial 1: https://www.youtube.com/watch?v=sZIX3CJkWBE
  2. ANSYS Meshing Method Basic Tutorial 2: https://www.youtube.com/watch?v=ATSnvFUbEk4
  3. ANSYS Meshing Refinement Mesh Basic Tutorial 3: https://www.youtube.com/watch?v=QVfJiydCu2g
  4. ANSYS Meshing Sizing Mesh Basic Tutorial 4: https://www.youtube.com/watch?v=ODJ11YDC6ac
  5. ANSYS Meshing Spheres of Influence Basic Tutorial 5: https://www.youtube.com/watch?v=AQlGccK_X-0
  6. Body of Influnce 6: https://www.youtube.com/watch?v=3E1p1w32jt0
  7. Section Plane 7: https://www.youtube.com/watch?v=hhZOJSESJEU
  8. Mapped and Match Control 8: https://www.youtube.com/watch?v=3pogywR8Vgw
  9. MultiZone and Inflation 9: https://www.youtube.com/watch?v=vxVBvANHuB0
  10. Sizing Soft/Hard 10: https://www.youtube.com/watch?v=D8UzdC5Gmhs
  108   Fri Jul 29 14:33:41 2016 Joy WestlandMechanicsGeneralA Tutorial in Importing SolidWorks Files and Applying a Gaussian Force in ANSYS for a LIGO Test Mass

Here is a tutorial in importing SolidWorks into ANSYS and the steps needed to apply a Gaussian Force to the LIGO test mass that’s imported.

Using SolidWorks:

  1. Download the SolidWorks Zip Folders from the DCC
    1. https://dcc.ligo.org/login/index.shtml?entityID=https%3A%2F%2Fdcc.ligo.org%2Fshibboleth-sp&return=https%3A%2F%2Fdcc.ligo.org%2FShibboleth.sso%2FLogin%3FSAMLDS%3D1%26target%3Dss%253Amem%253A15b6c314d87e3fa8b3768d89cb6b9836fe39c754 (LIGO-D1000760-v4)
    2. Open/export the zip folders into a different file
  2. For only the test mass and ears, open the file called: D0902456 ITM OPTIC WITH EARS ASSEMBLY. This will represent the bottom mass.
  3. An assembly can be imported into ANSYS but it’s easier to convert an assembly into a part before importing it into ANSYS. Also this will allow you to perform a split line.
    1. To do this, resave the assembly but before pressing save, change the file to PRT extension instead of Assembly.
  4. Open the PRT file of the test mass, make a new sketch on the face of the test mass.
  5. Draw a circle (Or whatever shape) on the test mass. In this case a 0.1 m circle was made in the middle of the test mass.
  6. From there go to “Insert” and then press “Curve” and then go to “Split Line”. Split line allows a user to project the sketch onto the part. This makes it possible to press the sketch separately from the part. This is a useful tool because in ANSYS, once split line has been used on a sketch, a user can press that shape separately on the face of the object. The importance of this tool allows a user to click on the circle separately and apply a force in that particular area.
    1. In “Split Line” press the “sketch” you want to project onto the “face” of the object.

Importing SolidWorks into ANSYS:

  1. Once the part has been made, choose the “static structural analysis”.
  2. In the geometry: Import the test mass that was converted into a part
  3. Once imported, press the geometry again and “edit”
  4. In the edit module, “generate” the part and once the part has loaded notice the circle can be pressed such that it is a different part of the test mass surface.
  5. Exit geometry module

Applying Gaussian Force:

  1. Enter the “model” module to edit the setup of the analysis
  2. Assign the material that you want for the test mass under “Geometry” à “Solid” under the “Material” and “Assignment”
  3. Under the “Coordinate Systems”, make a new coordinate system
  4. “Coordinate System:
    1. Click on the face, in this case the circle made on the test mass for the “geometry” and “apply”
    2. Under “Definition” à “Type” à Change from Cartesian to Cylindrical. This allows the Gaussian Force to be distributed correctly
  5. Under “Analysis Settings” Turn on “Large Deflection”
  6. Add “gravity” under “Static Structural”
    1. Treat the test mass is if it were the bottom mass so make sure gravity is pointing the correct way, such that the wires on the ears would extend upward.
  7. Insert “Fix” to the top surfaces of the ears
  8. Insert “Pressure”
    1. Press on the circle made to apply for the pressure
    2. Under “Definition” change the “Magnitude” to “Function” (the arrow at the end of the Magnitude entry bar)
    3. Once function has been activated: A Gaussian force of: 1/((3.141592654*0.0156^2)*2.718281828^((x/0.0156)^2)) was used in SI units. Note that ANSYS does not use symbols. Once that’s entered into the “Magnitude”, under “Function” the “Coordinate System” will appear. Change that to the coordinate system that was made in step 4 of the cylindrical system.
    4. Under “Graph Controls” Make sure that the X-Axis is changed from Time to X.
    5. Pick a range for the graph and the number of segments that you want to look at.
  9. “Solve” the system and right click on “Solution” and evaluate the different results you need.
Attachment 1: Bottom_Test_Mass_With_Ears.PNG
Bottom_Test_Mass_With_Ears.PNG
Attachment 2: Applying_Parameters.PNG
Applying_Parameters.PNG
Attachment 3: Applying_Gaussian_Force.PNG
Applying_Gaussian_Force.PNG
  115   Thu Nov 2 17:23:56 2017 AaronMechanicsPonderSqueezeModelling suspension noise

aLIGO Suspensions Toy Model

On Wednesday I started making my own model of the aLIGO suspensions, with the top of the silica fibers attached to ears that are fixed rather attached to an additional suspension stage (so this will be a one stage suspension).

I grabbed the aLIGO ear design from the DCC: LIGO-D080751-v4

I am almost done with the model, should have it working tomorrow and will add it to the experimental gravity github in an appropriate place.

  116   Fri Nov 3 15:03:10 2017 AaronMechanicsPonderSqueezeModelling suspension noise

Model Geometry

Test Mass

I found the dimension of the test mass flat in the drawings of the mock test mass design here: LIGO-D080687.

Fibers

I modelled the fibers with the profile described in LIGO-D080751, fig 3.7.

Ears

I grabbed the d-values from LIGO-T1000545, but since the d-value is defined as the distance from the center of mass (of the penultimate mass (PUM) or the E/ITM) to the bend point (BP) of the fiber (I believe the point on the fiber with maximal flexure in the fundamental mode), I did not go through the effort of figuring out where the bend point is but rather grabbed the horn-CM distance from LIGO-T1100407

I wanted to get the real aLIGO parameters for the first version of this model, and have parametrized the model in such a way that I can define all of the parameters that need to change (surface area of the ear-TM bond, length of the fibers, thickness and profile of fibers, d-value, etc) and scale them with mass in some way for future iterations on this design.

I need to pare down the number of parameters, because I started by fully defining the ears and now am importing a 3D model of the ears and planning to scale these with mass.

Materials

For the material of the entire test mass and suspension, I used the fused silica that is specified as [solid,NIST SRM 739 - Type I]. I wasn't sure the difference between the types of silica, but this one said SRM so I thought it might have been defined on my distribution of COMSOL by a LIGO person. A quick google search showed me that person may have been rana?? https://labcit.ligo.caltech.edu/~rana/research/etm.html

Physics

I'm using a solid mechanics model.

Fixed Boundary Constraint

I fixed the position of the bonding surfaces for the PUM ears, so it is as if they are contacted to a completely fixed PUM (the PUM is not included in the model, but the upper ears are included, so the constraint is on the ears not the fiber. See drawing).

Gravity

I added gravity to all parts of the model. Apparently, it is not trivial to use gravity in a frequency domain study in COMSOL, as described in this presentation here. Fortunately, the presentation in the link is interested in the transfer function for a mass on a string also, so I follow the simulation steps they describe below.

Boundary Load

I add a boundary load that will vary sinusoidally for the frequency domain study.

Mesh

I have not yet messed with the meshing for these models. Obviously the points with more flexture and smaller parts (like at the horns of the ears, the tapering parts of the fibers, etc) will require a finer mesh.

Study

I need to incorporate the advice on how to build this study described in the link above. The following might also be useful, though I haven't looked through them yet:

https://www.comsol.com/model/dynamics-of-double-pendulum-14021

https://www.comsol.se/forum/thread/4843/pendulum-response?last=2010-04-27T01:48:26Z

  117   Wed Nov 15 14:05:12 2017 AaronMechanicsPonderSqueezeModelling suspension noise

Model Geometry

I pared down the number of parameters in the model to only the necessary ones. These are the ones that remain:
 
TM_radius: Radius of the test mass
TM_width: Width of the test mass
TM_flats: length of TM flats
ear_length: length of the ear
horn_spacing: length of the ear
horn_gap: gap between the top of the horn and the TM on the near side
d_val: distance from the CM to the bend point
horn_BP: distance from the horn to the bend point
ear_height_tot_nominal: nominal total height of the ear and horn for the unscaled (aLIGO) design (this name made more sense in a previous version of the model)
fiber_stock_length: length of the fiber's stock
fiber_neck_length: length of the fiber's neck
fiber_thick_length: length of the thick section of the fiber
fiber_main_length: length of the main section of fiber (the thinnest part)
fiber_taper_length: length of the tapering section of fiber
fiber_stock_radius: radius of the fiber stock
fiber_thick_radius: radius of the thick section of fiber
fiber_main_radius: radius of the main section of fiber
F_load: force of the boundary load used for the excitation in the frequency domain study
ear_scale_height: scale the height of the ear
ear_scale_length: scale the length of the ear
ear_scale_width: scale the width of the ear

Materials

For the material of the entire test mass and suspension, I used the fused silica that is specified as [solid,NIST SRM 739 - Type I]. I wasn't sure the difference between the types of silica, but this one said SRM so I thought it might have been defined on my distribution of COMSOL by a LIGO person. A quick google search showed me that person may have been rana?? https://labcit.ligo.caltech.edu/~rana/research/etm.html

Physics

Rana suggests that for the purpose of this study, it is not necessary to actually have COMSOL handle gravity as a restoring force... I'm not sure if I understand why this is yet? It seems that if we are interested in the relative strain energy in different parts of the wire compared to other parts of the system, it is important that the wire be under tension. If we have no gravity, the wire is effectively not under tension.

Mesh

I have not yet messed with the meshing for these models. Obviously the points with more flexture and smaller parts (like at the horns of the ears, the tapering parts of the fibers, etc) will require a finer mesh.

Study

I need to incorporate the advice on how to build this study described in the link above. The following might also be useful, though I haven't looked through them yet:

https://www.comsol.com/model/dynamics-of-double-pendulum-14021

https://www.comsol.se/forum/thread/4843/pendulum-response?last=2010-04-27T01:48:26Z

  118   Mon Dec 4 16:27:13 2017 aaronMechanicsPonderSqueeze 

Meshing Surface Layers

Defining New Selections

I don't know why I wasn't seeing this problem with previous models (perhaps because I wasn't importing any CAD or STEP files), but my latest attempts at meshing and selecting specific domains of my model were being thwarted by inconsistent domain definitions. I was previously always manually selecting domains, which is confusing because all domains just get assigned a number when they are created. Worse, sometimes the numbers assigned to domains change when the model rebuilds, especially if there has been a significant change in the model geometry. This results in later steps selecting the wrong set of domains (or boundaries, etc).

To fix this problem, I created new selections (sets of domains, boundaries, or etc that receive their own label and can be selected as a group later in the model). The new selections include:

-Domain selections that separate surface and bulk layer for all domains in the fibers (fiber stock, neck, thick section, taper, and main section, in order from the horn to the center of the fiber)
-Boundary selections for all domain selections described above

This is probably also a necessary step for getting reliable results when interfacing with MATLAB, and might explain some weird problems I was running in to a while ago and just made haphazard fixes for.

Mesh Steps

I use the following mesh steps to get what seems like a pretty reliable meshing:

  1. Mesh the upper tapers (bulk and surface separately) with a free tetrahedral mesh
    1. I mostly use the defaults for an extremely coarse mesh, but the only parameter that seems to make a large difference is the minimum mesh size. I set this to the skin_depth for the surface layer, and (fiber_main_radius-skin_depth) for the bulk. fiber_main_radius-skin_depth should be the radius of the bulk domain, and the skin_depth characterizes the smallest length scale in the surface layers. I have some limited ability to tweak the minimum mesh sizes when the surface layer is comparable in size to the fiber_main_radius (so the surface layer is comparable to the entire fiber radius), but it seems to be best to keep the mesh this fine when the surface layer becomes small.
  2. Mesh the main fiber (bulk and surface separately) with a swept mesh
    1. I create a distribution with a fixed number of elements, at ceil(fiber_main_length/fiber_stock_radius/2). This is somewhat arbitrary--the stock has the largest radial length scale, so I figured I'd divide up the main fiber into units that tall. A better thing would be to know how high a mode we are interested in studying, and break up the fiber into enough pieces to observe that mode. Seems fine for now, might want this distribution to be a bit coarser though.
  3. Mesh the lower tapers (bulk and surface separately) with a free tetrahedral mesh
    1. Use the same size settings as on the upper tapers
  4. Mesh the thick sections (bulk and surface separately) with a swept mesh
    1. The distribution uses a fixed number of ceil(fiber_thick_length/fiber_stock_radius/2) elements. Again perhaps this can be coarser; it also doesn't attempt to make a finer mesh at the thermoelastic cancellation region.
  5. Mesh the necks (bulk and surface separately) with a free tetrahedral mesh.
    1. I use an extremely coarse mesh with the minimum mesh element size set to fiber_thick_radius-skin_depth for the bulk; for the surface it is an extra coarse mesh instead, because the extremely coarse mesh gave low quality mesh elements. I'm not sure why this is, I can't see much difference and the problem only seems to happen to one of the 8 neck sections (why not in all sections?). The problem only arises when the skin depth is less than 20um, for the other parameters at the nominal aLIGO values.
  6. Mesh the stocks (bulk and surface separately) with a swept mesh
    1. Distribution is again ceil(fiber_stock_length/fiber_stock_radius/2), resulting in 2 vertical divisions of the stock.
  7. Mesh the horns with a free tetrahedral mesh
    1. Size is an extremely coarse mesh, where the minimum element size is set to the skin_depth (because the horn sees the boundary of the stock surface, so its smallest adjacent element can have a scale down to the skin_depth), and maximum growth rate is increased to 9 (any higher results in low quality mesh elements; I set it high because most of the horn does not need to be meshed as finely as the part directly in contact with the fiber stock).

Note on the skin depth: I defined the thickness of the surface layer by the parameter skin_depth. Since we are interested in how the energy in the surface layer changes with the radius of the fiber, and expect that the surface depth doesn't change much with fiber radius (a very thick fiber will have a few micron surface layer; so will a fiber half that diameter). I managed to get the mesh working with a 15um surface layer, but was having trouble getting a 10um layer. I wasn't completely sure how small a surface layer would be necessary, but for a test mass 1/100 the size of the aLIGO masses, the main part of the fiber would have a radius of 20um, so if we want to study the surface layer for masses down to that size I figure the skin depth should be at most around 10-15um. It would be better to be motivated by the actual scale of the physics going on at the surface--how deep to micro cracks and other lossy imperfections go?

Study

I'm running the frequency study on a small number (2) of frequencies in each decade from 60Hz to 10kHz (so there should be 5-8 total frequencies in the study). I started it at 4:25pm on my local machine, and it hasn't gotten very far in 30 minutes, so I may abort the study and try to make the mesh coarser--especially the distribution settings, which are easy to change.

  119   Mon Dec 4 17:42:53 2017 gautamMechanicsPonderSqueezeFEA on optimus

We could run the simulations on the 32 core machine in the 40m lab (optimus)? I think Mariia was running some of her studies on optimus, and even though we had some problems with the licensing initially, I think she resolved these and has detailed the procedure in her elogs...

Quote:

 

Study

I'm running the frequency study on a small number (2) of frequencies in each decade from 60Hz to 10kHz (so there should be 5-8 total frequencies in the study). I started it at 4:25pm on my local machine, and it hasn't gotten very far in 30 minutes, so I may abort the study and try to make the mesh coarser--especially the distribution settings, which are easy to change.

 

  120   Mon Dec 4 19:49:32 2017 gautamMechanicsPonderSqueezeFEA on optimus

It's a good idea, I'll check out her elogs and get it started tonight.

I found that I had the relative tolerance set too low (0.001, while the iterations were not converging further than 0.02 or so); I changed it to 0.1, which let me run over a few modes relatively quickly once I reduced the number of sections on the main part of the fiber to 10. This is not sufficient--at 600Hz the fiber looks completely jagged because the mesh is not fine enough, so using optimus will be necessary. 

Quote:

We could run the simulations on the 32 core machine in the 40m lab (optimus)? I think Mariia was running some of her studies on optimus, and even though we had some problems with the licensing initially, I think she resolved these and has detailed the procedure in her elogs...

 

  121   Tue Dec 5 10:50:54 2017 aaronMechanicsPonderSqueezeFEA on optimus

I had some trouble running this on optimus.

Optimus has COMSOL 5.1 installed, but I made these files in 5.3. I downloaded the comsol 5.3 dvd.iso file last night, but on install I'm now getting the error "No locks available." I wasn't sure if this is a file permissions issue (sometimes the file has been 'locked' when it is saved and I open it in the GUI), an svn issue (most of the information I find online is from people running into some problem with lockd or their NFS), or a licensing issue (I don't think this is likely, COMSOL on my laptop shows more licenses available).

Has anyone seen this? Do I need to do something else with /cvs after install?

  122   Tue Dec 5 19:50:47 2017 aaronMechanicsPonderSqueezeFEA on optimus

Gautam advised me against trying to install version 5.3, lest it break version 5.1--I had already gone through the install, but looking at the install manual it says it shouldn't affect previous installs except that the default behavior when double clicking .mph files will always choose the latest version of COMSOL. Since we mostly run on the command line we should be fine. That said I haven't tested files with COMSOL 5.1.

Zach pointed me to sandbox1.ligo.caltech.edu at the lunch meeting, as well as these notes from Rana about running comsol on a remote server. I couldn't run the file myself on sandbox1 because I don't yet have a home directory there, but I asked Larry if he could set one up for me so I can use the clusters. Gautam helped me run them on sandbox1 by having me scp the file onto the shared 40m drive, then he moved it in to his home directory on sandbox1 and ran it. A version of the model asking for an analysis at 60Hz and 600Hz ran in 10-20 min, and looked good (I was able to scp the output to my laptop and open it in the gui, though I have a script not quite to the point where Matlab can do this for me, it's not ready yet and anyway I wanted to directly see that the output was normal). I modified the study to take 100 steps between 60Hz and 5kHz, and Gautam has now started that run, which if it scales linearly (hopefully better than linear, since it won't have to remesh and etc) will take about 8-16 hours.

Thank you Gautam for the extended mattermost session helping me run these!

ELOG V3.1.3-