ID 
Date 
Author 
Type 
Category 
Subject 
15

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

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

Attachment 1: Iron_Bar_test1_model.m

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

93

Thu Jul 10 16:51:14 2014 
Sam Moore  Optics  General  Duan and Heinert Comparison  (See Plots in attached document)
My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V). I am trying to match Duan's analytical expression with Heinert's analytical expression. This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. ) Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".
It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation. The calculation in Mathematica was very slow, so the upper limits on the integral were truncated. The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steadystate method for Duan's numerical case (Heinert's plot runs much faster). 
Attachment 1: 7_10_14.pdf


96

Mon Jul 14 19:14:31 2014 
Sam Moore  Optics  General  Duan and Heinert Comparison 
Quote: 
(See Plots in attached document)
My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V). I am trying to match Duan's analytical expression with Heinert's analytical expression. This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. ) Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".
It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation. The calculation in Mathematica was very slow, so the upper limits on the integral were truncated. The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steadystate method for Duan's numerical case (Heinert's plot runs much faster).

I have now plotted the DuanHeinert comparison for the case of an infinite upper bound. It turns out that the curves differ by a maximum of 10 percent for low frequencies. Such a discrepancy has been attributed to lack of experimental investigation into this regime (according to Duan). For our purposes, such a discrepancy is acceptable. We will therefore use the Heinert curve for subsequent calculations due to its faster computation time.

Attachment 1: duan_heinert_comparisonInfiniteepsconvertedto.pdf


101

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

Mon Dec 4 17:42:53 2017 
gautam  Mechanics  PonderSqueeze  FEA 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 58 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 coarserespecially the distribution settings, which are easy to change.


120

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

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

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

Tue May 12 14:16:28 2020 
Koji  General  General  FEA 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 
97

Thu Jul 31 20:55:38 2014 
Sam Moore  Optics  General  Finding the Right Meshing for the TIR cavity  In this document, I try to identify I good mesh by comparing the numerical solution from that mesh with my analytical model. Since there are problems with carrying out the analytical calculation, it is still not entirely clear which mesh should be used.

Attachment 1: 7_30_14.pdf


98

Sat Aug 2 00:22:34 2014 
Sam Moore  Optics  General  Finding the Right Meshing for the TIR cavity 
Quote: 
In this document, I try to identify I good mesh by comparing the numerical solution from that mesh with my analytical model. Since there are problems with carrying out the analytical calculation, it is still not entirely clear which mesh should be used.

I have refined the analytical calculation procedure, as outlined in this new document. The procedure indicates that the discrepancy between the analytical and numerical solutions are more likely attributed to meshing inaccuracies. 
Attachment 1: 8_1_14.pdf


141

Thu Jun 22 14:12:04 2023 
Raj  General  aLIGO  Finesse Modelling for aLIGO  The aLIGO model available on IFOSim uses the first surface of the mirror as Port 1 and the other surface as Port 2. This requires us to use negative radii of curvatures among others to get the correct modeling. However, this makes adding surface maps more difficult and was causing issues with getting the simulations to follow what is expected. The surface_map attribute for mirrors doesn't allow us to specify the mirror port in Finesse.
Following a discussion with Kevin, it seems that this approach is being deprecated. His suggestion is to use node 1 as the front of a mirror and node 2 as the back of the mirror. And likewise nodes 1 and 2 as the front of a beamsplitter and nodes 3 and 4 as the back.
I am rewriting the IFOSim aLIGO model with this convention to allow for later simulations not to be plagued by issues and allow for easier simulations of the beam. 
3

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

74

Wed Jul 31 15:39:11 2013 
Deep Chatterjee  Optics  General  First try with paramter optimization for TE and TR noise profiles  After the simulations have been found to match to a fair extent with the analytic results by Heinert and Liu and Thorne, the attempt is check out the parameters for which the TE and TR noise are
close to each other. This was done with the analytic results. The frequency range 10  1000 Hz was looked into. In between this range, the quantity that was minimized was the absolite value
of the logarithm of the ratio between the TR and TE noise. The fminsearch function was used to minimize the mentioned quantity. The parameters that were changed were  conductivity, thermo
optic coefficient and coefficient of linear expansion. The reason for choosing these three were 
> TR noise is independent of coefficient of linear expansion
> TE noise is independent of thermo optic coefficient
> The power law dependence on conductivity is different for the TE and TR cases as can be seen from the analytic expressions
once the code returned the optimized parameters, these values were plugged in and the results were plotted.
**Note that the minimization was done for frequencies between 10 to 1000 Hz 
Attachment 1: Jul31_param_opt.pdf


83

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

84

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

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

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

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

99

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

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

Attachment 1: threeD_cylinderTRnoise.png


Attachment 2: oneDcylinderTRnoise_copy.png


92

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

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

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


Attachment 2: heinert_analyticTest_residual_oneD.eps


53

Wed Jul 10 17:52:18 2013 
Koji  General  General  How 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 "TimeDependent Solver" in Model Builder
 You get bunch of options in TimeDependent 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.) 
21

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


82

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

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

Sun Jun 23 14:00:05 2013 
Emory Brown  General  General  Manipulating 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.tuegriddoc.unituebingen.de/dokuwiki/lib/exe/fetch.php?id=hpcuni%3Asoftwaredocs%3Acomsol%3Acomsol&cache=cache&media=hpcuni:softwaredocs: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. 
10

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

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

48

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

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

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

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

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

Attachment 1: SimpleCylinder.mph

51

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

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


Attachment 1: Eigenfrequency.png


37

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

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

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

Umax (*10^10 J)

1224

1.48859

7011 (used above)

1.51589

13306

1.51891

24181

1.52169

61772

1.52281

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

Thu Nov 2 17:23:56 2017 
Aaron  Mechanics  PonderSqueeze  Modelling 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: LIGOD080751v4
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 
Aaron  Mechanics  PonderSqueeze  Modelling 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: LIGOD080687.
Fibers
I modelled the fibers with the profile described in LIGOD080751, fig 3.7.
Ears
I grabbed the dvalues from LIGOT1000545, but since the dvalue 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 hornCM distance from LIGOT1100407
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 earTM bond, length of the fibers, thickness and profile of fibers, dvalue, 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/dynamicsofdoublependulum14021
https://www.comsol.se/forum/thread/4843/pendulumresponse?last=20100427T01:48:26Z 
117

Wed Nov 15 14:05:12 2017 
Aaron  Mechanics  PonderSqueeze  Modelling 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/dynamicsofdoublependulum14021
https://www.comsol.se/forum/thread/4843/pendulumresponse?last=20100427T01:48:26Z 
80

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

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

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

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

Attachment 1: aluminum2.png


Attachment 2: aluminum3.png


38

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

Attachment 1: ITMPlusRingPlusPID.m

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

Attachment 2: PID.m

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

43

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

20

Mon Jun 17 11:06:17 2013 
Emory Brown  General  General  Plans for the week of 61713 through 62213  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. 
17

Tue Jun 11 10:21:05 2013 
Deep Chatterjee  General  General  Re: 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.TT_amb');
model.param.descr('dT', '');
An empty quote in the description implies that no description is given. The above statements define the parameters mentioned alongside them according to Heinerts paper.
> [line 5361]. the string 'var1' is used to tag all the global variables that are created. variables are expressions created out of the parameters.
Just like the above case, 'set()' is used to create variables and give them an expression, here as one can see 4 variables are created T_amb, beam_size, beam_power. and f_mod.
model.variable.create('var1');
model.variable('var1').set('T_amb', [num2str(param.material.temperature) '[K]']);
model.variable('var1').descr('T_amb', '');
model.variable('var1').set('beam_size', [num2str(param.beam.radius) '[m]']);
model.variable('var1').descr('beam_size', '');
model.variable('var1').set('beam_power', [num2str(param.COMSOL.beam_power) '[W]']);
model.variable('var1').descr('beam_power', '');
model.variable('var1').set('fmod', num2str(f_mod));
model.variable('var1').descr('fmod', '');
> The next section deals with the geometry. In this case the model is a 1D axisymmetric model. Hence one has the the digit 1 which specifies the dimension and axisymmetric is set to 'true'
model.geom.create('geom1', 1);
model.geom('geom1').axisymmetric(true);
> The run() function is to build the geometry
model.geom('geom1').run;
> The 'Interval' feature is present in the case of 1D models. The following line creates an interval feature called 'i1'
model.geom('geom1').feature.create('i1', 'Interval');
model.geom('geom1').feature('i1').set('p2', num2str(param.mirror.radius));
The above line is used to set the value of the right endpoint as the radius of the mirror. The string 'p2' is used to denote the right end point..
So this step creates a cylinder of radius equal to the radius of the mirror. In case of an infinite mirror an extra interval is added from the the radius to twice the radius.
It would help if more description on "interval" is left as comments and on why the extra interval was added and how does it make the calculation different.
> The interval is run using the 'run()' command.
> The material definition is added in the following section. Probably the data seems to have been taken from an external file rather than the COMSOL material library. Details of the file would be of help from the author.
The following lines cannot be understood. Comsol help on"propertyGroup" gives no results found
model.material('mat1').propertyGroup('def').set('heatcapacity', [num2str(param.material.specific_heat_per_volume) '[J/(kg*K)]']);
model.material('mat1').propertyGroup('def').set('density', [num2str(param.material.density) '[kg/m^3]']);
model.material('mat1').propertyGroup('def').set('thermalconductivity', [num2str(param.material.thermal_conductivity) '[W/(m*K)]']);
> The meshing is done on auto.
> As far as the Physics is concerned i.e. what kind of results we expect out of this geometry, the 'HeatTransfer' is selected.
A 1D heat source is created called 'hs1'
model.physics('ht').feature.create('hs1', 'HeatSource', 1);
and the Heat Source is time varying and given by beam_intensity*sin(2*fmod*pi*t). This can be seen in lines 121122
model.physics('ht').feature('hs1').selection.set([1]);
model.physics('ht').feature('hs1').set('Q', 1, 'beam_intensity*sin(2*fmod*pi*t)');
> If the mirror is infinite, something called 'InfiniteElements' is applied to the outer interval. More details on "infiniteElements" and why was it used would be helpful.
> A study involving the feature "Transient" was used. However, it is required to know what study was implemented from the GUI since the term "Transient" was not found under Study Steps in the GUI
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 
30

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

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

Attachment 1: Custom_Closed.JPG


Attachment 2: custom_open.JPG


Attachment 3: Stock_Closed.JPG


Attachment 4: Stock_Open.JPG


16

Mon Jun 10 18:37:09 2013 
Matt A.  General  General  Response 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.

56

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


110

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

111

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

112

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