Simulation of Stress-Induced Birefringence in FEniCSx

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

from mpi4py import MPI
from import model_to_mesh
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import gmsh
import warnings

cad = gmsh.model.geo
lc = 0.002
cad.addPoint(0, -0.05, 0.01, lc, 1)
cad.addPoint(0, 0.05, 0.01, lc, 2)
cad.addPoint(0, 0.05, -0.01, lc, 3)
cad.addPoint(0, -0.05, -0.01, lc, 4)
cad.addPoint(0,0,0, lc, 5)
cad.addLine(2, 3, 1)
cad.addLine(4, 1, 2)
cad.addCircleArc(1, 5, 2, 3)
cad.addCircleArc(3, 5, 4, 4)
cad.addCurveLoop([3, 1, 4, 2 ], 10)
cad.addPlaneSurface([10], 11)
vol = cad.extrude([(2, 11)], 0.03,0,0 )
gmsh.model.addPhysicalGroup(3, [vol[1][1]], 41, "volume")
gmsh.model.addPhysicalGroup(2, [vol[0][1],vol[2][1], vol[3][1], vol[4][1], vol[5][1]], 42, "surface")
model_rank = 0
domain, cell_tags, facet_tags = model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, model_rank,3)

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:

E = 130e9
nu = 0.27
mu = E / (2.0 * (1.0 + nu))
lambda_ = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
rho = 2329
g = 9.8
V = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
def clamped_boundary(x):
    return np.logical_or(np.isclose(x[1], -0.05), np.isclose(x[1], 0.05))

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
    return ufl.sym(ufl.grad(u))  # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L =, v) * ufl.dx +, v) * ds
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

The deformation can then be plotted in pyvista which gives:


# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
#actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
actor_1 = p.add_mesh(warped, show_edges=False)
if not pyvista.OFF_SCREEN:
    figure_as_array = p.screenshot("deflection.png")

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×3 tensor.

gdim = domain.geometry.dim
Stress = fem.VectorFunctionSpace(domain, ("Discontinuous Lagrange", 0, (gdim,)), 9)
stress = fem.Function(Stress)
stress_expr = fem.Expression(sigma(uh), Stress.element.interpolation_points())
stress_array = stress.x.array.reshape((int(stress.x.array.shape[0]/9), 3, 3))

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.

def rot_mat(rot_z, rot_y, rot_x):
    a = rot_z
    b = rot_y
    y = rot_x
    mat = np.array([[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)],[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)], [-1*np.sin(b), np.cos(b)*np.sin(y), np.cos(b)*np.cos(y)]])
    return mat

def h_eigenvalues(mat):
    a = mat[0, 0]
    b = mat[1, 1]
    c = mat[0, 1]

    l1 = (a + b - np.sqrt(4 * c ** 2 + (a - b) ** 2)) / 2
    l2 = (a + b + np.sqrt(4 * c ** 2 + (a - b) ** 2)) / 2
    return l1, l2

def r_index(mat):
    l1, l2 = h_eigenvalues(mat)
    return (1/np.sqrt(l1)) - (1/np.sqrt(l2))

def stress(s1, s2, s3, s4, s5, s6):
    arr = np.array([s1, s2, s3, s4, s5, s6])
    return arr

def photoelastic(p11, p12, p44):
    arr = np.array([[p11, p12, p12, 0, 0, 0], [p12, p11, p12, 0,0,0], [p12, p12, p11, 0,0,0],[0,0,0,p44,0,0],[0,0,0,0,p44,0],[0,0,0,0,0,p44]])
    return arr

def elastic(PR, YM):
    arr = np.array([[1, -1*PR, -1*PR,0,0,0],[-1*PR, 1, -1*PR, 0,0,0], [-1*PR, -1*PR, 1,0,0,0],[0,0,0,1+PR,0,0],[0,0,0,0,1+PR,0],[0,0,0,0,0,1+PR]])
    arr = arr * (1/YM)
    return arr

def delta_b(PR, YM, p11, p12, p44, s1, s2, s3, s4, s5, s6):
    arr =, p12, p44),, YM), stress(s1, s2, s3, s4, s5, s6)))
    return arr

def b_0(nx, ny, nz):
    arr = np.array([1/(nx**2), 1/(ny**2), 1/(nz**2), 0,0,0])
    return arr

def b_new(PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz):
    delB = delta_b(PR, YM, p11, p12, p44, s1, s2, s3, s4, s5, s6)
    B0 = b_0(nx, ny, nz)
    B = delB + B0
    arr2 = np.array([[B[0], B[5], B[4]], [B[5], B[1], B[3]], [B[4], B[3], B[2]]])
    return arr2

def birefringence(rot_z, rot_y, rot_x, PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz):
    b = b_new(PR, YM, s1, s2, s3, s4, s5, s6, p11, p12, p44, nx, ny, nz)
    rot = rot_mat(rot_z, rot_y, rot_x)
    rotated =,, np.linalg.inv(rot)))
    small_b = rotated[1:, 1:]
    delta_n = r_index(small_b)
    return delta_n

Now we just have to find the birefringence for each element in the interpolated mesh. The first three elements of the birefringence function are the rotation of the optic around the z, x and y axes. The next two are the Poisson ratio and the Young’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.

npoints = stress_array.shape[0]
dn = np.zeros(npoints)
for x in range(npoints):
    a = stress_array[x]
    dn[x] = birefringence(rot_z=0, rot_x=0, rot_y=0, PR=0.27, YM=130e9, s1=a[0,0], s2=a[1,1], s3=a[2,2],s4= a[1,2], s5=a[0,2],s6=a[0,1], p11=-0.094, p12=0.017,p44= -0.051, nx=3.48, ny=3.48, nz=3.48)

Having calculated the birefringence we can now plot it onto the mesh. The simulation of birefringence gives a result comparable with past measurements (

pyvista.global_theme.cmap = 'jet'
warped.cell_data["Birefringence"] = dn
p = pyvista.Plotter()
p.add_mesh(warped, clim=[0, 5e-8])
if not pyvista.OFF_SCREEN:
    stress_figure = p.screenshot(f"stresses.png")

Leave a Reply