1. Introduction

The main topic in this work is the implementation, benchmarking and numerical analysis of a combined CFD-Optics solver for problems that involve a fluid component, a rigid body component and an optics component that usually is a laser used as i.e. an optical trap. Optical traps or optical tweezers are devices that use a laser to manipulate an object (usually on a micro-scale). The range of application for these devices is from biological and medical applications to their use in micro technology. As one can easily infer from these applications a simulation with any real world relevance should be able to handle geometries that are non-spherical and can in principle be of any geometrically complex shape. Therefore, in the following sections we show how non-spherical geometry can be included in the proposed solver. A typical application example of optical tweezers with non-spherical geometry is the combination of several differently shaped parts in a micro-assembly process. Furthermore, optical traps can be used to induce rotation of micro rotors or to power micro pumps. The analyzed phenomena in this work appear on a micro-scale to up to 20 microns. In this kind of small scale environment it is advisable to take a look at the validity of the underlying theory of the employed Navier-Stokes fluid solver. In recent studies [1] it could be shown that Navier-Stokes solvers can still produce accurate results in a range of 80 nanometers given certain side conditions regarding the Reynolds number and the domain dimensions. In another publication [2] the FEATFLOW solver has already proven its ability to accurately predict fluid phenomena on a scale of several hundreds of microns. Concerning the general fluid flow configurations we are dealing with low Reynolds number, laminar flow setups for which the FEATFLOW solver also has been shown to produce accurate results for various fluid phenomena [3–5] as well as for the prediction of hydrodynamic forces flow in the context of sedimentation problems or particulate flow simulations [6]. The coupling of the CFD solver with the ray tracing optics solver was done in order to enhance existing numerical simulation procedures for optical traps or tweezers. Often simulations for these kind of problems are done in a vacuum whereas in real applications of the technique the objects are always embedded in a surrounding medium (i.e. a fluid such as water). The downside of this approach is that the dampening friction forces of the fluid are not taken into account which would result in such a simulation overestimating i.e. the velocities of the particles involved. A simple improvement technique for this simulation approach would be to estimate the hydrodynamic forces based on the Stokes model which in case of laminar flows with spherical particles expresses the hydrodynamic friction as $\mathbf {F}_h = 6\pi \eta R v$, where $\eta$ is the fluid viscosity, $R$ the radius of the sphere and $v$ the flow velocity relative to the particle. This method of introducing hydrodynamic forces into the simulation has severe limitations though. The restriction to spherical geometries prohibits the application to advanced configurations where more complex geometries are common and where not only translational displacement is induced by the optic forces, but also rotation. In case of a rotating non-spherical geometry the resulting velocity will be different over the surface of the geometry so that the simple Stokes model is not applicable in this case. Furthermore, a hydrodynamics simulation component that is based on the Stokes model would be a one-way interaction where the fluid force has an influence on the particle, but the particle motion would have no influence on the fluid flow. This is why in our approach a two-way coupling is established that is based on the direct numerical evaluation of the hydrodynamic forces. Apart from the two-way interaction of fluid and laser forces other advantages of using the Fictitious-Boundary Method (FBM) are the ability to include non-spherical geometry and to resolve flow features which allows to gain a better understanding of the involved microfluidics.

2. Methodology

2.1 Coupling of the CFD and ray tracing solvers

At the current time a considerable amount of options exists for the choice of the solver for the fluid part of the problem. Concerning the underlying method of the solver the FEM, Finite Volume Method (FVM) or purely particle-based methods i.e. multi-particle collision dynamics (MPC) could be used for this task. We chose our in-house FEM solver FEATFLOW which has been developed and maintained at our chair for several decades. The FVM would be a fitting alternative approach, but having gathered more experience, deeper insights into the solver and very robust components for computations with embedded geometries like the micro rotors we preferred to conduct the computations with our own solver. The first step in the implementation of the combined solver is to define an interface that allows us to use the result of either solver in the other so that they are coupled in both directions. Thus, we chose an interface that communicates the result of the optics solver to the CFD solver where it is included as an additional force term $\mathbf {F}_o$ in the Newton-Euler equations which are solved in the particle motion step of the CFD solver. In a single step of this coupling approach we first step the fluid solver, then we update the hydrodynamic forces $\mathbf {F}_h$ with the new fluid solution, then in a step of the ray tracing solver the current optic forces are computed so that we can apply the total force from this time step to the particle and update its velocity. This velocity we then set as a Dirichlet boundary condition in the CFD solver which establishes the two-way coupling. Formally, the approach can be described as follows:

2.2 Computational aspects of the coupling

In this particular case of coupling two separate software packages, we were dealing with a MPI-parallel code based on domain decomposition in case of the CFD-solver and in case of the ray tracing solver for optic forces – a serial code. Although there exist several efficient methods for parallelization of ray tracing codes we opted to use the serial state of the software. The reason for this choice being that although plenty of parallelization schemes exist it is hard to find a scheme that truly harmonizes well with the domain decomposition approach used by the CFD-solver which means a scheme that uses the given resources in an effective manner where the used compute resources do not remain (partially) idle while one of the other components does its work. Furthermore, we realized that in a step of the combined solver the CFD-solver clearly requires the bigger share of the total computation time. Quantitatively this would be at least $90\%$ of the total time for a time step, of course depending on the amount of cores used in the parallel computation, but for the usual size of a compute job in the scenarios investigated in this work this relation holds true. Because of this we considered the use of a serial computation step in the combined solver an acceptable circumstance that does not pose a significant performance hit to the overall simulation procedure. In the practical realization this was done by introducing the serial computation step after the hydrodynamics calculation and then distributing the solution of the ray tracing solver to the MPI-processes of the CFD-solver.

2.3 Test case for spherical particles and comparison with experimental data

As a validation case for the combined solver described in the preceding section we decided to start with a simple case of a spherical particle. We conducted an experiment and measurements as basis for the validation case. In the experiment we analyzed the behavior of a particle on a scale of several microns under the influence of optical forces that are applied by an optical trap. We trapped a particle in the optical trap and then released it by turning off the laser. Without the optical trap holding the particle in place it begins to slowly move away from the center of the trap because of several factors of the surrounding fluid medium i.e. Brownian motion. As soon as the particle had moved away a micron from the center of the trap, the laser was again turned on and we measured the time it takes until the particle reaches the center of the optical trap again. This setup we then replicated in the combined CFD-Optics solver and we compared the simulated results with the experimental data (see section 3).

2.4 Ray tracing for non-spherical particles

Apart from configurations with spherical micro-particles we aimed at presenting a method that is capable of handling configurations involving geometrically more complex objects because for many practically relevant applications a restriction to spherical shapes is not sufficient. Although it is possible to describe a lot of geometries used in i.e. micro-assemblies by an analytical description a more general approach to handling geometries is desirable. Because of this we chose to use surface triangulations as geometry description which are employed widely in ordinary (graphics related) ray tracing application and are well studied. Apart from a basic algorithm for optic force computation via the ray tracing approach, we added an extended procedure that makes use of an octree decomposition of the geometry/domain as an acceleration structure. With an octree data structure it is possible to reduce the number of expensive ray-triangle intersection tests by first testing the ray for intersection with the surrounding hull box of the geometry. In case of an intersection the hull box is recursively refined in an octree fashion until we only need to test a minimal number of triangles for intersection with the ray. For our purposes a recursive refinement of 5 levels was found to be most effective, for details we refer to section 3.

2.5 Non-spherical test case: micro rotors

Finally, we put our solver to the test in a realistic application with non-spherical geometry. As with the first test configuration we set up an experiment to that we could compare our simulation results. In the experiment a micro-rotor was placed in a surrounding fluid medium. An optical tweezer was used to induce a rotation of the micro-rotor. We then measured the final angular velocity $\boldsymbol {\omega }_{final}$ of the rotor. The final angular velocity $\boldsymbol {\omega }_{final}$ will vary depending on the power of the laser. Another way to influence the final angular velocity is by changing the rotor geometry. To introduce this geometry aspect into the experimental configuration, micro rotors with different rotor blade angles were used. We measured for different power settings of the optical tweezer the final angular velocities that the rotor would achieve. The process was then repeated for a micro rotor with a different rotor blade angle. For the simulation micro rotors were modeled by surface triangulations and the angle of the rotor blades was adjusted in accordance to the experiment.

3. Results

3.1 Spherical particle test case and comparison with experiments

For the first test case introduced in the previous section we conducted an experiment with a spherical particle to which we applied an optical force. The particle was made of polymethyl methacrylate (PMMA), a material which has a density $\rho _p=1.18\left [\frac {g}{cm^3} \right ]$. The particle was embedded in a surrounding medium of water which was modeled in the simulation as a fluid with density $\rho _f = 10^3 \left [ \frac {kg}{m^3} \right ]$ and viscosity $\mu = 10^{-3} \left [ Pa \cdot s \right ]$. The radius of the particle was 2 $\left [ \mu m \right ]$ and the dimensions of computational domain were $20 \times 20 \times 20 \left [ \mu m, \mu m, \mu m \right ]$. As boundary conditions for the simulation we applied symmetry boundary conditions to all sides of the domain. For the simulation the computational domain was covered with 262144 hexahedral cells which in the FEM-discretization are realized using the Q2P1-Element pair. In accordance to the experiment the particle was placed one micron away from the center of the optical trap which in the simulation coincides with the center of the optical trap. The exact setup of the case is summarized in Table 1. The particle traveled a total distance of 1.0978 microns in the experiment in a time of 0.0292 seconds where the 1 micron distance was passed in 0.0259 seconds. The laser used in the experiment was modeled as a Gaussian ray with a wave length $\lambda$ of 532 $\left [ nm \right ]$, a numerical aperture of 1.2, an index of refraction of 1.33 for the surrounding medium and the power of 0.2 $\left [ W\right ]$. A summary of the optics parameters is given in Table 2. The time step size was set to $\Delta T = 1e-7$ in both simulation components.

Table 1. Fluid and particle parameters

View Table | View all tables in this article

Table 2. Optics parameter

View Table | View all tables in this article

In Fig. 1 the simulation trajectory of the spherical particle over time is shown. The point of arrival at the center of the optical trap is highlighted with a marker. While in the simulation the particle arrived at the center of the trap after 0.0266 seconds whereas this time in the experiment was 0.0259 which gives a relative error of $2.7\%$. Although we were satisfied with the accuracy of the result, we considered improvements with regard to the computation time necessary. The difficulties regarding the computation time can be explained by the micro-scale and the low Reynolds number configuration of the case. These circumstances made it necessary to use a very small time step of $\Delta t = 1e-7$ in order to prevent convergence problems of the solver and inaccuracies in the hydrodynamic force computation. This means that although the simulated time in this case was only $\approx 0.03s$ it took 300000 time steps of the solver to complete the task. This actual simulation time in our parallel computation depends on the number of MPI-processes used and on the number of elements in the computational domain. When using a job size of 16 MPI-processes this would result in a final computation time of over a week. As an improvement to the usual time stepping method we tried a macro time-stepping scheme. The laminar flow property of the case and similar cases of microfluidics allows us to assume that there will not be rapid jumps in the velocity profile of the flow field and consequently of the hydrodynamic forces. This makes it possible to take a larger time step with the current particle velocity for the particle simulation in the Newton-Euler equations and then computing updated hydrodynamics for the new position of the particle. Thus, the macro time step is a step of just the particle solver component of the CFD-solver with a larger time step. Following the macro time step is a series of iterations of the regular combined CFD-Optics solver with the smaller time step in order to compute the current forces acting on the particle and its velocity more accurately and then doing another macro step. In the iterations of the combined CFD-Optics solver the velocity of the particle gets corrected, but another question is how many steps of the regular solver need to be performed before another macro step can be made. We employed two criteria in our simulations, the first one is a fixed number of iterations between the macro time steps. As a second criterion, we kept the solver iterating until the change in the hydrodynamic forces from one time step to the next fell below a certain tolerance. The results of the simulation using the macro time step are shown in Fig. 2. In the macro time step case the particle reached the center of the optical trap after 0.0253 seconds using 77 macro time steps and a number of 200 iterations of the regular solver with a time step of $\Delta T=1e-7$. With the knowledge that the particle moves faster in the beginning of the simulation and slower towards the end, we adapted a smaller macro time step in the beginning of the simulation and a larger one during the flatter movement curve of the particle. The range of the macro step size was chosen adaptively in the range of $[0.0001, 0.0005]$. A total number of 23100 time steps were computed using the macro time step scheme which resulted in a significant improvement with regard to the overall computation time (using the same number of MPI-processes) compared to the standard time stepping scheme where 300000 time step were needed. Concerning the accuracy of the macro time stepping scheme we refer to Fig. 2. We see that in the curve of the macro time stepping scheme the center of the trap is reached in 0.0253 seconds. With regard to the experiment this would be a relative error of $2.3\%$. This however warrants a closer inspection as we know that with decreasing size of the macro time step the solution converges against the solution of the regular time stepping scheme. As the regular time stepping scheme showed a slightly higher value than the experiment, we know the nominally smaller error of the macro time stepping scheme is not due to superior accuracy, but due to the fact that a larger time step is taken for performance reasons, so that the macro time step approach results in a solution that is smaller than the experimental value and incidentally closer. It remains to note that both approaches produced results within $2-3\%$ of the experiment. The macro time stepping scheme can be classified as an explicit Euler stepping scheme and thus shares the advantages and disadvantages of these kind of methods. As can be seen from the trajectory curve the method overestimates the actual velocity in the steep part of the curve which can be compensated by reducing the macro time step size and it also shows a steeper drop of velocity in the dropping part of the curve. Furthermore, despite of the particle remaining at the center of the optical trap with the macro time stepping scheme the particle slightly oscillates around the center of the trap.

 

Fig. 1. Position change induced by optical forces

Download Full Size | PDF

 

Fig. 2. Position change induced by optical forces (macro step and regular stepping method)

Download Full Size | PDF

3.2 Ray tracing solver for non-spherical geometries

Before the routines for force calculation on non-spherical objects with triangulated surfaces could be integrated into the overall concept, the underlying routines first had to be verified. For this purpose, a comparison with wave optical solutions [7] was carried out. In addition, the results were compared with a geometric optics approach in which the intersections between beam and surface were calculated analytically. A comparatively small particle size was chosen to investigate the validity of the method at its limits . The result is shown in Fig. 3. For better comparison, the x-component of the trap efficiency

 

Fig. 3. Force calculation for surface triangulations (spherical particle with $\lambda =1064$ nm, radius $r_P=2.66$ µm $=2\lambda$, $n_P=1.59$ in water ($n_U=1.33$), NA=1.02) for a varying number of triangles showing the wave optical solution (OTT, blue: Nieminen et al. [7]), our own geometric optics solution with analytic intersection computation (analytical, dashed line)

Download Full Size | PDF

 

Fig. 4. Torque force on a rotor for varying $x$-positions (incoming laser ray from $z$-direction

Download Full Size | PDF

3.3 Numerical simulation of micro rotors

An experiment was conducted for the final test case of our combined CFD-Optics solver. In the experiment we triggered rotation of a micro rotor by using optic forces. The corresponding experimental setup is shown in Fig. 5.As lightsource a fiber-laser (Manlight Gevel, 5W output power) with a wavelength of 1064 nm was used. The optical trap was established and manipulated with help of a spatial light modulator (Holoeye Photonics). The rotor is created with help of two photon polymerization process [8]. As Polymer Femtobond (Laserzentrum Hannover) was used. In the measuring process, the micro rotor was embedded in a surrounding medium of a water/surfactant mixture inside a movable object holder. The properties of the mixture are almost identical with those of water and its viscosity was determined to be $\mu$ = 1.1 $\left [ mPa \cdot s \right ]$ and in the corresponding simulation we set the fluid density to $\rho _f = 1 \left [g/cm^3 \right ]$. The material of the micro rotor as well had a density of $\rho _{rotor}= 1 \left [g/cm^3 \right ]$. We used rotor blades that are attached at different angles in order to get a dependency on the geometry into the configuration.

 

Fig. 5. Experimental setup and an exemplary picture of a microrotor

Download Full Size | PDF

A real micro rotor is shown in Fig. 5(b) and a surface triangulation of a micro rotor used in the simulation is shown in Fig. 6.We used attachment angles of $45^{\circ}$ and $60^{\circ}$ in the experiment as well as in the simulation. The dimensions of the rotor are $22.7\times 21.0\times 11.0 \left [\mu m, \mu m, \mu m \right ]$ so that it could be embedded in a computational domain of size $30 \times 30 \times 15 \left [ \mu m, \mu m, \mu m \right ]$. The micro rotor was placed in the center of the computational domain and on all domain boundaries we prescribed symmetry boundary conditions in the CFD-solver. The computational domain was covered with a hexahedral hierarchical mesh which on the highest refinement level had a resolution of about one million elements. The fictitious boundary approach that we use is a variant of the fictitious domain method where there is only one mesh used for the fluid domain and the embedded solid domain. Mesh convergence and mesh independence was established in test calculations prior to the benchmark. In order to determine our target mesh resolution that provides mesh independent results, we at first computed the final angular velocity on a mesh $\mathcal {T}_1$. We then refined mesh $\mathcal {T}_1$ to obtain $\mathcal {T}_2$ and prolongated the solution to the refined mesh $\mathcal {T}_2$ which is easily doable in our setup with hierarchical meshes and regular refinement. We then resumed the computation on the mesh $\mathcal {T}_2$ and observed whether there was a significant change in the value for the angular velocity. We then chose the mesh resolution so that from one refinement level to the next there was no significant change, the final mesh is depicted in Fig. 7. Thus, an accurate resolution of the rotor geometry could be achieved which is crucial to the accuracy of the hydrodynamics force computation and in capturing the influence of the different rotor blade angles. We show the results of the simulation compared to the experimental data in Fig. 8 where the final angular velocity $\boldsymbol {\omega }_{final}$ depending on the power of the laser is displayed. We observe that there is a linear relation between the laser power and the final angular velocity $\boldsymbol {\omega }_{final}$. We can see that this linear relation is reproduced correctly by the simulation. Furthermore, the relation between the $60^{\circ}$ rotor and the $45^{\circ}$ rotor where the $60^{\circ}$ rotor consistently attains a higher final angular velocity for all power settings of the laser is predicted correctly by the simulation. When looking at the value of the final angular velocity, we see that the simulation computes a higher value than was recorded in the experiment. In the experiment the rotor was guided by a support structure. This may give rise to further frictional forces which were not considered in the simulation and led to a deceleration of the movement. Another possible explanation for the discrepancy is that the surface triangulation of the rotor was not exactly matching the rotor geometry in the experiment. It is worth noticing though, that the computed values are consistently higher in the simulation and the shear relation between the $60^{\circ}$ rotor and the $45^{\circ}$ curves show a good agreement. We can conclude from these simulations that in this form a highly accurate numerical solution could not be computed, but that qualitatively a good predictor for the relevant phenomena could be computed that was able to reproduce geometry dependent correlations.

 

Fig. 6. Computational models of the rotors (left: $45^{\circ}$, right: $60^{\circ}$)

Download Full Size | PDF

 

Fig. 7. Computational mesh used for the micro rotor calculation

Download Full Size | PDF

 

Fig. 8. Comparison between experiment and simulation

Download Full Size | PDF

4. Conclusion

In the presented work we developed a combined CFD-Optics solver for configurations with spherical as well as non-spherical geometries. The coupling between the CFD-Solver and the ray tracing optics solver was done by injecting a ray tracing step after the hydrodynamics computation step of the CFD-Solver and then adding up the forces. The combined solver was put to the test in different benchmark configurations for which we generated experimental data. As a consequence of our findings in the optic trap setup with a spherical particle we enhanced our simulation procedure by a macro time stepping scheme which could greatly reduce computation time while maintaining good accuracy.

Furthermore, we pursued our second goal which was the inclusion of non-spherical geometries in order to simulate practically relevant configurations. For the optics component of the combined solver an enhanced ray tracing solver based on the octree algorithm was implemented. This approach enabled us to calculate optic forces for surface triangulations with a large number of triangles in an efficient and accurate manner. As test case for non-spherical geometries we conducted an experiment with micro rotors. In this benchmark configuration the solver proved to be able to predict the relevant features and phenomena of the case. We thus have a simulation procedure at hand that can predict the influence of geometric details on the outcome of the simulation. In previous studies we have shown that our solver is capable of accurately simulating the effects introduced by particle size and density [4,9], viscosity and the rheology model [2] which were done without the influence of optic forces. Further studies given corresponding experiments with participation of optical forces and fluids of different viscosity and rheology as well as varying particle parameters could be a natural extension to the current state of our work.