For the writeup, click here.
Thin film interference seen in soap-bubbles is typically a time-variant process. Environmental forces like wind or gravity generate velocity fields in the bubble's membrane causing its thickness to vary spatially and temporally. As a result, in nature we often see a splendor of chromatic patterns evolving over a fluid surface.
Inspired by such beauty, we are interested in marrying fluid dynamics with thin film interference to render a simulation of this phenomenon. 2D fluid simulators are standalone applications, while thin-film interference calculations are typically embedded directly in the rendering pipeline through specialized BSDFs. Therefore, we decompose our proposal into two parts: The first describes our intended implementation of a 2D fluid simulator on a spherical topology. The second describes our intended implementation of BSDFs for rendering thin-film interference effects.
How these two elements connect is straightforward. At a high-level, our fluid simulator will propagate a scalar quantity through the fluid's velocity vector field. We will interpret this scalar as the thickness of the film at each point on the spherical bubble surface. As we render each frame, the rendering pipeline will query this thickness map at the points of ray intersection. Using the thickness, the effect of interference for each RGB color channel can be calculated and conveyed to the ray. We intend to employ and modify the rendering pipeline from Assignment 3 for our final project.
The Navier-Stokes equations for incompressible fluids are a system of PDEs that describe the rate of change of a fluid's velocity vector field \( \mathbf{u} \).
\begin{equation} \tag{1} \label{eq: NS1_vec} \frac{\partial{\mathbf{u}}}{\partial{t}} = -\mathbf{u} \cdot \nabla \mathbf{u} + \nu \nabla^2\mathbf{u} + \mathbf{f} \end{equation} \begin{equation} \tag{2} \label{eq: NS2_vec} \nabla \cdot \mathbf{u} = 0 \end{equation}where \( \nu = \frac{\mu}{\rho_0} \) is the kinematic fluid viscosity, \( \nabla^2 \mathbf{u} \) is the vector laplacian, and \( \mathbf{f} \) is a general force vector field encapsulating both internal and body forces. Generally, \( \mathbf{f} = -\nabla p + \mathbf{f_{ext}} \) where \( \nabla p \) is the density-normalized pressure gradient within the fluid and \( \mathbf{f_{ext}} \) is a density-normalized external force like gravity (e.g. \( \mathbf{f_{ext}} = \mathbf{g} \) ). Density-normalized simply means we've divided out both sides of equation by \( \rho_0 \) the fluid density so that the units are compatible. (Navier-Stokes Wiki). Equation two is simply a statement of incompressibility.
In practice, to visualize flow we need to consider the motion of a separate substance within the fluid like smoke or ink. The important feature of this substance is its density in every region of the fluid. This important scalar quantity is transported by the fluid velocity field and diffusion. Its evolution is characterized through the symmetric equation below.
\begin{equation} \tag{3} \label{eq: NS1_scal} \frac{\partial{\Psi}}{\partial{t}} = -\mathbf{u} \cdot \nabla\Psi + \nu \nabla^2\Psi + s \end{equation}Note that the example of density transport is a tautology for the transport of any scalar field. For the purposes of our project, we consider the scalar field to be a thickness map for a thin film object. In other words, we consider the transport of "thickness" within the bubble membrane.
Following [1] we have already implemented a simple planar fluid simulator with various boundary conditions.
Karman Vortex Sheets!
Kelvin-Helmholtz Instability!
Wrapped boundary conditions!
We implement this by maintaining both a scalar density field and a vector velocity field. Both are updated with a diffusion and advection solver. Diffsion is simulated by solving the poisson equation with gauss-seidel and advection is simulated by backtracing. In addition, at every step we ensure the incompressibility constraint by "projecting" the velocity field \( \mathbf{v} \). We do this by taking advantage of the Helmholtz decomposition, which states that any vector field (under sufficient conditions) can be written as the sum of a solenoidal (divergence-less) field and a conservative (curl-less) field:
\begin{equation} \mathbf{F} = -\nabla \Phi + \nabla \times \mathbf{A} \end{equation}where \(\mathbf{A}\) is the vector potential and \(\Phi\) is the scalar potential. Applying the divergence operator to both sides and using the fact that the divergence of a solenoidal field is 0 we have
\begin{equation} \nabla \cdot \mathbf{F} = -\nabla^2 \Phi \end{equation}This gives us a poisson equation that we can again solve using gauss-seidel, giving us the scalar potential. Subtracting off the gradient of this potential gives us the desired solenoidal projection.
Different boundary conditions can be implemented by modifying the gauss-seidel solvers and the backtracing step.
The incentives for creating a fluid simulator for a spherical manifold are two-fold: First, the sphere is the canonical shape of a bubble. Second, it poses an additional challenge that extends from the standard 2D planar fluid simulation for which there are many complete implementations readily available.
We have two ideas for how to go about implementing a spherical fluid simulator. Both approaches closely follow the planar implementation detailed in [1] which realizes a simple yet unconditionally-stable grid-based (Eulerian) fluid simulator. The first approach involves implementing [1], changing the boundary conditions to enable wrapping (since the sphere is a closed surface), and texture-mapping the fluid plane onto the sphere. While this approach is straightforward, it is not physically accurate and would induce warping and unrealistic spontaneous accelerations of the fluid along different regions of the spherical manifold. Additionally, this implementation would demand that flow over the poles be removed. Whether or not this first method is pursued depends on how visibly inconsistent the simulator appears compared to reality, and whether the second approach proves to be too difficult to complete within the allotted time-frame.
The second method involves refactoring the solvers for the Navier-Stokes equations \ref{eq: NS1_vec} and \ref{eq: NS1_scal} so that they are defined in spherical coordinates and account for additional weights/terms born from the spherical metric. Additionally, a careful implementation of new grid boundary conditions will permit fluid flows over the poles. Previous literature [2], [3] describe successful implementation of a fluid dynamics simulator and thoroughly describe resolutions to the singularity issues found at the north and south poles.
Numerical Integration MethodsIn both approaches described previously we appeal to the projection method. This is a two-step process: The first step advances the velocity field as prescribed by the convection and viscosity terms of \ref{eq: NS1_vec}.
\begin{equation} \tag{4} \label{eq: Advance Vel} \mathbf{u^*} = \mathbf{u_0} + \Delta t (\nu\nabla^2 \mathbf{u_0}^2 - \mathbf{u_0}\cdot\nabla\mathbf{u_0}) \end{equation}This invariably results in an updated vector field \( \mathbf{u}^* \) with non-zero divergence, violating equation \ref{eq: NS2_vec}. Therefore, the second step is to numerically solve the following Poisson Equation for some a scalar pressure field \( p \).
\begin{equation} \tag{5} \label{eq: Poisson} \nabla^2 p = \frac{1}{\Delta t} \nabla \cdot \mathbf{u^*} \end{equation}Recall from vector calculus that the gradient of a scalar field is curl-free. Thus we can subtract the gradient of the pressure field from our advanced velocity field to remove the divergence without destroying the vorticity. This removal is known as the projection step.
\begin{equation} \tag{6} \label{eq: Project Vel} \mathbf{u_1} = \mathbf{u^*} - \Delta t \nabla p \end{equation}Consistent with the approaches in the literature, we will solve equation \ref{eq: Advance Vel} using a Semi-Lagrangian backtrace method. This method determines the future fluid velocity at each grid cell by backtracing the current velocity field to a different location on the cell grid and performing a bilinear interpolation on the four cells surrounding the backtraced location.
Also consistent with precedent approaches, we will solve equation \ref{eq: Poisson} using the Implicit Euler method. This method defines a system of linear equations that requires a matrix inversion in order to find \(p\). Instead of inverting the matrix we employ Gauss-Seidel relaxation, an iterative scheme for approximating the inversion.
A sphereical simulator would be nice to have, but may be hard to implement. Alternatives to the sphereical simulator include, as mentioned above, a simple texture mapping from 2D on to the sphere, which we have done here:
Naive texture mapped sphere
An alternative, possibly better hack is to simulate the fluid on a cube. These boundary conditions are fairly easy to implement and may be close enough to be visually appealing. We might suffer from artifacts due to the corners of the cubes.
Thin film interference is fundamentally a result of the wave-like nature of light. Light incident on a transmissive thin film surface undergoes both reflection and refraction. The refracted portion has an additional opportunity to reflect at the second interface of the thin film. The phase difference between these two reflections induces either constructive or destructive interference. The fractional intensities of the total transmitted and reflected light for given wavelength are defined as:
\begin{equation} \tag{7} \label{eq: I_T and I_R} I_T = 1 - I_R = \frac{n_2 \cos(\theta_2)}{n_0 \cos{\theta_0}} \frac {|\beta|^2} {|\alpha|^2 - 2\alpha \cos(\varphi) + 1} \end{equation}where \(\alpha\) and \(\beta\) are related to the Fresnel coefficients as described in [4] and \(\varphi\) is the phase difference between the two reflections.
\begin{equation} \tag{8} \label{eq: phase diff} \varphi = \frac{2 \pi}{\lambda}(2 n_1 \delta \cos(\theta_1) + \epsilon) \end{equation}The phase difference \(\varphi\) is the crux of this visual phenomena. Importantly, it is dependent on the wavelength \(\lambda\) of the incident light, as well as the thickness \(\delta\) of the thin film. This quantity \(\delta\) is precisely what we will query from the 2D fluid simulator at each ray intersection with the bubble. To reemphasize, the thickness of the film is a function of position over the bubble surface for a given instant in time. As the thickness map is a discrete grid, we will use bilinear interpolation to distill the appropriate thickness for ray intersections that do not map precisely onto grid points. Finally, we will be restricting our implementation of thin film interactions to systems containing only two material interfaces. Hence, there are a maximum of three refractive indices \(n_1, n_2,\) and \(n_3\) involved in any calculation.