ID 
Date 
Author 
Type 
Category 
Subject 
72

Thu Jul 25 15:54:58 2013 
Deep Chatterjee  Optics  General  TR results for different dimensions 
Quote: 
Quote: 
In this post I simulate the procedure of calculating the TR noise for finite cavities as proposed by Heinert and check for a
match.
The technique of performing all necessary calculations in COMSOL and exporting the results was applied to the TR codes.
It was seen that the codes gives similar output as the technique of extraction of Fourier coefficients in place of time averaging
as has been done in the codes of Koji Arai. One can see the output as the present code runs to be similar to the previous ones
found in the SVN.
However, the results in the present case were off by a constant factor close to 100. This maybe due to some 'm'  'cm' or similar difference between
analytic calculations and COMSOL values of parameters. Although, it has not been found yet, the correction is hopeful to be
found soon.
The codes give results similar to the analytic result for other values of the mirror radius and beam radii (apart from the constant
factor I have mentioned above). One may have a look at the trend of the graphs between analytic and simulated values in the plots
attached. These plots are for the case when the mirror radius = 25m while the beam radius = 9 cm i.e. the original radii were 0.25m
and 9cm respectively i.e. the ratio has been changed by a order of 2.
As mentioned before the reason for the constant factor difference will be looked into.

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

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


Attachment 2: Jul_25_2.pdf


Attachment 3: relative_plot.pdf


73

Mon Jul 29 22:42:57 2013 
Deep Chatterjee  Optics  General  Avoiding transient solutions in the Computation of TE/TR noise 
An error was being encountered in the computation of the TR noise lately. It was observed that while running the simulations in the case of the materials which have a lower value of the thermal diffusivity (silicon / sapphire at room temperature), the simulated result were slightly off from the analytic result. On the other hand, if the simulation was run with a material of higher diffusivity(same materials at lower temperature), the match would be better. The reason being the transient solution not dying off significantly during the period of the simulation. Since a time average was being taken of the quantity integral{grad_T ^2}, the transient contributed to the integral. To get the correct value, the fourier coefficient of the time signal of integral{grad_T ^2}, was extracted at twice the frequency of the pressure oscillation. The reason being that the signal was squared. Extracting the response at this frequency after the integration is logical since the integration is over space while the response we extract is over time.
The same procedure was also applied to the TE noise calculation. However, this time we obtained similar result as the case where this procedure was not applied but a simple time averaging was performed. The tail of the plot, at high frequency, is still seen to deviate from the analytic result of Liu and Thorne as was the case previously. A plot is attached showing the spectrum for Fused silica at 290K. The conclusion being that the transients do not affect the TE noise calculation  the plot stayed the same even after filtering them out. This is probably because unlike the TR case which has a heat source present along the cylinder axis, the TE noise calculation involves applying pressure only on the face of the cylinder, and the transient do not contribute much to the volume integral.

Attachment 1: Jul29_2.pdf


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


75

Thu Aug 8 17:17:19 2013 
Deep Chatterjee  Optics  General  Something like cancellation 
For the material parameters of Sapphire at 300K, the TE and TR Noise profiles, though not very close, lie close to within an order of magnitude. Sapphire has a positive coefficient of linear expansion. We just inverted the sign of this quantity
and the ran the codes that puts heat and pressure simultaneously to the test mass. The total noise looks to be lower than the TR noise which is greater.
Mirror radius 0.25[m]
Mirror height 0.46[m]
Beam Radius 0.09[m]
If we have physical parameters which make the two Noise sources come closer to each other and then flip the sign of alpha, we may be able to see some noise reduction to a greater extent. 
Attachment 1: suspected_cance_AUg8l.pdf


100

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

2

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

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

5

Mon Jun 3 21:01:46 2013 
Emory Brown  Optics  General  Analytic Calculation of Thermal Noise due to Brownian Motion using Levin and Thorne's Methods 
Motivation:
Reduction of Brownian thermal noise in future gravitational wave detectors is of significant interest. It has been suggested that changing the shape of the mirrors used may reduce the Brownian thermal noise. Before we can study how alterations in mirror shape effect Brownian thermal noise, we must be able to calculate the thermal noise analytically in order to compare it to our finite element model.
Methods:
I made Mathematica notebooks to perform calculation of the thermal noise. The first notebook implemented Levin's method directly.
Sx[f_] := (4*Kb*T/f)*(1  \[Sigma]^2)*1.87322*\[Phi]/(\[Pi]^3*E0*r0)
To test this against Levin's paper, the same values were used as in the paper, such that
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below. Additionally, the value for 100 Hz was explicitly computed and agreed with Levin's value of 8.7*10^40 m^2/Hz computed from his equation 15.
The notebook was modified to perform Thorne's method of computing the thermal noise for a finite test mass. This calculation was performed using equations 28, 35a, 54, 56, and 57. Equations 56 and 57 require the approximation made in eqn 37 which assumes a relatively low mass mirror. According to Liu and Thorne, this is a "rather good approximation for realistic parameter values." Performing the calculations again using this method gave Sx[100 Hz]=7.80081*10^40 m^2/Hz. So, for these parameters, the finite test mass provides a correction factor of about 10%. Another loglog plot of Sqrt[Sx[f]] against f was made using this method.
We can now modify the parameters above to match the values in our finite element model to verify the results of our finite element model.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.

Attachment 1: ComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.png


Attachment 2: DifferencePlot.png


Attachment 3: ThorneComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.nb

(* Contenttype: application/vnd.wolfram.mathematica *)
(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)
(* CreatedBy='Mathematica 9.0' *)
(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
... 2754 more lines ...

6

Tue Jun 4 17:25:11 2013 
Emory Brown  Optics  General  Comparison of First Results from COMSOL Model to the Analytic Solution 
Motivation: Before we can attempt to modify the mirror designs to reduce the thermal noise due to Brownian motion in them, we must verify that our model works and that its results match with the analytically calculated ones.
Methods: I had previously constructed a model of a cylindrical mirror in COMSOL: http://nodus.ligo.caltech.edu:8080/COMSOL/3. I updated the values used in the COMSOL model to match those from Levin's paper so that they would match the values I had already computed analytically. These values are listed below:
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
Using the simplest of the boundary conditions I had attempted to implement, fixing the mirror face opposite the applied force, I ran a stationary solver on the model. After the solver had completed, I ran a volume integration of the strain energy in the mirror and obtained U_{max}=1.52887*10^10 J. I also ran a surface integral of the force applied to the surface to confirm that the total force, F_{0}, was the 1N that the COMSOL model had applied to the mirror's face, which it was.
Levin's equations 10 and 12 were combined to give S_{x}(f)=[(K_{b}*T) / (Pi*f)] * [U_{max} / F_{0}^2] * Phi
With the applied force of 1N and the value of U_{max}=1.52887*10^10 J, S_{x}(100 Hz)=2.0157*10^40 m^2 / Hz which agrees to within a factor of 4 with the results of the calculation based on Liu and Thorne's paper which gave a value of 7.80081*10^40 m^2 / Hz. A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below.
Future work:
The obvious next step to take in the project is to attempt to get the better boundary conditions to work, in particular the "gravitational" body load suggested by Liu and Thorne. We noted while working today that a much simpler case where we applied a force of 2N to one surface of a cylinder and an opposing force of 2N in the opposite direction did not converge. It may be worth working with this case and attempting to get it to converge in order to inform how we can make the more complicated case converge. If we can get that case to converge and it agrees with the analytic results, then we will be ready to start varying the relative sizes of the two mirror faces and determining the effect on thermal noise due to Brownian motion.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.

Attachment 1: ComsolComparisonPowerSpectralDensitySurfaceFluctuationsVsFrequency.png


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


22

Tue Jun 18 15:26:50 2013 
Emory Brown  Optics  General  Attempting to Implement Liu and Thorne's Boundary Conditions 
Liu and Thorne demonstrated in their paper that the optimal set of boundary conditions for a problem like the one I am working on is the balance the force applied to one of the mirror faces by an opposing force of equal magnitude distributed throughout the body of the mirror. When we first attempted to implement this boundary condition, we obtained an error from COMSOL indicating that a solution could not be generated as "The relative error is greater than the relative tolerance." I attempted to find a solution to this problem so that we can run our simulations using the correct boundary conditions.
I tried a number of settings changes and tweaks suggested on the COMSOL forums and in official documentation on the error (https://community.cmc.ca/docs/DOC1453), but without success. Most of the suggestions were not relevant to our model, and confusingly when the number of elements in the mesh was increased from 13628 to 24570 the relative residual increased contrary to our initial attempts in http://nodus.ligo.caltech.edu:8080/COMSOL/3. It seems like the next action to take is to attempt to generate a simpler model and attempt reproduce and fix the same error.
I also decided to run both an eigenfrequency and frequency domain solver on our model using improved meshing with 24570 elements. I first ran these solvers using the rigid back constraint and both of them returned Umax=1.49931*10^10 which gives S_{x}(100 Hz)=1.97673*10^40 m^2/Hz which agrees quite well with the solution obtained from our previous simulations which gave S_{x}(100 Hz)=2.0157*10^40 m^2 / Hz. This still significantly differs from the analytically obtained value of 7.80081*10^40 m^2 / Hz using Liu and Thorne's calculation.
I then changed the boundary condition to use Liu and Thorne's gravitational body load boundary condition and both of these methods returned Umax=1.17164*10^9 J which corresponds to S_{x}(100 Hz)=1.54472*10^39 m^2/HZ. This only differs from the analytic result by a factor of 1.98, which is a significant improvement over the previous simulations, but still a much more significant difference than I would hope for given that they now use the same boundary conditions and it seems unlikely that using a stationary state solver will give a different result when both the eigenfrequency and frequency domain solvers agree to 6 digits.
I will spend some time attempting to figure out why the stationary solution is returning the relative residual error. Afterwords, I will start to learn how to use Matlab in conjunction with COMSOL in order to allow us to vary the mirror's size as we want in a systematic manner. Also, Yanbei Chen suggested that we may want to add in a mirror coating later in the project and see how its Brownian noise responds to changes in the mirror substrate shape.
edit: Taking the advice of the error documentation, I also ran a time dependent solver which returned the same result as the eigenfrequency and frequency domain solvers.

23

Wed Jun 19 14:59:13 2013 
Emory Brown  General  General  The 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 2D 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 timedependent solver and simply view times after the system will have settled to a stationary state (#2 https://community.cmc.ca/docs/DOC1453). 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 suboptimal 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. 
26

Thu Jun 20 12:51:44 2013 
Emory Brown  General  General  The 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 nonstationary 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 timedependent solver for the rigid back boundary condition and computed the strain energy for the final two timesteps. 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: TimeDependent 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 2D 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 timedependent solver and simply view times after the system will have settled to a stationary state (#2 https://community.cmc.ca/docs/DOC1453). 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 suboptimal 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 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. 
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. 
33

Tue Jun 25 15:18:41 2013 
Emory Brown  Optics  General  A new boundary condition to attempt to avoid the relative residual error 
I contacted COMSOL support about our difficulties with the relative residual error and was told "Typically for solid mechanics this error is because you have not constrained the body enough. So any solution + a rigid body transformation is also a solution. To remove this nonuniqueness, you need a solutions modulo the rigid body transformations. So try and constrain the body somehow." This is what the gravitationally balanced body load was supposed to do, but using the stationary solver has been nonconvergent and the eigenfrequency solver has generated odd modes, with the first real modes being sheer modes. These outputs indicate a problem with our model. Matt noticed yesterday that the meshing in our model was slightly nonsymmetrical, but we initially dismissed it as not being significant since it was a minor difference which we did not think could account for the large errors we are facing. At further consideration though, if some asymmetry arising from any source, even a numerical or rounding error, were present, that could cause a slight rotation of the object which would cause the object to experience a torque and minor sheer force. The stationary solution would not converge in this case, but the eigenfrequecy solution might, and if it did it would make sense to see sheer modes for some of its low frequency eigenmodes, as we have. One possible solution to this problem is to change the way we mesh the material to ensure a symmetrical distribution of nodes in the xy plane, probably by extruding lower dimensional meshed systems into our model. I am unsure if we would be able to implement this solution once we start to change the size of the radii of the object's faces. An alternate solution is to find another set of boundary conditions which should be equivalent to the gravitational body load constraint, but which are stable relative to minor perturbations of the system's conditions. I think that I have found another set of boundary conditions which should work and not be too difficult to implement in COMSOL.
The sides of the object should not move in the direction orthogonal both to their displacement from the center of the circular plane of the object even with them in the z direction or in the z direction, so the edge displaced from the center in the y direction should be fixed in the x direction.
The object's center of mass should not move which because of the previous condition can be reduced to not moving in the z direction.
I think that these boundary conditions should either be compatible with, or replace Liu and Thorne's boundary conditions in our model. I am going to spend some time attempting to implement these boundary conditions, see if they converge, and see if adding Liu and Thorne's gravitationally balanced force with them makes any difference, either in results (which it shouldn't) or the amount of time required to run. I am not sure how long this will take to implement, but I don't expect it to take more than a day or two. 
36

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

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

Tue Jul 2 13:32:30 2013 
Emory Brown  Optics  General  Using Matlab to run our COMSOL simulations 
We want to be able to test a number of values for the ratio of the radii of the frustum which will be much easier once we can get Matlab to automate running our COMSOL model. I have been able to output the COMSOL model to a .m file and was able to figure out a way to export the Umax value we want to Matlab (see the attached .m file).
The model ran in Matlab after debugging it for a short while and returned the same result as the COMSOL model had by itself. I setup a script to loop over some values of RatioRadii and feed them to the model. The first value it fed to the model worked and returned 1.2099e10 J for RatioRadii=0.1, but for RatioRadii=0.2 the model had a
meshing error. I am going to have to see if there is some way we can modify the mesh to fix this error. When I briefly looked into it today, I found that for several RatioRadii which I tried, the object could be meshed for certain finenesses, but not others. I am going to spend tomorrow working on the first progress report, but after that is done I will
work on fixing this meshing error. 
Attachment 1: model_7_1_13.m

function [Umax] = model(ratioRadii)
% ratioRadii must be a scalar
% model_7_1_13.m
%
% Model exported on Jul 1 2013, 17:08 by COMSOL 4.3.2.189.
import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
... 220 more lines ...

Attachment 2: looper.m

function [RatioRValues,UmaxValues] = looper
% Iterates over a range of RatioRValues
RatioRValues=[0.1:0.1:10]
UmaxValues=[]
for r=RatioRValues,
UmaxValues = [UmaxValues model_7_1_13(r)];
UmaxValues
end
... 4 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. 
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


52

Wed Jul 10 14:15:20 2013 
Emory Brown  Optics  General  Comparison of Umax and the frequency of the lowest eigenmode against ratioR 
I was able to run the looper program designed last week over the interval ratioR=[0.5, 0.9] with a step size of 0.01 and obtained the plot below with a minimal value of Umax=1.4723e10 J obtained for ratioR=0.71, which is about 3.2% less than the value of Umax=1.5202e10 J when ratioR=1. This is not very significant, but the shifting upward of the frequencies of the lowest eigenmodes as described in the previous post is. Additionally, these curves both have their optimal points in the same region, with a ratioR value of about 0.7.
I also ran a higher resolution test to see how the frequency of the lowest real eigenmode responded to changes in ratioR over this region. The plot generated is also attached with a maximal (most favorable) value of 7210.0 Hz at ratioR=0.74.
Data from both of these runs is included in the attached tarball stored as txt files delineated by tabs.
The next few things to try are substituting silicon for fused silica as the material in the model and seeing how the curves change and seeing the response from reversing the orientation of the optimal test mass (ratioR~0.72). I expect that the orientation change will result in much higher noise and a lower minimal real eigenfrequency. 
Attachment 1: UmaxValuesVsRatioR.png


Attachment 2: FrequencyLowestEigenmodes.png


Attachment 3: data.tar.gz

55

Fri Jul 12 00:09:17 2013 
Emory Brown  Optics  General  Analysis of substrate Brownian noise in a silicon test mass with changing ratioR value 
I performed the same analysis we previously did on a test mass made of silicon (using the builtin silicon within COMSOL) and obtained the attached plots. The plots show more of a difference between the two solutions optimal points than for fused silica, but they are still quite close. All data is included in the attached tarball. The values for ratioR=1 or a cylindrical test mass place the frequency of the lowest real eigenmode at 8491 Hz and the Umax value of 6.5484*10^11 J. 
Attachment 1: SiliconEigenfrequencyPlot.png


Attachment 2: UmaxPlotSilicon.png


Attachment 3: Silicon.tar.gz

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


59

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


76

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

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


Attachment 2: Simple_Cantilever_200_mm_Long.PNG


Attachment 3: Comparison_between_ANSYS_and_Calculations_Fused_Silica_200_mm_L.PNG


107

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

108

Fri Jul 29 14:33:41 2016 
Joy Westland  Mechanics  General  A 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:
 Download the SolidWorks Zip Folders from the DCC
 https://dcc.ligo.org/login/index.shtml?entityID=https%3A%2F%2Fdcc.ligo.org%2Fshibbolethsp&return=https%3A%2F%2Fdcc.ligo.org%2FShibboleth.sso%2FLogin%3FSAMLDS%3D1%26target%3Dss%253Amem%253A15b6c314d87e3fa8b3768d89cb6b9836fe39c754 (LIGOD1000760v4)
 Open/export the zip folders into a different file
 For only the test mass and ears, open the file called: D0902456 ITM OPTIC WITH EARS ASSEMBLY. This will represent the bottom mass.
 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.
 To do this, resave the assembly but before pressing save, change the file to PRT extension instead of Assembly.
 Open the PRT file of the test mass, make a new sketch on the face of the test mass.
 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.
 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.
 In “Split Line” press the “sketch” you want to project onto the “face” of the object.
Importing SolidWorks into ANSYS:
 Once the part has been made, choose the “static structural analysis”.
 In the geometry: Import the test mass that was converted into a part
 Once imported, press the geometry again and “edit”
 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.
 Exit geometry module
Applying Gaussian Force:
 Enter the “model” module to edit the setup of the analysis
 Assign the material that you want for the test mass under “Geometry” à “Solid” under the “Material” and “Assignment”
 Under the “Coordinate Systems”, make a new coordinate system
 “Coordinate System:
 Click on the face, in this case the circle made on the test mass for the “geometry” and “apply”
 Under “Definition” à “Type” à Change from Cartesian to Cylindrical. This allows the Gaussian Force to be distributed correctly
 Under “Analysis Settings” Turn on “Large Deflection”
 Add “gravity” under “Static Structural”
 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.
 Insert “Fix” to the top surfaces of the ears
 Insert “Pressure”
 Press on the circle made to apply for the pressure
 Under “Definition” change the “Magnitude” to “Function” (the arrow at the end of the Magnitude entry bar)
 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.
 Under “Graph Controls” Make sure that the XAxis is changed from Time to X.
 Pick a range for the graph and the number of segments that you want to look at.
 “Solve” the system and right click on “Solution” and evaluate the different results you need.

Attachment 1: Bottom_Test_Mass_With_Ears.PNG


Attachment 2: Applying_Parameters.PNG


Attachment 3: Applying_Gaussian_Force.PNG


4

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


29

Sun Jun 23 14:29:04 2013 
Koji  Optics  General  Differences on using a coarser mesh in calculating TR noise 
What about the dependence on the size of the time slice?
It is unclear what the unit fo the frequency. Hz? ot logHz (which should be expressed as 10^1, 10^2, ...)?
Can you combine the result of the caluculation and the error together?
i.e. Combine the plot similar to http://nodus.ligo.caltech.edu:8080/COMSOL/4 and your plot with the horizontal axies lined up. 
45

Wed Jul 3 00:04:46 2013 
Koji  Optics  General  Using Matlab to run our COMSOL simulations 
Please use zipped files instead of uploading extremely long text files on the elog.
Or use wiki or svn instead of elog for uploading them.
And appropriate word wrapping so that we can read the entries. Please. 
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.) 
90

Sun Jun 29 20:25:44 2014 
Koji  Optics  General  Difficulty with the COMSOL stationary module; Test Cases 
What about this example? The result is easier to understand intuitively.
Consider a bar with the length of L.
Let's say there is no body heat applied, but the temperature of the bar at x=L is kept at T=0
and at x=0 is kept at T=T0 Exp[I w t].
The equation for the bar is
...(1)
Consider the solution with the form of T(x, t) = T(x) T0 Exp(I w t), where T(x) is the position dependent transfer function.
T(x) is a complex function.
Eq.1 is modified with T(x) as
With the boundary condition of
This can be analytically solved in the following form
where alpha is defined by
So kappa/Cp is the characteristic (angular) frequency of the system.
Here is the example plot for L=1 and alpha = 1 (red), 10 (yellow green), 100 (turquoise), 1000 (blue)
If the oscillation is slow enough, the temperature decay length is longer than the bar length and thus the temperature is linear to the position.
If the oscillation is fast, the decay length is significantly shorter than the bar length and the temperature dependence on the position is exponential.
Now what we need is to solve this in COMSOL 
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.) 
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

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

113

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

7

Wed Jun 5 13:31:41 2013 
Matt A  Optics  General  Analytic Calculation of Thermal Noise due to Brownian Motion using Levin and Thorne's Methods 
Nice Work Emory,
Can you make a plot with both of the lines on the same plot, with a legend, units, and everything? I'd like to see these compared sideby side.
Quote: 
Motivation:
Reduction of Brownian thermal noise in future gravitational wave detectors is of significant interest. It has been suggested that changing the shape of the mirrors used may reduce the Brownian thermal noise. Before we can study how alterations in mirror shape effect Brownian thermal noise, we must be able to calculate the thermal noise analytically in order to compare it to our finite element model.
Methods:
I made Mathematica notebooks to perform calculation of the thermal noise. The first notebook implemented Levin's method directly.
Sx[f_] := (4*Kb*T/f)*(1  \[Sigma]^2)*1.87322*\[Phi]/(\[Pi]^3*E0*r0)
To test this against Levin's paper, the same values were used as in the paper, such that
Kb=1.38065*10^23 (Boltzmann's constant)
T=300 (Temperature)
Sigma=0.16 (Poison's Ratio)
E0=71.8*10^9 (Young's Modulus)
r0=0.0156 (Gaussian Beam Size)
Phi=10^7 (Loss Angle)
RMirror=0.17 (Mirror Radius)
H=0.2 (Mirror Height)
A loglog plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below. Additionally, the value for 100 Hz was explicitly computed and agreed with Levin's value of 8.7*10^40 m^2/Hz computed from his equation 15.
The notebook was modified to perform Thorne's method of computing the thermal noise for a finite test mass. This calculation was performed using equations 28, 35a, 54, 56, and 57. Equations 56 and 57 require the approximation made in eqn 37 which assumes a relatively low mass mirror. According to Liu and Thorne, this is a "rather good approximation for realistic parameter values." Performing the calculations again using this method gave Sx[100 Hz]=7.80081*10^40 m^2/Hz. So, for these parameters, the finite test mass provides a correction factor of about 10%. Another loglog plot of Sqrt[Sx[f]] against f was made using this method.
We can now modify the parameters above to match the values in our finite element model to verify the results of our finite element model.
Levin, Y. (1998). Internal thermal noise in the LIGO test masses: A direct approach. Physical Review D, 57(2), 659.
Liu, Y. T., & Thorne, K. S. (2000). Thermoelastic noise and homogeneous thermal noise in finite sized gravitationalwave test masses. Physical Review D, 62(12), 122002.


9

Wed Jun 5 16:56:53 2013 
Matt A.  Optics  General  Difference in the Levin and Liu & Thorne results for thermal noise 
Good Work Deep,
Can you include the equations that you used to calculated these expressions?
Quote: 
The analytical expressions for thermal noise, as calculated by Liu & Thorne and Levin, was plotted as a function of frequency in a log  log plot using Matlab.
The value of the parameters were used from Levin's paper on thermal noise.
phi = 1*10^{7} (loss angle)
r0 = 1.56*10^{2 }(beam radius)
T = 300 K (temperature)
E0 = 7.18*10^{10 } (Young's modulus)
sigma = 0.16 (Poisson ratio)
The value of S_{x}(100) that given in Levin's paper (8.7*10^{40} m^{2} /Hz) while the LiuThorne value is 9.1*10^{40} m^{2} /Hz.

The Liu Thorne expression being
S = 4*kb*T/(pi*f) * (1  sigma^2) / (2*sqrt(pi)*E*sqrt(2)*r0) * phi
That of Levin goes as
S = 4*kb*T/f * (1  sigma^2) / (pi^3*E*r0) *I*phi
sigma : Poisson ratio
phi : loss angle
r0 : Beam spot radius
E : Young's modulus
I : This is a sum the value of which is approx. 1.873.22 [Eqn.(A6) of Levin  Internal Thermal Noise] 
10

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

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

13

Thu Jun 6 12:40:07 2013 
Matt A.  Optics  General  Conventional Thermal noise (Sec V) from Liu & Thorne 
Good work Deep,
Can you write more on what this is and why you're doing it? We want our elog entries to be easy to understand for anybody who reads it.
Quote: 
The power spectrum of thermal noise in the case of Finite test masses differ slightly from that of the Infinite Test Mass. The expression for the PSD in the infinite test mass case is not a closed form solution but an infinite sum. In this post, a comparison has been made to check how fast does the sum (and as a result, the PSD) converge. The numerical values used for the calculation have been taken from the paper by Levin and that by Liu and Thorne.
>The linearized PSD plots are created for the case of thermal noise in finite test mass (Sec. V) of Liu and Thorne.
S = 8*kb*T*phi*(U0 + delta_U)/omega
>The maximum energy due to stress is considered by an infinite sum here. A comparison has been made regarding the convergence of the sum.
>The two plots correspond to the cases of considering 10 and 100 terms in the sum respectively.
>The plot shows that the difference is not much and hence convergence is fast.
>The relative difference is plotted w.r.t S_100, the PSD considering 100 terms in the sum.
>The relative difference goes like abs(S_100  S_10)/S_100, where 10 and 100 represent the number of terms considered in the sum.
>The algorithm used to evaluate the sum involving Bessel functions was the one by GWINC. (http://nodus.ligo.caltech.edu:8080/COMSOL/10).

I had made some minor errors in the expressions previously while typing the expressions. I have made the corrections and uploaded the new plots in place of the older ones. 
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.

24

Thu Jun 20 10:26:38 2013 
Matt A.  General  General  The 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 2D 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 timedependent solver and simply view times after the system will have settled to a stationary state (#2 https://community.cmc.ca/docs/DOC1453). 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 suboptimal 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.


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
