40m QIL Cryo_Lab CTN SUS_Lab TCS_Lab OMC_Lab CRIME_Lab FEA ENG_Labs OptContFac Mariner WBEEShop
  COMSOL elog, Page 1 of 3  Not logged in ELOG logo
ID Date Author Type Category Subjectup
  118   Mon Dec 4 16:27:13 2017 aaronMechanicsPonderSqueeze 

Meshing Surface Layers

Defining New Selections

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

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

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

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

Mesh Steps

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

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

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

Study

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

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

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

1. Open the ANSYS workbench

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  108   Fri Jul 29 14:33:41 2016 Joy WestlandMechanicsGeneralA Tutorial in Importing SolidWorks Files and Applying a Gaussian Force in ANSYS for a LIGO Test Mass

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

Using SolidWorks:

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

Importing SolidWorks into ANSYS:

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

Applying Gaussian Force:

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

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

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.
  107   Wed Jun 29 14:14:44 2016 Joy WestlandMechanicsAnalysisANSYS Tutorials with Basic Meshing

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

  1. ANSYS Meshing Fine Mesh Basic Tutorial 1: https://www.youtube.com/watch?v=sZIX3CJkWBE
  2. ANSYS Meshing Method Basic Tutorial 2: https://www.youtube.com/watch?v=ATSnvFUbEk4
  3. ANSYS Meshing Refinement Mesh Basic Tutorial 3: https://www.youtube.com/watch?v=QVfJiydCu2g
  4. ANSYS Meshing Sizing Mesh Basic Tutorial 4: https://www.youtube.com/watch?v=ODJ11YDC6ac
  5. ANSYS Meshing Spheres of Influence Basic Tutorial 5: https://www.youtube.com/watch?v=AQlGccK_X-0
  6. Body of Influnce 6: https://www.youtube.com/watch?v=3E1p1w32jt0
  7. Section Plane 7: https://www.youtube.com/watch?v=hhZOJSESJEU
  8. Mapped and Match Control 8: https://www.youtube.com/watch?v=3pogywR8Vgw
  9. MultiZone and Inflation 9: https://www.youtube.com/watch?v=vxVBvANHuB0
  10. Sizing Soft/Hard 10: https://www.youtube.com/watch?v=D8UzdC5Gmhs
  55   Fri Jul 12 00:09:17 2013 Emory BrownOpticsGeneralAnalysis 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 built-in silicon within COMSOL) and obtained the attached plots.  The plots show more of a difference between the two solutions optimal points than for fused silica, but they are still quite close.  All data is included in the attached tarball.  The values for ratioR=1 or a cylindrical test mass place the frequency of the lowest real eigenmode at 8491 Hz and the Umax value of  6.5484*10^-11 J.

  5   Mon Jun 3 21:01:46 2013 Emory BrownOpticsGeneralAnalytic 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 log-log 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 log-log 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 gravitational-wave test masses. Physical Review D, 62(12), 122002.
 
 
 
 
  7   Wed Jun 5 13:31:41 2013 Matt AOpticsGeneralAnalytic 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 side-by 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 log-log 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 log-log 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 gravitational-wave test masses. Physical Review D, 62(12), 122002.
 
 
 
 

 

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

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

and computes quantities as given in Heinert's paper.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

over which the function is defined.

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

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

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

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

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

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

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

* The coefficients are squared and added in the line 88

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

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

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

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

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

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

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

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

  106   Tue Jun 28 16:54:23 2016 Nikhil MathurMechanicsAnalysisAnsys 14.5 Indroductory Tutorial: Modal Frequency Convergence

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

  41   Fri Jun 28 14:03:07 2013 Deep ChatterjeeGeneralGeneralApproach for TE noise for beam in cavity

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

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

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

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

  46   Thu Jul 4 15:23:47 2013 Deep ChatterjeeOpticsGeneralAttempt to compute TE noise following Liu and Thorne

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

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

Steps followed by Liu and Thorne:

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

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

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

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

> PSD is calculated using FDT due to Levin

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

Liu&Thorne_type.bmp

Steps followed in the COMSOL model:

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

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

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

> Thermal expansion was enabled in the Solid Stress physics

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

> Heat insulation wsa applied to all surfaces.

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

> A custom mesh was put.

> The model was exported as '.m'

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

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

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

> All the prefactors were ignored while plotting

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

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

something_like_TE_1.bmp

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

something_like_TE_1_10times.bmp

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

something_like_TE_1_1000times.bmp

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

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

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

  47   Fri Jul 5 14:06:24 2013 Deep ChatterjeeOpticsGeneralAttempt to compute TE noise following Liu and Thorne

Quote:

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

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

Steps followed by Liu and Thorne:

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

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

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

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

> PSD is calculated using FDT due to Levin

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

Liu&Thorne_type.bmp

Steps followed in the COMSOL model:

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

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

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

> Thermal expansion was enabled in the Solid Stress physics

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

> Heat insulation wsa applied to all surfaces.

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

> A custom mesh was put.

> The model was exported as '.m'

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

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

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

> All the prefactors were ignored while plotting

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

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

something_like_TE_1.bmp

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

something_like_TE_1_10times.bmp

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

something_like_TE_1_1000times.bmp

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

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

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

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

comparison.bmp

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

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

L&T_COMSOL_difference.bmp

  22   Tue Jun 18 15:26:50 2013 Emory BrownOpticsGeneralAttempting 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/DOC-1453), but without success.  Most of the suggestions were not relevant to our model, and confusingly when the number of elements in the mesh was increased from 13628 to 24570 the relative residual increased contrary to our initial attempts in http://nodus.ligo.caltech.edu:8080/COMSOL/3.  It seems like the next action to take is to attempt to generate a simpler model and attempt reproduce and fix the same error.

I also decided to run both an eigenfrequency and frequency domain solver on our model using improved meshing with 24570 elements.  I first ran these solvers using the rigid back constraint and both of them returned Umax=1.49931*10^-10 which gives Sx(100 Hz)=1.97673*10^-40 m^2/Hz which agrees quite well with the solution obtained from our previous simulations which gave Sx(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 Sx(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.

 

  73   Mon Jul 29 22:42:57 2013 Deep ChatterjeeOpticsGeneralAvoiding 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.

 

  12   Wed Jun 5 20:54:59 2013 Deep ChatterjeeGeneralGeneralBessel Function roots

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

The difference between this and the algorithm by

Greg von Winckel goes as
Difference_between_algorithms.png

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

 

  42   Mon Jul 1 13:01:58 2013 Arnaldo RodriguezOpticsGeneralBrief Parameter Study on P-active Control Loop

 The following is a plot of the defocus as a function of time with only the proportional gain being nonzero. Only three cases were run before the program crashed.

K_p.png

The dotted line indicates the setpoint.

One can clearly see the loop becomes unstable and oscillates out of control at a proportional gain of 240, and can see the characteristic steady-state offset of the P-controller from the set point develop in the other cases.

  76   Mon Sep 16 10:07:32 2013 EvanGeneralConfigurationCOMSOL 4.3 on the OS X command line

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

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

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

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

comsol server > server_log.txt &

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

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

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

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

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

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

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

On the server:

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

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

On my laptop, from the MATLAB prompt:

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

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

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

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

  140   Tue Aug 11 16:35:07 2020 aaronGeneralConfigurationCOMSOL: remote server w/ matlab

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

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

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

export COMSOL_MATLAB_INIT=’matlab -nosplash’

 
4. You can run whatever comsol script you need. Make sure that in your script you import the comsol functions by calling the following
import com.comsol.model.*
import com.comsol.model.util.*
  1   Thu Dec 13 18:47:18 2012 ranaGeneralGeneralCOMSOL: who's using up all the licenses?

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

Specifically, with Mac COMSOL, I do:

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

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

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

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

etc.

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

Quote:

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

Specifically, with Mac COMSOL, I do:

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

 

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

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

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

etc.

 

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

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

CADIMPORT licenses currently in use: (2/2)

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

  61   Tue Jul 23 18:26:04 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


  67   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  66   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  65   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  64   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  63   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  62   Tue Jul 23 20:53:37 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  68   Tue Jul 23 20:53:45 2013 Deep ChatterjeeOpticsGeneralComparison between Liu and Thorne Results and COMSOL results for TE noise

Quote:

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

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

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

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

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

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

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

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

 

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

Jul_23.bmp


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

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

  6   Tue Jun 4 17:25:11 2013 Emory BrownOpticsGeneralComparison 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 Umax=1.52887*10^-10 J.  I also ran a surface integral of the force applied to the surface to confirm that the total force, F0,  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 Sx(f)=[(Kb*T) / (Pi*f)] * [Umax /  F0^2] * Phi

With the applied force of 1N and the value of Umax=1.52887*10^-10 J, Sx(100 Hz)=2.0157*10^-40 m^2 / Hz which agrees to within a factor of 4 with the results of the calculation based on Liu and Thorne's paper which gave a value of 7.80081*10^-40 m^2 / Hz.  A log-log plot of Sqrt[Sx[f]] with f ranging from 0.1 to 5000 Hz was plotted and is displayed below.

Future work:

The obvious next step to take in the project is to attempt to get the better boundary conditions to work, in particular the "gravitational" body load suggested by Liu and Thorne.  We noted while working today that a much simpler case where we applied a force of 2N to one surface of a cylinder and an opposing force of 2N in the opposite direction did not converge.  It may be worth working with this case and attempting to get it to converge in order to inform how we can make the more complicated case converge.  If we can get that case to converge and it agrees with the analytic results, then we will be ready to start varying the relative sizes of the two mirror faces and determining the effect on thermal noise due to Brownian motion.

 

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 gravitational-wave test masses. Physical Review D, 62(12), 122002.

 

 

 

 

 

 

 

 

 

 

 

  14   Thu Jun 6 17:00:11 2013 Deep ChatterjeeOpticsGeneralComparison of ITM and FTM geometries (FIG. 3 Liu & Thorne)

In section V of their paper on thermoelastic and homogenous thermal noise, Liu and Thorne have corrected the result of power spectral density of a Finite Test Mass from the result of BHV. Although the result of the infinte test mass has a closed form solution, that for the finite test mass is not closed and has to be approximated numerically. From the infinite series (Eq. (56), Eq.(57)) 100 terms were taken to approximate the sum. Considering that the convergence is fast (my last ELOG entry has a plot which shows little difference between considering 10  and 100 terms), 100 terms are a fair approximation.

The spectral density of thermal noise in the finite cavity is slightly lower than the corresponding infinite case. To highlight the difference between the two, they have plotted a quantity C_FTM which is a ratio of the linearized spectral densities of the Finite Test Mass and the Infinite Test Mass.

In this post, the same quantity has been numerically computed and plotted.

C_FTM.png

The above plot maybe compared to the one given on Fig. 3 of Liu & Thorne. Here r0 is the beam spot radius and C-FTM is the quantity mentioned previously. A snapshot of the figure from Liu & Thorne is shown below

Fig3_Liu&Thorn.bmp

The match between the figures seem fair.

Since C_FTM is a ratio, its frequency independent (both S_ITM and S_FTM have 1/f dependence on frequency which cancels on taking the ratio). This is verified in the above plot where plots are made for 10, 100 and 1000 Hz and they fairly coincide.

  52   Wed Jul 10 14:15:20 2013 Emory BrownOpticsGeneralComparison 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.4723e-10 J obtained for ratioR=0.71, which is about 3.2% less than the value of Umax=1.5202e-10 J when ratioR=1.  This is not very significant, but the shifting upward of the frequencies of the lowest eigenmodes as described in the previous post is.  Additionally, these curves both have their optimal points in the same region, with a ratioR value of about 0.7.

I also ran a higher resolution test to see how the frequency of the lowest real eigenmode responded to changes in ratioR over this region.  The plot generated is also attached with a maximal (most favorable) value of 7210.0 Hz at ratioR=0.74.

Data from both of these runs is included in the attached tarball stored as txt files delineated by tabs.

The next few things to try are substituting silicon for fused silica as the material in the model and seeing how the curves change and seeing the response from reversing the orientation of the optimal test mass (ratioR~0.72).  I expect that the orientation change will result in much higher noise and a lower minimal real eigenfrequency.

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

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

The method:

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

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

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

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

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

 

 

  11   Wed Jun 5 20:39:48 2013 Deep ChatterjeeOpticsGeneralConventional Thermal noise (Sec V) from Liu & 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.

Thorne_thermal_noise.png

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

Percentage_difference.png

>The algorithm used to evaluate the sum involving Bessel functions was the one by GWINC. (http://nodus.ligo.caltech.edu:8080/COMSOL/10).

  13   Thu Jun 6 12:40:07 2013 Matt A.OpticsGeneralConventional 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.

Thorne_thermal_noise.png

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

Percentage_difference.png

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

  39   Thu Jun 27 17:12:33 2013 Arnaldo RodriguezOpticsGeneralDefocusing in P-Controlled Manual Simulation

 To test the effectiveness of the PID loop, I've only given a nonzero value to the K_p constant and varied it over a series of runs to verify the response by the controller.

The following is a plot of the baseline compared to the P-active simulations for different proportionality gain constants. I've yet to determine whether all the trends in behavior are physical.

Each different case has an evident droop (steady-state error), and seem to be destabilizing as K_p increases. However, most cases reach steady-state quicker than the baseline.

Defocus_(P_Active).png

 It should be noted that the ring heater is not correcting any effects from the beam, as it is turned off. The ring heater is correcting itself; or rather, its own initial heat pulse of 5 W * delta_t (600s) = 3000 J before the controller becomes active.

  19   Wed Jun 12 13:52:44 2013 Deep ChatterjeeOpticsGeneralDifference in COMSOL and Analytic solution for TR noise

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

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

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

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

Also the difference is higher at the higher frequencies.

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

analytic~COMSOL.png

  8   Wed Jun 5 13:49:13 2013 Deep ChatterjeeOpticsGeneralDifference in the Levin and Liu & Thorne results for thermal noise

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*1010   (Young's modulus)

sigma = 0.16 (Poisson ratio)

The value of Sx(100) that given in Levin's paper (8.7*10-40 m2 /Hz) while the Liu-Thorne value is 9.1*10-40 m2 /Hz.

Levin_throne_comparison.png

  9   Wed Jun 5 16:56:53 2013 Matt A.OpticsGeneralDifference 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*1010   (Young's modulus)

sigma = 0.16 (Poisson ratio)

The value of Sx(100) that given in Levin's paper (8.7*10-40 m2 /Hz) while the Liu-Thorne value is 9.1*10-40 m2 /Hz.

Levin_throne_comparison.png

 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]

  27   Thu Jun 20 16:29:38 2013 Deep ChatterjeeOpticsGeneralDifferences on using a coarser mesh in calculating TR noise

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

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

difference_in_mesh.png

 

  29   Sun Jun 23 14:29:04 2013 KojiOpticsGeneralDifferences 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.

  88   Sat Jun 28 21:59:11 2014 Sam MooreOpticsGeneralDifficulty with the COMSOL stationary module; Test Cases

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

 

 

  89   Sun Jun 29 15:37:18 2014 Sam MooreOpticsGeneralDifficulty with the COMSOL stationary module; Test Cases

Quote:

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

 

 

 The following document shows the relative difference between these two plots.

  90   Sun Jun 29 20:25:44 2014 KojiOpticsGeneralDifficulty 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

  15   Fri Jun 7 17:41:41 2013 Deep ChatterjeeGeneralGeneralDisscussion of the code of a 'comsol with matlab' model file

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

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

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

that compute the analytical solution or the COMSOL simulations.

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

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

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

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

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

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

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

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

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

 

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

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

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

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

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

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

model = ModelUtil.create('Model');

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

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

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

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

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

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

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

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

model.param.set('beam_shape', '2/(pi*beam_size^2)*exp(-2*r^2/beam_size^2)');
model.param.descr('beam_shape', '');
model.param.set('beam_intensity', 'beam_power*beam_shape');
model.param.descr('beam_intensity', '');
model.param.set('dT', 'mod1.T-T_amb');
model.param.descr('dT', '');

An empty quote in the description implies that no description is given. The above statements define the parameters mentioned alongside them according to Heinerts paper.

>  [line 53-61]. the string 'var1' is used to tag all the global variables that are created. variables are expressions created out of the parameters.

Just like the above case, 'set()' is used to create variables and give them an expression, here as one can see 4 variables are created T_amb, beam_size, beam_power. and f_mod.

model.variable.create('var1');
model.variable('var1').set('T_amb', [num2str(param.material.temperature) '[K]']);
model.variable('var1').descr('T_amb', '');
model.variable('var1').set('beam_size', [num2str(param.beam.radius) '[m]']);
model.variable('var1').descr('beam_size', '');
model.variable('var1').set('beam_power', [num2str(param.COMSOL.beam_power) '[W]']);
model.variable('var1').descr('beam_power', '');
model.variable('var1').set('fmod', num2str(f_mod));
model.variable('var1').descr('fmod', '');

> The next section deals with the geometry. In this case the model is a 1D axisymmetric model. Hence one has the the digit 1 which specifies the dimension and axisymmetric is set to 'true'

model.geom.create('geom1', 1);
model.geom('geom1').axisymmetric(true);

> The run() function is to build the geometry

model.geom('geom1').run;

> The 'Interval' feature is present in the case of 1D models. The following line creates an interval feature called 'i1'

model.geom('geom1').feature.create('i1', 'Interval');

model.geom('geom1').feature('i1').set('p2', num2str(param.mirror.radius));

The above line is used to set the value of the right end-point as the radius of the mirror. The string 'p2' is used to denote the right end point..

So this step creates a cylinder of radius equal to the radius of the mirror. In case of an infinite mirror an extra interval is added from the the radius to twice the radius.

It would help if more description on "interval" is left as comments and on why the extra interval was added and how does it make the calculation different.

> The interval is run using the 'run()' command.

> The material definition is added in the following section. Probably the data seems to have been taken from an external file rather than the COMSOL material library. Details of the file would be of help from the author.

The following lines cannot be understood. Comsol help on"propertyGroup" gives no results found

model.material('mat1').propertyGroup('def').set('heatcapacity', [num2str(param.material.specific_heat_per_volume) '[J/(kg*K)]']);
model.material('mat1').propertyGroup('def').set('density', [num2str(param.material.density) '[kg/m^3]']);
model.material('mat1').propertyGroup('def').set('thermalconductivity', [num2str(param.material.thermal_conductivity) '[W/(m*K)]']);

> The meshing is done on auto.

> As far as the Physics is concerned i.e. what kind of results we expect out of this geometry, the 'HeatTransfer' is selected.

A 1D heat source is created called 'hs1'

model.physics('ht').feature.create('hs1', 'HeatSource', 1);

and the Heat Source is time varying and given by  beam_intensity*sin(2*fmod*pi*t). This can be seen in lines 121-122

model.physics('ht').feature('hs1').selection.set([1]);
model.physics('ht').feature('hs1').set('Q', 1, 'beam_intensity*sin(2*fmod*pi*t)');

> If the mirror is infinite, something called 'InfiniteElements' is applied to the outer interval. More details on "infiniteElements" and why was it used would be helpful.

> A study involving the feature "Transient" was used. However, it is required to know what study was implemented from the GUI since the term "Transient" was not found under Study Steps in the GUI

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

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

 

  93   Thu Jul 10 16:51:14 2014 Sam MooreOpticsGeneralDuan and Heinert Comparison

(See Plots in attached document)

 

My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V).  I am trying to match Duan's analytical expression with Heinert's analytical expression.  This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. )   Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".

It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.  
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation.  The calculation in Mathematica was very slow, so the upper limits on the integral were truncated.  The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
 
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steady-state method for Duan's numerical case (Heinert's plot runs much faster).
  96   Mon Jul 14 19:14:31 2014 Sam MooreOpticsGeneralDuan and Heinert Comparison

Quote:

(See Plots in attached document)

 

My plan has been to replicate Duan's numerical thermoconductive (TE + TR) phase noise plot presented in his paper (section V).  I am trying to match Duan's analytical expression with Heinert's analytical expression.  This requires some rescaling of Heinert's TR displacement noise. (I also needed to divide Heinert's expression by 4 pi to match the Fourier Transform convention. )   Duan's analytical expression for the phase noise is obtained by evaluating the triple integral given in equation 13 of the Duan paper "General Treatment of Thermal Noise in Optical Fibers".

It turns out that an additional factor of 2 multiplies the phase noise because Duan's Fourier Transform only takes into account positive frequencies; there are also negative frequencies that occur in equal amplitude.  
This integral was evaluated in Mathematica due to numerical noise in MATLAB's calculation.  The calculation in Mathematica was very slow, so the upper limits on the integral were truncated.  The following plots in the attached document show the resulting noise profile agreements for two different upper limits.
 
If the residual for the highest upper limit is considered acceptable for a match between the two plots, then I will use Heinert's plot as a reference when using the COMSOL steady-state method for Duan's numerical case (Heinert's plot runs much faster).

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

 

ELOG V3.1.3-