Geometric optimization of microreactor chambers to increase the homogeneity of the velocity field

In this work microfluidic flow-through chambers are investigated. They are filled with magnetic nanoparticle (MNP) suspension in order to facilitate enzymatic reactions. The enzyme is immobilized on the surface of the MNPs. These reactions have been found to be flow rate dependent. To overcome this issue various chamber geometries have been examined and optimized geometries have been designed and tested experimentally. The investigation is supported with dedicated CFD simulations using the open source software OpenFOAM. The paper presents the theoretical background and the results of the simulations. The simulations have been verified with measurements and these too are presented in the paper.


Introduction
The microfluidic application of immobilized enzymes offers novel possibilities in diagnostics, synthetic, and analytical applications [1]. Catalytic biomolecules such as enzymes covalently attached to magnetic nanoparticles (MNPs, see figure 1(a) are applicable as a biocatalyst in microfluidic systems in a variety of ways. Biofunctionalised nanoparticles thus showing biocatalytic activity are able to flow in microfluidic cartridges together with the liquid. Magnetic microreactor cartridges (such as Spinsplit MagneChip) are microfluidic devices consisting of several reactor chambers ( figure 1(b)) within which the MNPs are entrapped using an external magnetic field. The accumulated particles form a thin biocatalyst layer in the microchambers while the magnet placed beneath the chamber keeps the particles in the chamber volume despite the liquid flow through the chambers (figure 1(c)). The inhomogeneous phase packed-bed microreactor formed in this way creates a unique opportunity to develop modular microsystems with the ability to allow flexible variation of the biocatalysts [2]. Performing enzymatic reactions in such microreactors has many advantages due to their reduced (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
size. The small reactor volume requires a smaller amount of expensive reagents which are often available only in small quantities [3]. The measurement of a high number of samples can also be facilitated, e.g. by microreactor parallelization or by automation, which also reduces the overall measurement time and manpower compared to traditional batch chemistry methods. Once the S substrate solution is introduced in the microreactor chip, the enzymes transform the substrate molecules to P product molecules (see figure 1(b)) according to the kinetics of the given enzyme. Such a reaction, i.e. the change of product concentration over time, may be detected by an appropriate detector. The ability to have precise environmental control in the microliter scale enables safe, clean and highly reproducible reactions and parameter sweeping with repeated measurements [1]. Therefore a given enzyme can be investigated under various conditions, resulting in enzyme kinetics data ( figure 1(d)). High surface-to-volume ratio, and precise temperature control are the further benefits of small volume, as well as reduced cross-effects during enzyme characterisation [4][5][6].
It is claimed that the enzyme-substrate reactions are dependent on the flow rate in flow-through reactors [7]. The dependence was also shown in microfluidic reactors [3,4]. Different theoretical models of the flow rate dependence are given in [3,8]. In spite of many advantages of using microscale reactor chips, the flow rate dependence seems to be a key limitation of such devices as inhomogeneous velocity distribution in the chambers (see figure 1(c)) may result in inhomogeneous reaction rates and unreliable measurement data, as enzymes in different locations are subjected to different flow conditions. The goal of the paper is to create a microfluidic chip in which the flow rate dependence of an enzyme-substrate reaction can be measured precisely. The geometry is optimized based on our previous microfluidic chip, in which enzymatic reactions have been successfully performed [1] (see figure 2). In the current work the original chamber geometry was modified in order to address the above issues i.e. to achieve more homogeneous velocity distribution. The chip under invest igation consisted of serially connected microchambers. Note that the reason for realizing the chambers as channel widenings is that this means more MNPs can be kept against the flow with a single magnet. More suspension means a higher yield in a reaction, which is advantageous (e.g. the product concentration will be higher, thus it can be measured more quickly or more precisely).
First, we accomplished a CFD simulation to model the velocity field of the MNP filled chamber used in the original measurements in [1]. After this we designed new chamber geometries and compared the homogeneity of their simulated flow field to the original geometry. This design process is presented in detail in our recent work [9,10]. In the following sections we only show the simulations of the original and two new chamber geometries, which have been selected for the measurements. Results of simulations with further geometries from [9] can be found in the appendix. After the simulations the measurements are shown with the new chamber geometries. Finally, we compare our simulation model to the measurement data.
Enzyme reactions were successfully performed in the microreactor chip; however, the analysis of these measurements is not in the scope of this current paper.

Flow type
The flow is considered laminar at the used flow rates (Q 200 μl min −1 ). The fluid in the measurements was water and isopropanol (IPA), both of which are incompressible. We investigated a steady-state flow.

Porous media
The MNP suspension is modelled as porous media due to the previous measurement results [1]. In the porous media the gradient of the pressure is proportional to the flow velocity: where μ is the dynamic viscosity, and D visc is the viscous resistance. The viscous resist ance generally can be direction dependent, but in our case it is assumed to be isotropic due to the spherical shape of the nanoparticles. The magnitude of the resistance was determined from the results of the previous measurements [1]:

Comparison of the geometries based on velocity
To compare the different chamber geometry cases we compute the following quantities in each simulation. The volume  weighted average velocity of the MNP suspension is calculated as:v where V sus is the total volume of the MNP suspension, while V i and v i are the ith computational cell's volume and velocity in the MNP suspension, respectively. The volume weighted standard deviation is then calculated as In the CFD simulations the volume weighted average velocity and deviation in the MNP suspension are calculated in each case. Our goal is to minimize the ratio of the deviation to the average velocity.

CFD software
We used the open source software OpenFOAM 5.0 for the simulations [11]. The simulations were run on an Intel Core i5 4 core CPU with 4 GB RAM under 64 bit Ubuntu GNU/Linux (kernel: 4.4.0-87-generic). The post-processing software was the open source software ParaView 5.4 [12].

Geometry
Instead of modelling the whole microfluidic structure, we model led only one microchamber with the inlet and outlet channels. In this case it was enough to model the quarter of the chamber volume using the two symmetries of the geometry and the flow. We investigated different chamber geometries with simulations.

Mesh
The mesh was created with the blockMesh and the snap pyHexMesh programs, which both belong to OpenFOAM. For the mesh generation we also used the m4 macro processor to parametrize the geometry and the mesh. For the simulations prepared after the measurements, first we drew the MNP shape based on the pictures obtained during the measurement. This was done in the open source software FreeCAD 0.16 [13]. Using the software we imported the selected measurement picture and then, based on the picture, we drew the shape of the chamber and MNP suspension precisely. The shapes were exported to stl files, which contained the triangulated surface of each shape. We then built the mesh with snappyHexMesh using the stl files.

Simulation settings and boundary conditions
We selected OpenFOAM's simpleFoam solver, which is capable of modellin incompressible, steady-state flow. The solver was modified to also calculate the volume of the computational cells from the MNP suspension for the volume weighted average velocity and deviation calculation (see equations (3) and (4)). The MNP suspension was model led as a porous media based on equation (1). Porous media settings were done in a fvOptions file using explicitPorosi tySource to the MNP suspension, which was represented as a cell zone (certain defined volume) in the simulations.
We set an homogeneous velocity inlet boundary condition for the inlet, which was calculated from the applied flow rate. At the walls the velocity was 0, while at the outlet the gradient of the velocity was 0 (zeroGradient boundary condition). The pressure was 0 at the outlet and zeroGradient at the walls and at the inlet.
For the laminar simulations we have to give the fluid kinematic viscosity, as a material property. We set the kinematic viscosity of the water to ν = 8.928 · 10 −7 m 2 s −1 at T = 25 °C, and the isopropanol's viscosity to ν = 2.576 · 10 −6 m 2 s −1 at the same temperature (ρ IPA = 781.2 kg m −3 and the dynamic viscosity is µ IPA = 2.012 · 10 −3 Pa s at T = 25 °C [14]).

Simulation results
3.5.1. Original geometry. The chamber geometry used in the original measurement was cylindrical. The height of the structure was uniformly h = 110 μm. The diameter of the chamber was D = 3600 μm, which means that the volume of one microchamber was approximately V ≈ 1.12 μl. The applied flow rate in the simulation was set to Q = 28.61 μl min −1 , one of the values which was used in the measurements. The fluid was water, the material properties were set at T = 25 °C. The MNP suspension totally filled the chamber. The simulated flow field is shown in figure 3. As expected, the velocity in the centre parts of the MNP filled chamber was higher than at the side parts. The volume weighted average velocity and deviation were determined in each simulation with the different meshes. The results are presented in table 1. The ratio of the calculated velocity deviation to the average velocity is 61.80% in the case of the finest mesh.
In the following sections we investigate the new enhanced geometries based on [9,10] to decrease this ratio between the deviation and the average velocity, i.e. to increase the velocity field homogeneity in the MNP suspension. Velocity magnitude in the middle of the chamber at the symmetry plane with the finest mesh. The fluid flows from the left to the right. The maximum of the colorbar is decreased to 3.5 · 10 −3 m s −1 to highlight the velocity field in the chamber. The MNP suspension totally fills the chamber, its shape is highlighted with the dashed black line.

3.5.2.
Chamber with three inlet and three outlet channels. In our previous work we investigated cylindrical chambers with more than one inlet and one outlet channel [9]. Simulations were created for a chamber with three inlet and three outlet channels, where the angle between the adjacent channels were set in the different cases to α = 30 • , 45 • , 60 • and 75 • . A chamber with three inlet and one outlet channel was also investigated, where the angles between the adjacent inlet channels were set to the previously listed angles in the different cases. A chamber with five inlet and five outlet channels was also simulated. These cases can be found in the appendix or in [9]. From these cases the best candidate was a chamber with three inlet and three outlet channels, where the angle between the adjacent channels was 30 • . Based on these simulation results, we designed a new geometry where the structure height and the channel width were the same as in the original case, but the chamber diameter is decreased to D = 2800 μm, and therefore the volume of one chamber is V ≈ 0.68 μl [10]. The three inlet channels branch from a main inlet channel before the chamber (the same is true for the outlet), see figure 4. The flow rate was again Q = 28.61 μl min −1 . The simulated flow field in the middle of the chamber is shown in figure 4. The volume weighted average velocity and deviation with the different meshes are presented in table 2. The ratio of the deviation to the average velocity is 29.24% with the finest mesh, which is approximately 2.1 times lower than in the original case.

Chamber as a widened channel.
We also investigated a case where a chamber was similar to a widened channel [9,10]. The chamber had a maximum diameter approximately of D = 2300 μm at the centre. The structure height, channel width and flow rate are the same as in the previous simulations. The border of the MNP suspension is assumed to be as in figure 5. The volume weighted average velocities and deviations with the different meshes are presented in table 3. The ratio of the deviation to the average velocity is 26.59% with the finest mesh, which is approximately 2.3 times lower than in the original case.

Goal of the measurement
The main goal of the measurement was to investigate how the new chamber geometry is being filled with the MNP suspension. The pressure drop of the structures at different flow rates were measured when the chambers were not filled with MNPs and when they were filled with MNPs. After this, the simulation model was compared to the measurement data based on the pressure measurements.

Microfluidic chip
Two new different microfluidic structures were designed and manufactured based on the simulation results. Each structure consisted of six serially connected chambers. In the first structure the chambers had the widened channel shape, which was investigated in section 3.5.2. The second structure had chambers with three inlet and three outlet channels (section 3.5.3). The layouts of the two chips are shown in figure 6. In the following sections the first chip in figure 6 is noted with C1 and the second with C2. The chips were created using the PDMS molding technique (PDMS is from Dow Corning Inc., Sylgard 184). The molding master was created from SU-8 photoresist and contained both structures. The PDMS was poured into the master and remained at rest until the crosslinking (∼one day at room temperature). All chips were designed to be matched to a rectangular glass plate with the dimensions 75 mm · 25 mm.

MNP suspension creation
The core of the nanoparticle is made of magnetite (Fe 3 O 4 ). The surface of the nanoparticle is covered with a SiO 2 layer. Then PAL (phenylalanine ammonia-lyase) enzyme is immobilized onto the surface (for more detail, see [15]). The MNPs were created by SynBiocat Ltd (Budapest, Hungary).
To create the MNP suspension for the C1 chip, initially we added 2.6 mg MNP powder and 4.2 mg PEG 4000 (polyethylene-glycol) in a vial. The PEG was important to avoid the aggregation of the nanoparticles. Then, we poured 600 μl    ultra-clean water and 200 μl isopropanol into the vial and the mixture was sonicated for 30 min at T = 30 °C. After that the suspension was diluted with an additional 1800 μl water and 400 μl isopropanol which means that the total added volume of water was 2400 μl and the total volume of isopropanol was 600 μl. Finally, the vial containing the diluted MNP suspension was continuously shaken with a shaker until usage.
In the case of the measurements with the C2 chip (chambers with the three inlet and three outlet channels) we used slightly different settings. The core and size of the MNPs were the same as in the previous case, but they were coated with PfXAL enzyme to perform different enzymatic measurements than in the previous case (these measurements are not in the scope of this paper). 2 mg MNP and 2 μl PEG 400 was poured into a vial with 200 μl ispropanol and 800 μl water. After 30 min sonicating at T = 30 °C the suspension was diluted with an additional 3200 μl water and 800 μl isoproranol. The suspension was then continuously shaken with the shaker until usage.

Measurement with the C1 chip
We used the MagneFlow instrument supplied by Spinsplit Ltd (Budapest, Hungary) in order to take the measurements. The instrument is capable of accommodating microreactor chips and actuating magnets at predesigned positions to anchor MNPs within the chambers of the reactor chip, dispensing fluids into the reactor chip and measuring any pressure drop that develops in the reactor. The schematic of the measurement layout is shown in figure 7. Using the MagneFlow instrument, one can connect more inlet sources to the chip. The software controlled syringe can dispense fluid through the selected valve via the valve control unit (also software controlled). In our case we used different inlet pipes for the MNP suspension, isopropanol and water.
Initially, we connected every pipe to the chip and the MagneFlow instrument, as shown in figure 7. Then, with the software controlled syringe pump the pipes were filled, as well as the microfluidic structure, with isopropanol, ensuring gas-bubble free loading of the reactor volume. Bubbles are disturbing, because they can prevent the fluid flow and also modify the shape of the MNP suspension. The MagneFlow device has a temperature controller next to the chip, which was set to T = 25 °C, while in the room the temperature was about 23 °C-24 °C. Using a camera (also part of the instrument), which was attached under the chip, we continuously monitored the chip with the chambers. The pressure was measured with the pressure sensor of the instrument which generates a   Each chip consists of six serially connected chambers. The first structure has chambers with the widened channel shape while the second has chambers with three inlet and three outlet channels. The angles here between the adjacent channels are α = 30 • . In the following sections the first chip is denoted by C1 while the second is denoted by C2.   proportional voltage signal with the pressure. The continuous flow was applied with the PC controlled syringe. As can be seen in figure 7, the chip had three different inlet ports. First the pressure drop was measured at the different flow rates without MNP. We filled the structure and the pipes with isopropanol. Using this liquid we were able to measure the pressure drop at the different flow rates, and we did not notice any gas bubbles. At every flow rate we measured the pressure for at least 200 s. The reason for this long measurement time was that during the measurement the pressure did not increase sharply. It started to increase slowly until it reached its final value, which needed approximately two to three minutes (see figure 11). The error of the measurement was assumed by the fluctuation of the final pressure value. Note that, although the pressure sensor can have an error between ±3.45 kPa based on its datasheet, assuming that this error was constant during the measurements, we did not have to deal with it. The reason for this is that we were interested only in the pressure difference between the different measurements (see later).
After these measurements had been taken, we filled the chambers with MNP. To achieve this, we dispensed the MNP syringe to the chip, and at the same time we switched on the magnets over the chambers. Initially the magnets were not close to the chip, therefore they did not have any effect. The switching procedure means that they were moved next to the chambers. This procedure was also PC controlled. First we switched on the magnet over the last chamber along the flow. As the magnetic field anchored the nanoparticles, the chamber was filled up with MNPs. After it was filled we switched on the magnet at the previous chamber, one after another. Using this procedure we could fill up all six chambers. The photo of the MNP-filled chambers after the filling procedure is shown in figure 8. After accomplishing the filling procedure we measured the pressure drop at the same flow rates as in the previous measurements. Using isopropanol as the liquid, we experienced no bubbles, and the nanoparticle suspension shape has remained constant during the whole measurement. We have taken photos after every flow rate measurement with the camera in order to verify this.
In line with expectations, the measured pressure drop seemed to be linearly proportional to the flow rate in both cases, see figure 9. The pressure drop values were higher in the second case, because the MNP suspension represents an extra resistance to the fluid flow. As shown in figure 9, we fit a line with zero offset to the pressure difference dependence against the flow rate. The fitting was based on the theoretical assumption that the MNP suspension can be treated as porous media.
Finally, we estimated the structure's uniform height with a microscope (microscope: Olympus BX-51M; objective: MPlanFL-N 10x with a WHN10X eyepiece). We measured the objective movement between two states: at the first state we adjusted the upper side of the structure in focus, and then, in the second case, we set the lower side. We measured 12 different sections of the chip. Although the SU-8 mask was designed with a height of 100 μm, the measured height of the channels and chambers in the structure was 118 ± 5 μm.

Measurement with the C2 chip
Measurements have also been performed with the C2 chip (with chambers with three inlet and three outlet channels). The picture of the MNP filled structure is presented in figure 10. As one can see, the MNP suspension filling of the chambers was different from what was previously expected. The MNPs filled most of the chamber but, in addition, they have also filled some parts of the outlet channels. The pressure drop of the system was measured at different flow rates. These measurements were taken in a similar way to those in the previous case, i.e. first the pressure drop was measured with isopropanol without MNPs and then it was measured with the MNP-filled chambers. The only different setting compared to the previous measurements was that in this case the output signal of the pressure sensor was monitored with one of the analogue inputs of the MagneFlow instrument. In this way the output signal, and therefore the pressure, could be logged over the measurement time. The graph based on this log file is shown in figure 11. The pressure drop at the different flow rates can be identified from this graph easily. The pressure at the different flow rates in the 'MNP free' and 'MNP filled' cases is shown in figure 12, while the difference between the two graphs is also shown there. The pressure values were identified with a moving average method with 50 samples, while the error was calculated as the standard deviation of the current sample selection. This method gave a more precise determination of the error than the previous error estimations with the other chip.
Finally, we measured the structure height at 12 different points in the same way as for the previous chip. The average height was 118 ± 3 μm.

C1 chip
Next we estimated how the MNP suspension filled the chambers. First we focus on the chip with the widened channel shaped chambers. The exact chamber fillings with MNPs were not straightforward to determine because the border of the chambers cannot be identified from the measurement pictures. The reason for this is that the liquid-filled PDMS structure is transparent (see e.g. figure 8). To make the structure more  -shaped black area). The PDMS structure is transparent; the white circles under the chambers are the neodymium magnets (highlighted with a white circle at the magnified picture) while the metal area with the six holes is part of the chip holder. Note that the structure is flipped vertically compared to the lower structure in figure 6 because the camera saw the structure from below.
visible we filled them with air after the measurements, which highlighted the channel and chamber borders. Of course, in this case there were no MNPs in the chambers. We denoted the borders with red lines based on the air-filled case picture, and then we merged these lines with one of the original measurement pictures. The result of this process is shown in figure 13, where the MNP-filled chamber five is shown with its borders, counted from the inlet. It can be seen in figure 13, that the given chamber was not totally filled with MNP and this was also true for the other chambers. In all chambers the MNP suspension had a leaf shape, but its diameter was smaller than the chambers. This is an undesired effect, because some parts of the fluid bypasses the suspension in the chamber, which decreases the efficiency of the enzyme-substrate reaction, as parts of the substrate bypasses the suspension, i.e. the quantity of the generated product will be smaller.
To investigate this partially-filled case, next we tried to simulate the measurement in OpenFOAM. For the simulation we imported the picture of figure 13 into the FreeCAD software and then we drew the shape of the MNP suspension based on the picture, see figure 14. Although we already had the designed chamber shape from the previous simulations (see section 3.5.3), we also drew the chamber based on the picture to get more precise results ( figure 14). In the CAD software both the chamber and MNP suspension was extruded perpendicular to the image plane surface. The faces then were exported to stl files. The objects drawn in FreeCAD were bigger than the real structures, therefore the stl files had to be scaled down. This was achieved by a simple script which was written in the AWK programming language. In the simulations we set the chamber's height to half the measured h = 118 μm value, i.e. it was set to 59 μm using the flow's vertical symmetry. Although from the pictures no information could be identified about the vertical filling, we assumed that the MNP vertically filled the chamber totally. Note that-in contrast to the previous simulations-here we did not use the vertical symmetry plane, because the MNP suspension was not axially symmetric. The mesh was generated with the snappyHexMesh program using the modified stl input files. The mesh consisted of 347 544 elements.
It can be seen in figure 14 that the simulated flow field had high velocities in the gap between the MNP suspension and the chamber border. To estimate the amount of the flow which bypassed the suspension, we investigated the flow rates through an artificial internal surface. This surface had a normal in the x direction (in which the inlet and outlet channels were directed), and is shown by a dashed green line in figure 14. We calculated the flow rate through the intersection of this surface and the gap, and we also calculated the flow rate through the intersection of the surface and the suspension. The results showed that 51.11% of the total flow rate flowed through the gap, while the remaining flowed through the MNP suspension. It can be seen that even this small gap highly reduced   the flow rate through the suspension, therefore reduced the efficiency of the device. The calculated ratio of the volume weighted deviation to the volume weighted average velocity in the MNP suspension was 50.68%. Although this value was smaller than in the original measurements, it was rather high to be acceptable for the further measurements.
We were interested in comparing the simulations and the measurement data, therefore we checked the pressure drop in the simulation. It was p filled ≈ 1237 Pa. We then repeated the simulation with the empty chamber, i.e. when there was no MNP suspension. Using the same flow rate the calculated pressure drop was p empty ≈ 651 Pa. The difference between the pressure drops in the two cases were caused by the suspension, and its value was p diff = 586 Pa for the simulated chamber. This means that for the whole chip, which consisted of six chambers, the pressure difference was 6 · 586 Pa = 3516 Pa based on the simulation. This result corresponds with the measured pressure of 3.59 kPa at the flow rate Q = 100 μl l −1 (see figure 9). This confirms that the porous model of the MNP suspension can be used in the simulations. The reason for the slight difference between the simulation and measurement results could be the fact that the shape of the MNP suspensions in the other chambers were not exactly the same as in the investigated chamber five.

C2 chip
The MNP filled second chip has been also investigated with a simulation. This was done in nearly the same way as for the previous chip. Although the filling of each chamber in figure 10 was not exactly the same, the MNP suspensions in the different chambers have roughly the same shape. Bigger differences in the MNP quantity can be found at the outlet channels next to the chambers.
For the simulation we focused on chamber four, counted from the inlet. Using the air-filled image after the measurements, the chamber layout and the MNP filled picture, the photo of the chamber with the highlighted chamber walls and MNP suspension are shown in figure 15. Unfortunately, the edges of the chamber and the suspension were not as sharp as in the previous case, which means that there is a bigger uncertainty in their exact positions. The CFD simulation model was created based on this picture with FreeCAD and snap pyHexMesh, the simulation results can be seen in figure 15.
The fluid was isopropanol and the flow rate was set to Q = 131 μl min −1 . Based on the simulation, only 26.68% of the fluid flowed through the MNP, while the remaining part bypassed it through the gap. This is undesirable as it reduces the efficiency of the microreactor in the case of an enzyme-substrate reaction.
The ratio of the volume weighted velocity deviation to the volume weighted average was 51.36% which was not much lower than in the original measurements. The pressure difference between the MNP-free and MNP-filled cases was 484.3 Pa for one chamber at Q = 131 μl min −1 in the simulations. This means that for the whole structure the total pressure difference was 6 · 484.3 Pa ≈ 2905.8 Pa. This value is acceptable compared to the fitted line in the measurements, which gives 3.18 kPa at the investigated flow rate (see figure 12). It should be noted, however, that the difference between the measurement and simulation can originate from several facts. First, the border of the MNP suspension and the chamber wall could not be identified as clearly as in the C1 chip case. The second reason for the difference could be the slightly different filling of each chamber. Third, the different type of enzyme and PEG solution could also modify the viscous resistance of the modelled MNP suspension.

Discussion and future work
The two investigated chips were successfully filled with MNPs. In both cases the different chambers were filled with MNPs in a roughly similar way, however the suspension shapes were different to what was previously expected. Based on the measurement and simulation results, we plan to redesign the investigated microfluidic chip with reduced chamber sizes for the further measurements. In this way, we expect to  avoid the incomplete filling of the chambers by using the same MagneFlow device with the same magnets. Another workaround would be to use stronger magnets.
These facts imply that there is a need for a simulation model for the MNP filling procedure. We aim to achieve this by using the multiphase particle-in-cell method (MP-PIC). This is already implemented in OpenFOAM in the MPPICFoam solver based on [16]. In our case, the effect of the neodymium magnets on the particles should also be considered. The magnetic field can be calculated with OpenFOAM's magn eticFOAM solver and then the magnetic force calculation will have to be merged into the MP-PIC solver. By achieving these tasks we try to simulate the filling process of the MNPs. The main benefit of this model would be to give a prediction about the final shape of the MNP suspension. In addition, this method could give a more accurate model for the MNP suspension than the current porous model. Accurate modelling of the MNP suspension shape will help us design of better chips more easily than the current trial-and-error way.

Summary and conclusion
In this paper we have investigated the velocity field of enzymatic microreactors in microfludidic structures with MNPs. Enzyme-substrate reactions are performed in the microreactors (chambers), where the reaction parameters have a flow rate, i.e. velocity dependence. To investigate this dependence more precisely, our goal in this paper was to design a simple chamber geometry in which the velocity field is more homogeneous than in the simple cylindrically-shaped chambers with one inlet and one outlet channel. To achieve that objective we investigated several new geometries.
The different chamber geometries were simulated with the open-source software OpenFOAM, with each case on four different meshes. The homogeneity was investigated by calculating the ratio of the volume weighted velocity deviation to the volume weighted average velocity in the MNP suspension. First, the original cylindrical chamber with one inlet and one outlet channel was simulated, where the ratio of the deviation to the average velocity was 61.80%. Then several new geometries were investigated, from which the two best candidates were fabricated and tested experimentally. The first chip consists of chambers with a widened channel shape while the second has cylindrical chambers with three inlet and three outlet channels, setting α = 30 • between the adjacent channels. The calculated ratio of the deviation to the average velocity was 26.59% in the widened channel chamber and 29.24% in the chamber with three inlet and three outlet channels.
The two fabricated chips were tested by filling them with MNP suspension. In each chip the chambers were filled roughly similarly, but the filling was incomplete as a small gap remained between the chamber wall and the MNP suspension. This was undesirable from the point of view of the microreactor efficiency. The pressure was measured at different flow rates with the MNP-free and MNP-filled chips to compare the simulation model to the measurements. In both cases the simulation and measurement results were in good agreement.
The fabricated chips were capable of performing enzymesubstrate reactions, but the shape of the MNP suspension was different from the previously expected one. In further work the chamber sizes should be decreased to avoid the gaps between the suspension and the wall. It is planned that the filling procedure of the chambers will be simulated to predict the shape of the MNP suspension.  Note that the magnet was precisely over the centre of the chamber; the apparent shift in the picture was caused by the perspective distortion. Bottom: the simulated velocity field of the isopropanol flow in the middle of the structure. The flow rate was set to Q = 131 μl min −1 .      calculated magnitude of the velocity field for the α = 45 • case is shown in figure A1. Each geometry was simulated on four different meshes, as in the original case. The results for the average velocity and deviation are shown in figure A2.
Results for the simulations with the finest meshes are collected in table A1. The ratio of the velocity deviation to the average velocity is the lowest in the α = 30 • case, it is only 31.50%, therefore approximately two times lower than in the original case.

A.2. Simulation with three-one channels
In this section, chambers with three inlet and one outlet channels are simulated. The angles between the adjacent inlet channels are set to the same α = 30 • , 45 • , 60 • and 75° values as in the previous case. The outlet channel is positioned symmetrically to the symmetry plane. All cases have been simulated on four different meshes. Results for α = 30 • can be seen in figure A3. Simulation results for the average velocity and deviation are shown in figure A4. Results with the finest meshes are collected in table A2. The ratio of the velocity deviation to the average velocity is calculated and found to be slightly less than at the original geometry (lowest value is 52.2% at α = 30 • ).

A.3. Chamber with five inlet and five outlet channels
Next, the flow field of a chamber with five inlet and five outlet channels is presented. The angle α between the adjacent channels is set to α = 30 • . To reach the same flow rate, the velocity at the inlet is decreased by one fifth as in the original case. The magnitude of the velocity field in the middle of the structure is shown in figure A5. The calculated average velocity and deviation on the different meshes are presented in table A3.
The ratio of the velocity deviation to the average velocity is 32.29% with the finest mesh.