You are viewing the documentation for an older COMSOL version. The latest version is available here.
PDF

Annular Ultraviolet Reactor with Particle Tracing
Introduction
In this example, the sterilization of water in a simple ultraviolet (UV) water purification reactor is modeled using a combination of optical ray tracing, computational fluid dynamics (CFD), and Lagrangian particle tracking methods.
The model consists of three studies. In the first study, the fluence rate in the annular region surrounding a UV lamp is solved for, using the Geometrical Optics physics interface. The fluence rate represents the amount of radiation that a small spherical detector would absorb per unit time, divided by the cross sectional area of such a detector. It has the same units as irradiance (SI unit: W/m2), but is a volume quantity and can be defined at any arbitrary point in space, not just on surfaces.
In the second study, the velocity profile of water flowing through the reactor is solved for, using the Turbulent Flow, k-ε interface. The reactor in this example is a U-shape with inlet and outlet tubes running perpendicular to the orientation of the lamp.
In the third study, particles are traced through the reactor using the Particle Tracing for Fluid Flow interface. The particle velocity is based on the velocity of the water solved for in the second study. As the particles move through the reactor, the accumulated dose along their trajectories is computed, based on the fluence rate solved for in the first study.
The results show that all particles arriving at the outlet have been exposed to a certain minimum dose. A small number of particles, mostly those passing very close to the lamp, are subject to a substantially higher dose. This model could be used as a starting point to develop more energy-efficient UV water purification reactors by seeking designs with a more uniform accumulated dose.
This model requires the Ray Optics Module, CFD Module, and Particle Tracing Module.
Model Definition
The model geometry consists of the annular region between a cylindrical UV lamp and the cylindrical reactor that surrounds it, along with inlet and outlet tubes. All of the domains contain flowing water at room temperature.
Ultraviolet Lamp Model
The surface of the lamp is treated as a diffuse (Lambertian) emitter using the Release from Boundary node. Rays are released at the lamp surface with uniform spatial density, and the initial direction of each ray is sampled according to the cosine law. The total power of the lamp is specified; in this example the power distribution over the lamp surface area is assumed to be uniform, although it is also possible to define a weighting factor.
As rays propagate through the water, away from the lamp surface, the power of each ray is attenuated based on the internal transmittance of the water, which in this example is taken to be 70% per centimeter of propagation. For distilled water, the internal transmittance of germicidal UV radiation is approximately 98%, so the value of 70% used here could represent that the water is less pure.
As rays propagate through the water, the volumetric fluence rate is computed using the dedicated Fluence Rate Calculation node. To obtain an accurate distribution of the fluence rate, it is important to release a sufficiently large number of rays and to use a sufficiently fine mesh in the domain. In this example, 100,000 rays were released in order to balance accuracy with solution time and file size considerations, but in some publications the number of rays could be orders of magnitude greater (Ref. 1).
Rather than releasing rays diffusely at the lamp surface, a possible extension of this model would be to release rays from positions inside the lamp. At the lamp surface, then, the light could be refracted into the water domain according to Snell’s law and the Fresnel equations. Another possible extension is to use the Mirror boundary condition on the exterior surfaces of the reactor, rather than an absorbing Wall boundary condition.
Fluid Flow
The Turbulent Flow, k-ε physics interface was used to solve for the fluid velocity and pressure in the reactor. At the Inlet boundary, the built-in Fully developed flow option was used to define the inflow velocity profile.
The default mesh is controlled by the Turbulent Flow, k-ε interface in this model, which causes boundary layers to grow on surfaces with the Wall boundary condition. Note that a Coarse mesh was used to reduce solution time in this example model, so the results should be verified using a finer mesh. The mesh also controls the resolution of the fluence rate computed by the Geometrical Optics interface, which uses constant shape functions.
Particle Tracing
The Particle Tracing for Fluid Flow interface was used to track particles as they flow past the lamp. The fraction of deactivated or killed bacteria was predicted using a simple exponential decay model with a fixed inactivation constant k (SI unit: m2/J),
where N is the number or concentration of active bacteria at any point along a particle trajectory and N0 is the initial number or concentration. The variable E0 appearing in the time integral is the fluence rate obtained from the Ray Tracing study.
To integrate the fluence rate along each particle trajectory, an Auxiliary Dependent Variable was defined with time derivative equal to E0. At any time, the value of E0 for any particle is found by evaluating the volumetric fluence rate at the particle’s instantaneous position.
Simulations of particles in a fluid are prone to numerical stiffness when the particles are small and the fluid velocity is large. This is because accelerations of small particles in response to drag by the surrounding fluid occur over the time scale τp given by
where
ρp (SI unit: kg/m3) is the density of the particle,
dp (SI unit: m) is the particle diameter, and
μ (SI unit: Pa s) is the dynamic viscosity of the surrounding fluid.
For a micron-sized particle in water, this time scale is on the order of 10-7 seconds. A full inertial treatment of the particle motion would require time steps close to this order of magnitude, whereas the maximum time for particles to cross through the reactor is on the order of 10 seconds. It would therefore be possible to include the effects of particle inertia but the model would be significantly more time-consuming. For a similar reason, the effects of turbulent dispersion have been neglected in this example; each model particle follows a streamline of the mean fluid velocity. In an inertial particle tracing simulation with very small time steps, it would also be possible to get a more realistic solution by enabling turbulent dispersion in the settings for the Drag Force.
Results and Discussion
A Slice plot of the fluence rate is shown in Figure 1. The fluence rate is greatest in the region adjacent the lamp, and significantly smaller at the outermost extents of the reactor volume. Thus, particles passing through the center of the reactor are subjected to more intense UV radiation compared to particles moving along the outer edges.
The fluid velocity is shown in Figure 2 and the particle trajectories in Figure 3. The Color Expression in the trajectory plot shows the fraction of bacteria yet to be inactivated, starting at 1 and with 0 representing complete inactivation.
In Figure 4 a Histogram plot shows the accumulated dose of particles that reach the outlet by the final time. All particles receive a dose of at least 10 mJ/cm2, but for some particles the accumulated dose is substantially higher.
Figure 1: Fluence rate distribution in a cross section of the UV reactor.
Figure 2: Slice plot of the fluid velocity norm, with arrows indicating the velocity direction.
Figure 3: Particle trajectories in the UV reactor. The color expression shows the fraction of surviving bacteria.
Figure 4: Histogram of the accumulated dose along particle trajectories. The plot was filtered to exclude particles that did not reach the outlet by the final time.
References
1. Y.M. Ahmed, M. Jongewaard, M. Li, and E.R. Blatchley III, “Ray Tracing for Fluence Rate Simulations in Ultraviolet Photoreactors,” Environ. Sci. Technol., vol. 52, no. 8, pp. 4738–4745, 2018.
2. D.A. Sozzi and F. Taghipour, “UV Reactor Performance Modeling by Eulerian and Lagrangian Methods,” Environ. Sci. Technol., vol. 40, no. 5, pp. 1609–1615, 2006.
Application Library path: Ray_Optics_Module/Ultraviolet_Sterilization/annular_ultraviolet_reactor_particle
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Optics>Ray Optics>Geometrical Optics (gop).
3
Click Add.
4
In the Select Physics tree, select Fluid Flow>Single-Phase Flow>Turbulent Flow>Turbulent Flow, k-ε (spf).
5
Click Add.
6
In the Select Physics tree, select Fluid Flow>Particle Tracing>Particle Tracing for Fluid Flow (fpt).
7
Click Add.
8
Click  Study.
9
In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Geometrical Optics>Ray Tracing.
10
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Reactor
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, type Reactor in the Label text field.
3
Locate the Size and Shape section. In the Radius text field, type r_reac.
4
In the Height text field, type L_reac.
Inlet Pipe
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, type Inlet Pipe in the Label text field.
3
Locate the Size and Shape section. In the Radius text field, type r_inl.
4
In the Height text field, type L_inl+r_reac.
5
Locate the Position section. In the z text field, type z_inl.
6
Locate the Axis section. From the Axis type list, choose x-axis.
Outlet Pipe
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, type Outlet Pipe in the Label text field.
3
Locate the Size and Shape section. In the Radius text field, type r_inl.
4
In the Height text field, type L_inl+r_reac.
5
Locate the Position section. In the z text field, type L_reac-z_inl.
6
Locate the Axis section. From the Axis type list, choose x-axis.
Union 1 (uni1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
2
Click in the Graphics window and then press Ctrl+A to select all objects.
3
In the Settings window for Union, locate the Union section.
4
Clear the Keep interior boundaries check box.
Lamp
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, type Lamp in the Label text field.
3
Locate the Size and Shape section. In the Radius text field, type r_lamp.
4
In the Height text field, type L_lamp.
5
Locate the Position section. In the z text field, type d_lamp.
Difference 1 (dif1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Difference.
2
3
In the Settings window for Difference, locate the Difference section.
4
Find the Objects to subtract subsection. Click to select the  Activate Selection toggle button.
5
6
In the Geometry toolbar, click  Build All.
Geometrical Optics (gop)
1
In the Model Builder window, under Component 1 (comp1) click Geometrical Optics (gop).
2
In the Settings window for Geometrical Optics, locate the Ray Release and Propagation section.
3
In the Maximum number of secondary rays text field, type 0.
4
Select the Only store accumulated variables in solution check box.
5
Locate the Intensity Computation section. From the Intensity computation list, choose Compute power.
Ray Properties 1
1
In the Model Builder window, under Component 1 (comp1)>Geometrical Optics (gop) click Ray Properties 1.
2
In the Settings window for Ray Properties, locate the Ray Properties section.
3
In the λ0 text field, type 254[nm]. A low-pressure mercury vapor lamp would emit most of its UV radiation at this wavelength (Ref. 1).
Medium Properties 1
1
In the Model Builder window, click Medium Properties 1.
2
In the Settings window for Medium Properties, locate the Medium Properties section.
3
From the n list, choose User defined. In the associated text field, type 1.38. This is approximately the refractive index of water at 254 nm (Ref. 1).
4
From the Optical attenuation model list, choose Internal transmittance, 10 mm sample thickness.
5
From the τi,10 list, choose User defined. In the associated text field, type 0.7. Different values of the internal transmittance of water can be used here, depending on the clarity of the water. For pure water, the internal transmittance of germicidal UV radiation is about 0.98 per centimeter (Ref. 1). The value of 0.7 shown here indicates that the water is less clear.
Release from Boundary 1
1
In the Physics toolbar, click  Boundaries and choose Release from Boundary.
2
3
In the Settings window for Release from Boundary, locate the Initial Position section.
4
From the Initial position list, choose Density.
5
In the N text field, type 100000.
6
Locate the Ray Direction Vector section. From the Ray direction vector list, choose Lambertian.
7
Select the Specify tangential and normal vector components check box.
8
In the Nw text field, type 1.
9
Specify the r vector as
10
From the Sampling from distribution list, choose Random.
These settings will cause rays to be released from the lamp surface with uniform spatial density, with a distribution of initial directions following the cosine law with respect to the surface normal at each release position. To ensure that the surface normal points in the outward direction rather than the inward direction, look for the arrows in the Graphics window.
11
Locate the Total Source Power section. In the Psrc text field, type P.
Fluence Rate Calculation 1
1
In the Physics toolbar, click  Domains and choose Fluence Rate Calculation.
2
Wall 1
1
In the Physics toolbar, click  Boundaries and choose Wall.
2
In the Settings window for Wall, locate the Boundary Selection section.
3
From the Selection list, choose All boundaries.
Turbulent Flow, k-ε (spf)
In the Model Builder window, under Component 1 (comp1) click Turbulent Flow, k-ε (spf).
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
3
In the Settings window for Inlet, locate the Boundary Condition section.
4
From the list, choose Fully developed flow.
5
Locate the Fully Developed Flow section. Click the Flow rate button.
6
In the V0 text field, type 25[gal/min].
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>Water, liquid.
4
Click Add to Component in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Particle Tracing for Fluid Flow (fpt)
1
In the Model Builder window, under Component 1 (comp1) click Particle Tracing for Fluid Flow (fpt).
2
In the Settings window for Particle Tracing for Fluid Flow, locate the Particle Release and Propagation section.
3
From the Formulation list, choose Newtonian, ignore inertial terms.
Drag Force 1
1
In the Physics toolbar, click  Domains and choose Drag Force.
2
3
In the Settings window for Drag Force, locate the Drag Force section.
4
From the u list, choose Velocity field (spf).
Particle Properties 1
1
In the Model Builder window, click Particle Properties 1.
2
In the Settings window for Particle Properties, locate the Particle Properties section.
3
From the ρp list, choose User defined. The default value of the particle density will be used.
Inlet 1
1
In the Physics toolbar, click  Boundaries and choose Inlet.
2
3
In the Settings window for Inlet, locate the Initial Position section.
4
From the Initial position list, choose Density.
5
In the N text field, type 1000.
6
In the ρ text field, type spf.U.
Outlet 1
1
In the Physics toolbar, click  Boundaries and choose Outlet.
2
Particle Counter 1
1
In the Physics toolbar, click  Boundaries and choose Particle Counter.
2
Absorbed Dose
1
In the Physics toolbar, click  Global and choose Auxiliary Dependent Variable.
2
In the Settings window for Auxiliary Dependent Variable, type Absorbed Dose in the Label text field.
3
Locate the Auxiliary Dependent Variable section. In the Field variable name text field, type dose.
4
In the R text field, type gop.frc1.E0.
5
Locate the Units section. Click  Custom Unit.
6
In the Dependent variable quantity table, enter the following settings:
Definitions
Variables 1
1
In the Home toolbar, click  Variables and choose Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Mesh 1
In this example a Free Tetrahedral mesh is used. To get higher resolution of the radial dependence of the fluence rate, a structured mesh could be created, although this may require the addition of mesh control surfaces to the geometry sequence.
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Physics-Controlled Mesh section.
3
From the Element size list, choose Coarse.
4
Click  Build All.
Study 1
Step 1: Ray Tracing
1
In the Model Builder window, under Study 1 click Step 1: Ray Tracing.
2
In the Settings window for Ray Tracing, locate the Study Settings section.
3
From the Time-step specification list, choose Specify maximum path length.
4
In the Lengths text field, type 0 0.2.
5
Locate the Physics and Variables Selection section. In the table, enter the following settings:
6
In the Home toolbar, click  Compute.
Results
Fluence Rate Slice Plot
1
In the Home toolbar, click  Add Plot Group and choose 3D Plot Group.
2
In the Settings window for 3D Plot Group, type Fluence Rate Slice Plot in the Label text field.
Slice 1
1
In the Fluence Rate Slice Plot toolbar, click  Slice.
2
In the Settings window for Slice, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Geometrical Optics>Heating and losses>gop.frc1.E0 - Fluence rate - W/m².
3
Locate the Expression section. In the Unit field, type mW/cm^2.
4
Locate the Plane Data section. From the Plane list, choose zx-planes.
5
In the Planes text field, type 1.
6
Locate the Coloring and Style section. Click  Change Color Table.
7
In the Color Table dialog box, select Thermal>Magma in the tree.
8
9
In the Settings window for Slice, locate the Coloring and Style section.
10
From the Color table transformation list, choose Nonlinear.
11
Set the Color calibration parameter value to -1.5.
12
Click to expand the Quality section. From the Resolution list, choose No refinement.
13
In the Fluence Rate Slice Plot toolbar, click  Plot. Compare the resulting plot to Figure 1.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Physics interfaces in study subsection. In the table, clear the Solve check boxes for Geometrical Optics (gop) and Particle Tracing for Fluid Flow (fpt).
4
Find the Studies subsection. In the Select Study tree, select General Studies>Stationary.
5
Click Add Study in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Stationary
1
In the Settings window for Stationary, click to expand the Values of Dependent Variables section.
2
Find the Values of variables not solved for subsection. From the Settings list, choose User controlled.
3
From the Method list, choose Solution.
4
From the Study list, choose Study 1, Ray Tracing.
5
From the Time (s) list, choose Last.
6
In the Home toolbar, click  Compute.
Results
Slice
1
In the Model Builder window, expand the Velocity (spf) node, then click Slice.
2
In the Settings window for Slice, locate the Plane Data section.
3
From the Plane list, choose zx-planes.
4
In the Planes text field, type 1.
5
Locate the Coloring and Style section. Click  Change Color Table.
6
In the Color Table dialog box, select Linear>Viridis in the tree.
7
Velocity (spf)
In the Model Builder window, click Velocity (spf).
Arrow Volume 1
1
In the Velocity (spf) toolbar, click  Arrow Volume.
2
In the Settings window for Arrow Volume, locate the Arrow Positioning section.
3
Find the x grid points subsection. In the Points text field, type 15.
4
Find the y grid points subsection. In the Points text field, type 3.
5
Find the z grid points subsection. In the Points text field, type 50.
6
Locate the Coloring and Style section. From the Arrow type list, choose Cone.
7
From the Color list, choose Black.
8
In the Velocity (spf) toolbar, click  Plot. Compare the resulting plot to Figure 2.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Physics interfaces in study subsection. In the table, clear the Solve check boxes for Geometrical Optics (gop) and Turbulent Flow, k-ε (spf).
4
Find the Studies subsection. In the Select Study tree, select General Studies>Time Dependent.
5
Click Add Study in the window toolbar.
6
In the Home toolbar, click  Add Study to close the Add Study window.
Study 3
Step 1: Time Dependent
1
In the Settings window for Time Dependent, locate the Study Settings section.
2
In the Output times text field, type range(0,0.1,10).
3
Click to expand the Values of Dependent Variables section. Find the Values of variables not solved for subsection. From the Settings list, choose User controlled.
4
From the Method list, choose Solution.
5
From the Study list, choose Study 2, Stationary.
6
In the Home toolbar, click  Compute.
Results
Particle Trajectories (fpt)
In the Model Builder window, expand the Particle Trajectories (fpt) node.
Color Expression 1
1
In the Model Builder window, expand the Results>Particle Trajectories (fpt)>Particle Trajectories 1 node, then click Color Expression 1.
2
In the Settings window for Color Expression, locate the Expression section.
3
In the Expression text field, type rs.
4
Locate the Coloring and Style section. Click  Change Color Table.
5
In the Color Table dialog box, select Rainbow>Prism in the tree.
6
Particle Trajectories 1
1
In the Model Builder window, click Particle Trajectories 1.
2
In the Settings window for Particle Trajectories, locate the Coloring and Style section.
3
Find the Line style subsection. From the Type list, choose Line.
4
In the Particle Trajectories (fpt) toolbar, click  Plot. Compare the resulting plot to Figure 3.
1D Plot Group 6
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Dataset list, choose Particle 1.
4
From the Time selection list, choose Last.
Histogram 1
1
In the 1D Plot Group 6 toolbar, click  Histogram.
2
In the Settings window for Histogram, locate the Expression section.
3
In the Expression text field, type dose.
4
In the Unit field, type mJ/cm^2.
5
Locate the Bins section. From the Entry method list, choose Limits.
6
In the Limits text field, type range(0,5,120).
7
Locate the Output section. From the Function list, choose Discrete.
8
From the Normalization list, choose Sum of values.
Filter 1
1
In the 1D Plot Group 6 toolbar, click  Filter.
2
In the Settings window for Filter, click Replace Expression in the upper-right corner of the Point Selection section. From the menu, choose Component 1 (comp1)>Particle Tracing for Fluid Flow>Particle Counter 1>fpt.pcnt1.rL - Logical expression for particle inclusion.
3
In the 1D Plot Group 6 toolbar, click  Plot. Compare the resulting plot to Figure 4.