{"id":354,"date":"2023-12-06T07:06:18","date_gmt":"2023-12-06T06:06:18","guid":{"rendered":"https:\/\/logbooks.ifosim.org\/finesse\/?p=354"},"modified":"2023-12-06T07:06:18","modified_gmt":"2023-12-06T06:06:18","slug":"simulation-of-stress-induced-birefringence-in-fenicsx","status":"publish","type":"post","link":"https:\/\/logbooks.ifosim.org\/finesse\/2023\/12\/06\/simulation-of-stress-induced-birefringence-in-fenicsx\/","title":{"rendered":"Simulation of Stress-Induced Birefringence in FEniCSx"},"content":{"rendered":"\n<p>I used fenicsx to create stress simulations of a silicon test mass which could then be converted to birefringence. The model used was based on the silicon test masses created for use in the Gingin 7m 2um cavity project. These have a diameter of 10cm and a thickness of 3cm.  All the values used in the code are in SI units. Gravity acts in the negative z direction and the optical axis is the x direction.  I created the mesh using the python api for gmsh and then imported it into dolfinx<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>from mpi4py import MPI\nfrom dolfinx.io.gmshio import model_to_mesh\nimport pyvista\nfrom dolfinx import mesh, fem, plot, io, default_scalar_type\nfrom dolfinx.fem.petsc import LinearProblem\nfrom mpi4py import MPI\nimport ufl\nimport numpy as np\nimport gmsh\nimport warnings\nwarnings.filterwarnings(\"ignore\")\ngmsh.initialize()\n\ngmsh.model.add(\"test_mass\")\ncad = gmsh.model.geo\nlc = 0.002\ncad.addPoint(0, -0.05, 0.01, lc, 1)\ncad.addPoint(0, 0.05, 0.01, lc, 2)\ncad.addPoint(0, 0.05, -0.01, lc, 3)\ncad.addPoint(0, -0.05, -0.01, lc, 4)\ncad.addPoint(0,0,0, lc, 5)\ncad.addLine(2, 3, 1)\ncad.addLine(4, 1, 2)\ncad.addCircleArc(1, 5, 2, 3)\ncad.addCircleArc(3, 5, 4, 4)\ncad.addCurveLoop(&#091;3, 1, 4, 2 ], 10)\ncad.addPlaneSurface(&#091;10], 11)\nvol = cad.extrude(&#091;(2, 11)], 0.03,0,0 )\ncad.synchronize()\ngmsh.model.addPhysicalGroup(3, &#091;vol&#091;1]&#091;1]], 41, \"volume\")\ngmsh.model.addPhysicalGroup(2, &#091;vol&#091;0]&#091;1],vol&#091;2]&#091;1], vol&#091;3]&#091;1], vol&#091;4]&#091;1], vol&#091;5]&#091;1]], 42, \"surface\")\ngmsh.model.mesh.generate(3)\nmodel_rank = 0\ndomain, cell_tags, facet_tags = model_to_mesh(\n    gmsh.model, MPI.COMM_WORLD, model_rank,3)<\/code><\/pre>\n\n\n\n<p>I then used dolfinx to run the finite element analysis. This produces a measurement of the deformation of each element. I set the boundary conditions so that the flat edges of the test mass were fixed. The code for the linear elasticity calculation is based on this dolfinx tutorial:<a href=\"https:\/\/jsdokken.com\/dolfinx-tutorial\/chapter2\/linearelasticity_code.html\" target=\"_blank\" rel=\"noreferrer noopener\"> https:\/\/jsdokken.com\/dolfinx-tutorial\/chapter2\/linearelasticity_code.html<\/a><\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>E = 130e9\nnu = 0.27\nmu = E \/ (2.0 * (1.0 + nu))\nlambda_ = E * nu \/ ((1.0 + nu) * (1.0 - 2.0 * nu))\nrho = 2329\ng = 9.8\nV = fem.VectorFunctionSpace(domain, (\"Lagrange\", 1))\ndef clamped_boundary(x):\n    return np.logical_or(np.isclose(x&#091;1], -0.05), np.isclose(x&#091;1], 0.05))\n\n\nfdim = domain.topology.dim - 1\nboundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)\nu_D = np.array(&#091;0, 0, 0], dtype=default_scalar_type)\nbc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)\nT = fem.Constant(domain, default_scalar_type((0, 0, 0)))\nds = ufl.Measure(\"ds\", domain=domain)\ndef epsilon(u):\n    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)\n\n\ndef sigma(u):\n    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)\n\n\nu = ufl.TrialFunction(V)\nv = ufl.TestFunction(V)\nf = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))\na = ufl.inner(sigma(u), epsilon(v)) * ufl.dx\nL = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds\nproblem = LinearProblem(a, L, bcs=&#091;bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\nuh = problem.solve()<\/code><\/pre>\n\n\n\n<p>The deformation can then be plotted in pyvista which gives:<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>pyvista.start_xvfb()\n\n# Create plotter and pyvista grid\np = pyvista.Plotter()\ntopology, cell_types, geometry = plot.vtk_mesh(V)\ngrid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n\n# Attach vector values to grid and warp grid by vector\ngrid&#091;\"u\"] = uh.x.array.reshape((geometry.shape&#091;0], 3))\n#actor_0 = p.add_mesh(grid, style=\"wireframe\", color=\"k\")\nwarped = grid.warp_by_vector(\"u\", factor=1.5)\nactor_1 = p.add_mesh(warped, show_edges=False)\np.show_axes()\nif not pyvista.OFF_SCREEN:\n    p.show()\nelse:\n    figure_as_array = p.screenshot(\"deflection.png\")<\/code><\/pre>\n\n\n\n<figure class=\"wp-block-image size-medium is-resized\"><img loading=\"lazy\" decoding=\"async\" width=\"300\" height=\"217\" src=\"https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/deformation-1-300x217.png\" alt=\"\" class=\"wp-image-356\" style=\"width:840px;height:auto\" srcset=\"https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/deformation-1-300x217.png 300w, https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/deformation-1-768x556.png 768w, https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/deformation-1.png 829w\" sizes=\"auto, (max-width: 300px) 100vw, 300px\" \/><\/figure>\n\n\n\n<p>Now we have to convert from the deformation to the stress tensor by interpolating onto a new vector function space. This gives the flattened list of the stress tensor components for each element. This is reshaped into a 3&#215;3 tensor.<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>gdim = domain.geometry.dim\nStress = fem.VectorFunctionSpace(domain, (\"Discontinuous Lagrange\", 0, (gdim,)), 9)\nstress = fem.Function(Stress)\nstress_expr = fem.Expression(sigma(uh), Stress.element.interpolation_points())\nstress_array = stress.x.array.reshape((int(stress.x.array.shape&#091;0]\/9), 3, 3))<\/code><\/pre>\n\n\n\n<p>Next we define the functions to calculate the birefringence. The calculation for the stress induced birefringence first calculates the change in the inverse dielectric tensor (optical index ellipsoid) and then calculates the eiegenvalues of the rotated tensor. The eigenvalues can then be used for a calculation of the birefringence.<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>def rot_mat(rot_z, rot_y, rot_x):\n    a = rot_z\n    b = rot_y\n    y = rot_x\n    mat = np.array(&#091;&#091;np.cos(a)*np.cos(b), np.cos(a)*np.sin(b)*np.sin(y) - np.sin(a)*np.cos(y), np.cos(a)*np.sin(b)*np.cos(y) + np.sin(a)*np.sin(y)],&#091;np.sin(a)*np.cos(b), np.sin(a)*np.sin(b)*np.sin(y) + np.cos(a)*np.cos(y), np.sin(a)*np.sin(b)*np.cos(y) - np.cos(a)*np.sin(y)], &#091;-1*np.sin(b), np.cos(b)*np.sin(y), np.cos(b)*np.cos(y)]])\n    return mat\n\n\ndef h_eigenvalues(mat):\n    a = mat&#091;0, 0]\n    b = mat&#091;1, 1]\n    c = mat&#091;0, 1]\n\n    l1 = (a + b - np.sqrt(4 * c ** 2 + (a - b) ** 2)) \/ 2\n    l2 = (a + b + np.sqrt(4 * c ** 2 + (a - b) ** 2)) \/ 2\n    return l1, l2\n\n\ndef r_index(mat):\n    l1, l2 = h_eigenvalues(mat)\n    return (1\/np.sqrt(l1)) - (1\/np.sqrt(l2))\n\n\ndef stress(s1, s2, s3, s4, s5, s6):\n    arr = np.array(&#091;s1, s2, s3, s4, s5, s6])\n    return arr\n\n\ndef photoelastic(p11, p12, p44):\n    arr = np.array(&#091;&#091;p11, p12, p12, 0, 0, 0], &#091;p12, p11, p12, 0,0,0], &#091;p12, p12, p11, 0,0,0],&#091;0,0,0,p44,0,0],&#091;0,0,0,0,p44,0],&#091;0,0,0,0,0,p44]])\n    return arr\n\n\ndef elastic(PR, YM):\n    arr = np.array(&#091;&#091;1, -1*PR, -1*PR,0,0,0],&#091;-1*PR, 1, -1*PR, 0,0,0], &#091;-1*PR, -1*PR, 1,0,0,0],&#091;0,0,0,1+PR,0,0],&#091;0,0,0,0,1+PR,0],&#091;0,0,0,0,0,1+PR]])\n    arr = arr * (1\/YM)\n    return arr\n\n\ndef delta_b(PR, YM, p11, p12, p44, s1, s2, s3, s4, s5, s6):\n    arr = np.dot(photoelastic(p11, p12, p44), np.dot(elastic(PR, YM), stress(s1, s2, s3, s4, s5, s6)))\n    return arr\n\n\ndef b_0(nx, ny, nz):\n    arr = np.array(&#091;1\/(nx**2), 1\/(ny**2), 1\/(nz**2), 0,0,0])\n    return arr\n\n\ndef b_new(PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz):\n    delB = delta_b(PR, YM, p11, p12, p44, s1, s2, s3, s4, s5, s6)\n    B0 = b_0(nx, ny, nz)\n    B = delB + B0\n    arr2 = np.array(&#091;&#091;B&#091;0], B&#091;5], B&#091;4]], &#091;B&#091;5], B&#091;1], B&#091;3]], &#091;B&#091;4], B&#091;3], B&#091;2]]])\n    return arr2\n\n\ndef birefringence(rot_z, rot_y, rot_x, PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz):\n    b = b_new(PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz)\n    rot = rot_mat(rot_z, rot_y, rot_x)\n    rotated = np.dot(rot, np.dot(b, np.linalg.inv(rot)))\n    small_b = rotated&#091;1:, 1:]\n    delta_n = r_index(small_b)\n    return delta_n<\/code><\/pre>\n\n\n\n<p>Now we just have to find the birefringence for each element in the interpolated mesh. The first three elements of the <code>birefringence <\/code>function are the rotation of the optic around the z, x and y axes. The next two are the Poisson ratio and the Young&#8217;s modulus of silicon. The next six are the elements of the stress tensor in voigt notation. Then the next three are the non-zero elements of the photoelastic tensor for cubic crystals. The last three are the refractive indicies in the x, y and z direction.<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>npoints = stress_array.shape&#091;0]\ndn = np.zeros(npoints)\nfor x in range(npoints):\n    a = stress_array&#091;x]\n    dn&#091;x] = birefringence(rot_z=0, rot_x=0, rot_y=0, PR=0.27, YM=130e9, s1=a&#091;0,0], s2=a&#091;1,1], s3=a&#091;2,2],s4= a&#091;1,2], s5=a&#091;0,2],s6=a&#091;0,1], p11=-0.094, p12=0.017,p44= -0.051, nx=3.48, ny=3.48, nz=3.48)\n<\/code><\/pre>\n\n\n\n<p>Having calculated the birefringence we can now plot it onto the mesh. The simulation of birefringence gives a result comparable with past measurements (<a href=\"https:\/\/dcc.ligo.org\/LIGO-P2200357\" target=\"_blank\" rel=\"noreferrer noopener\">https:\/\/dcc.ligo.org\/LIGO-P2200357<\/a>).<\/p>\n\n\n\n<pre class=\"wp-block-code\"><code>pyvista.global_theme.cmap = 'jet'\nwarped.cell_data&#091;\"Birefringence\"] = dn\nwarped.set_active_scalars(\"Birefringence\")\np = pyvista.Plotter()\np.add_mesh(warped, clim=&#091;0, 5e-8])\np.show_axes()\nif not pyvista.OFF_SCREEN:\n    p.show()\nelse:\n    stress_figure = p.screenshot(f\"stresses.png\")<\/code><\/pre>\n\n\n\n<figure class=\"wp-block-image size-full\"><img loading=\"lazy\" decoding=\"async\" width=\"813\" height=\"600\" src=\"https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/birefringence_sim.png\" alt=\"\" class=\"wp-image-357\" srcset=\"https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/birefringence_sim.png 813w, https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/birefringence_sim-300x221.png 300w, https:\/\/logbooks.ifosim.org\/finesse\/wp-content\/uploads\/sites\/5\/2023\/12\/birefringence_sim-768x567.png 768w\" sizes=\"auto, (max-width: 813px) 100vw, 813px\" \/><\/figure>\n","protected":false},"excerpt":{"rendered":"<p>I used fenicsx to create stress simulations of a silicon test mass which could then be converted to birefringence. The model used was based on the silicon test masses created for use in the Gingin 7m 2um cavity project. These have a diameter of 10cm and a thickness of 3cm. All the values used in [&hellip;]<\/p>\n","protected":false},"author":100,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"ssl_alp_hide_revisions":false,"footnotes":"","ssl_alp_hide_crossreferences_to":false},"categories":[1],"tags":[87,90,86,89,88],"ssl-alp-coauthor":[80],"class_list":["post-354","post","type-post","status-publish","format-standard","hentry","category-uncategorised","tag-birefringence","tag-dolfinx","tag-fenics","tag-fenicsx","tag-stress"],"_links":{"self":[{"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/posts\/354","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/users\/100"}],"replies":[{"embeddable":true,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/comments?post=354"}],"version-history":[{"count":6,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/posts\/354\/revisions"}],"predecessor-version":[{"id":363,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/posts\/354\/revisions\/363"}],"wp:attachment":[{"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/media?parent=354"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/categories?post=354"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/tags?post=354"},{"taxonomy":"ssl-alp-coauthor","embeddable":true,"href":"https:\/\/logbooks.ifosim.org\/finesse\/wp-json\/wp\/v2\/ssl-alp-coauthor?post=354"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}