Simulation of Sound Wave Propagation Inside a Spherical Ball Submerged in a Pipeline Wadie R. Chalgham*1, Abdennour C. Seibi1, and Mehdi Mokhtari1 1 University of Louisiana at Lafayette *P.O. Box 40029, Lafayette, LA 70504, [email protected]

Abstract: One of the limitations of pipelines performance and structural integrity assessment is the continuous inspection of possible leaks due to corrosion or other types of failure mechanisms. Efforts to develop new technologies started several decades ago where different inspection techniques were used to enhance pipelines structural integrity. Although available technologies present some advantages, they still have some limitations in practice. This paper presents preliminary numerical simulation results of a research effort aiming at developing a new inspection technique consisting of inserting an autonomous self-recharging spherical ball, equipped with multiple sensors, and moving inside the pipeline. The simulation initially focused on the fluid flow around the spherical ball and noise level propagating inside the pipeline. In addition, COMSOL was used to study the effect of fluid type, ball material, leak location and initial leak noise on the sound pressure level propagation inside the pipe and through the ball. The main goal of this paper is to use the simulation results to calibrate the control system embedded inside the mobile inspection ball. Keywords: Leak Detection, Fluid Flow, Pipeline Inspection, Safety, Environment

1. Introduction Undetected leaks present major financial, environmental, and human threats. Every day, more than 7 billion gallons of clean drinking water are lost due to pipeline failure leading to a loss of $11 billion per year from water leaks only (James, 2011). Moreover, the U.S. Department of Transportation reported that 623 gas and hazardous liquid pipeline incidents happened in 2013 resulting in 10 fatalities, 47 injuries and an estimated cost of $336 million in property damage (Calder, 2014). Thereby, leak detection is essential to mitigate any future incidents from occurring, which require the development of innovative and advanced inspection techniques.

2. Objectives and Scope This paper presents simulation results of the fluid flow around a mobile spherical ball inside a pipeline as well as the leak noise propagation inside the ball. This preliminary study aims at developing a new inspection tool consisting of a mobile ball detecting leaks inside pipelines using acoustic signals. The velocity and pressure profiles will be used to calibrate the control system inside of the ball. The control system relates sound pressure levels to leak detection and accounts for the perturbations caused by the fluid flow around the ball.

3. Use of COMSOL Multiphysics® Software 3.1 Purpose of Use COMSOL Multiphysics was first used to create the geometry of this autonomous ball, which presents a novel design in mobile inspection tools. Then, COMSOL Computational Fluid Dynamics (CFD) module was used to simulate the velocity and pressure propagation around the ball located inside the pipeline for two cases: i) no leak along the pipe and ii) induced leak. The goal of using the CFD module is to find the velocity and pressure profiles around the moving ball. Next, the Acoustics and Vibrations module was used to simulate the leak noise propagation generated when a leak occurs. This module was used to check how leak noise travels inside a pipeline and how it propagates inside the mobile ball. The intended results will be used to provide enough data to build a control system able to detect leak sizes and locations with higher precision. Moreover, COMSOL was used to conduct a sensitivity analysis that compared results for different fluid types, ball materials, leak locations and initial leak noises. The aim of the sensitivity analysis is to provide a better calibration of the control system.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

3.2 Initial and Boundary Conditions

5. Simulation Results

The model under study consists of a fluid flow around a stationary spherical ball placed inside a pipeline as well as the sound propagation generated from an induced leak. In order to accurately understand the collected results, several assumptions were made in the model design: i) single phase laminar flow was considered, ii) constant fluid temperature, iii) constant heat capacity of the fluid (as long as the measurements are at the same depth with the same flow rate then the heat transfer of the fluid will remain the same without random changes to the heat capacity of the fluid), iv) negligible friction-induced heat, and v) the pipe walls are assumed to be sound insulators since this study does not account for external noise.

5.1 Numerical Simulation of Velocity Profiles This simulation used the Computational Fluid Dynamics (CFD) module. Figure 2 shows the velocity profile of the fluid flowing inside the pipe with the leak. The simulation shows how the velocity magnitude of the fluid inside the pipe increases around the ball from the upper and lower sides. This increase is caused by the reduction in area of the flowing water. In fact, at a distance of 12 in. from the inlet, the fluid has a flow surface of π/4*(102 - 62) or 16π in.2 compared to an initial flow surface π/4*102 or 25π in.2. To better visualize the velocity profile around the ball, five vertical sections along the pipe are considered. Figure 3 shows five vertical pipe sections equally spaced starting from the location of 6 in. from the pump to study the change of fluid velocity before, around and after the ball.

4. Numerical Model In order to simulate the fluid flow around the ball, a section of a typical cylindrical pipeline model was developed having a 10 in. diameter as shown in Figure 1. Water was used as the flowing fluid, which was given an initial inlet velocity of 1 mm/s. The ball has a diameter of 6 in. A 0.5 in. diameter leak was induced at an initial location of 25 in. from the pump, which was varied to conduct a sensitivity analysis for a 40 in. pipe section. Vertical Pipe Section (in.)

Figure 4 describes precisely the change in velocity by plotting the magnitude of the velocity along the vertical direction at the five different sections.

Horizontal Pipe Section (in.)

Vertical Pipe Section (in.)

Figure 2. 2D Velocity Distribution of the Fluid around the Ball in Case of a Leak

1

2

3

4

5

Horizontal Pipe Section (in.)

Figure 1. Mobile Ball Flowing inside the Pipeline

Figure 3. Five Vertical Sections of Pipe

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

Figure 4. Average Horizontal Velocity Magnitudes of the Fluid along five Vertical Sections inside the Pipe

Figure 6. Pressure Distribution around the Mobile Ball without a Leak

The results show how the velocity peaks on the upper and lower sides of the spherical ball. As expected, the highest velocity peaks are observed at the Section 3 with magnitudes of 3.1 mm/s, more than 3 times the initial inlet velocity of 1 mm/s.

The plots show that the pressure difference is negligible inside the pipe except around the ball and near the leak location of 25 in. Also, the plots illustrate how the pressure magnitudes decrease significantly around the ball location, from 9 to 15 in., reaching negative magnitudes which could be explained by the high fluid velocity at this location. Moreover, the pressure tends to increase after the ball location reaching a positive magnitude again after 25 in.

5.2 Numerical Simulation of Pressure Profiles This simulation used the Computational Fluid Dynamics (CFD) module. Figure 5 shows the pressure profile of the fluid flowing inside the pipe in the case of a leak. The simulation shows how the pressure reduces significantly when the fluid flows around the spherical ball. In order to understand the effect of the leak on the pressure distribution of the fluid inside a pipeline, another simulation was performed without a leak.

Vertical Pipe Section (in.)

Figure 6 shows a plot of the pressure distribution along a horizontal section below the ball, at an elevation of 2 in.

Leak Location

Horizontal Pipe Section (in.)

The results from Figure 6 confirm Bernoulli’s equation. In fact, the increase in velocity observed in Figure 4 around the ball (from 1 to 3.1 mm/s) resulted in a significant reduction in pressure of more than 5 mPa. 5.3 Numerical Simulation of Sound Pressure Level Propagation inside the Pipeline To simulate a leak noise generated from a leak size of 0.5 in., a point source was created at a location of 25 in. where the leak centerline is located. The spherical ball has a steel layer of 0.5 in. thickness filled with air. An initial leak noise power of 8.5e-11 W was assigned at this location to study sound propagation along the pipe and inside the ball. Figure 6 shows the simulated results of acoustic propagation from the leak location and how the magnitude of the sound pressure level decreases as it gets away from the leak location.

Figure 5. Pressure Distribution in Case of a Leak

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

ball front location of 15 in. from the leak location of 25 in.

Figure 6. Sound Pressure Level in dB in the Presence of a Leak at a location of 25 in.

In addition, the calibration magnitude for this case scenario is ± 1 dB, calculated by subtracting the sound pressure level at the front side of the ball of 92 dB from the rear side having a magnitude of 91 dB. Thus, the expression that will be implemented in the control system for this case scenario is 20 ± 1 dB.

6. Sensitivity Analysis on Leak Noise Propagation

Figure 7. Sound Pressure Level Propagation Inside the Ball Only

In addition, this simulation shows how the sound pressure level increases in the steel layer and decreases inside the spherical ball. In fact, steel has higher density compared to water and air, which allows it to have higher pressure levels. Figure 7 presents a plot of the sound pressure level propagation in decibel (dB) at the centerline of the pipe. The results from Figure 7 show how the sound level increases at the leak location of 25 in. and then decreases and peaks again at the ball location. At the steel outer layer of the ball, the sound pressure level increased by 20 dB from 72 to 92 dB, an increase caused by the higher density of steel. It is important to also note that the sound pressure level at the rear side of the ball has a magnitude of 91 dB, a 1 dB lower than the value at the front side of the ball. These results are significant since they will be used to correctly calibrate the control system, which will be designed in the next phase of this research effort. For instance, for this particular case scenario, an increase of 20 dB detected at the ball location describes that the leak is at a distance of 10 in., calculated by subtracting the

The obtained results so far represent an initial study of one case and needs to be investigated for more scenarios where the fluid type, ball outer layer material, leak noise intensity and leak location are different. This section describes the effect of these parameters. 6.1 Effect of Fluid Type In order to understand the effect of fluid type on the sound propagation inside a pipe, a different fluid medium such as oil was used using the same conditions as for water. Figure 8 shows the effect of the flowing fluid type on the sound pressure level propagation. The results show that the sound pressure levels are lower when using oil as the flowing fluid than when using water. The difference is explained by the difference in density; oil is less dense than water. The maximum recorded difference between the water and oil plots is 1.2 dB. As a result, the calibration for fluid type is not negligible and has a magnitude of ± 1.2 dB.

Figure 8. Fluid Type Effect on Sound Pressure Level Propagation inside a Pipeline

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

6.2 Effect of Ball Material In order to understand the effect of the ball outer layer material on the sound propagation inside a pipe, different simulations were conducted where aluminum and polypropylene were used instead of steel. For comparison purposes, the fluid type was maintained as water. Figure 9 shows the effect of the ball material on the sound pressure level propagation. The results show that the sound pressure level varies with respect to the material type where the sound decreases in aluminum and polypropylene as compared to steel. The difference is explained by the difference in density; aluminum is less dense than steel and polypropylene is less dense than aluminum. The maximum recorded difference between the steel and aluminum/polypropylene is 2.7 and 17.2 dB; respectively. As a result, the calibration for ball material is not negligible and has a peak magnitude of ± 17.2 dB. This calibration parameter depends on the material type used for the ball.

Figure 10. Leak Noise Power Effect on Sound Pressure Level Propagation inside a Pipeline

Figure 10 shows the different sound pressure level magnitudes for various audible initial leak noises. The results show that when the initial leak noise power increases, the sound pressure levels increase as well. The maximum recorded difference in magnitude is 0.48 dB, describing an increase of leak power of 1 e-11 W (from 8.5 to 9.5 e-11 W). As a result, the calibration for initial leak noise is not negligible and has a peak magnitude of ± 0.48 dB per 1e-11 W change in initial leak power. 6.4 Effect of Leak Location The previous simulations used a fixed leak location of 25 in., which is located at a distance of 12 in. (1 ft.) from the center of the ball. In order to understand the effect of leak location on the sound pressure distribution, different leak locations were used.

Figure 9. Ball Material Effect on Sound Pressure Level Propagation inside a Pipeline

6.3 Effect of Initial Leak Noise Power The initial leak noise power used in the previous simulations was assumed to be constant at 8.5 e-11 W. To understand the effect of different initial leak noise powers on the sound pressure level, various simulations with different leak noise magnitudes were compared. A spherical steel ball flowing in water was used.

Figure 11 shows the simulation results of sound pressure level propagation inside the pipe when leak locations are at distances of 13, 25 and 37 in. from the inlet or at distances of 0, 1, and 2 ft. respectively from the front side of the ball. For this sensitivity analysis, the fluid type is water, the ball outer material is steel and the leak initial power is 8.5 e-11 W. The results show that when the leak location is 1 ft. farther from the front side of the ball (i.e., from 13 in. to 25 in.), the recorded sound pressure level magnitudes decrease by a peak value of 23.6 dB.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

9. References

Figure 10. Leak Location Effect on Sound Pressure Level Propagation inside a Pipeline

As a result, the calibration for leak location is very important and has a peak magnitude of ± 23.6 dB per 1 foot change in leak location. This sensitivity analysis showed that the calibration of the control system embedded inside the ball is primarily affected by the leak location since it has the highest calibration magnitude of 23.6 dB when compared to the other calibration parameters.

8. Conclusions The simulation results provided data to calibrate the control system of the ball with respect to the flowing fluid type, ball material, initial leak noise power, and leak location inside a pipeline. The following calibrations were found and represent peak values to have more accurate results:  

 

The fluid type has an effect of ± 1.2 dB on the sound pressure level. The ball outer layer material has an effect of ± 17.2 dB on the sound pressure level received at that outer shell of the ball where the acoustic sensors will be placed. The leak noise has an effect of ± 0.48 dB per 1e-11 W change in power on the sound pressure level. The leak location has an effect of ± 23.6 dB per 1 foot change distance on the sound pressure level.

These calibration values represent the noise level threshold below which any noise level will be detected by the sensors.

1. Chalgham, W. R., Seibi, A. C., and Boukadi, F., Simulation of Leak Noise Propagation and Detection Using COMSOL Multiphysics, ASME Proceedings of the International Mechanical Engineering Congress & Exposition, Phoenix, Arizona, USA (2016) 2. Chalgham, W. R., Seibi, A. C., and Lomas, M., Leak Detection and Self-Healing Pipelines Using Twin Balls Technology, SPE Annual Technical Conference and Exhibition, Dubai, UAE (2016) 3.Ariaratnam, S.T. and Chandrasekaran, M., Development of an Innovative Free-Swimming Device for Detection of Leaks in Oil and Gas Pipelines, Proceedings of the 2010 ASCE Construction Research Congress “Innovation for Reshaping Construction Practices”, ASCE, Banff, Alberta (2010) 4. Bond, A., Deployment of equipment into fluid containers and conduits. United States Patent and Trademark Office, U.S. Patent No. 5,889,703, USA (2005) 5. Bond, A., Mergelas, B. and Jones, C., Pinpointing Leaks in Water Transmission Mains, Proceedings of ASCE Pipelines 2004, San Diego, California, USA (2004) 6. Calder, B., Technologies Mimic the 5 Senses to Monitor Pipelines, Intel Free Press (2014) 7. Chapman, H., Development of a Successful Internal Leak Detection and Pipeline Condition Assessment Technology for Large Diameter Pipes. 6th Annual WIOA NSW Water Industry Engineers & Operators Conference (2012) 8. Chatzigeorgiou D., Khalifa A., Youcef-Toumi K. and Ben- Mansour R.. An in-pipe leak detection sensor: Sensing capabilities and evaluation. ASME/IEEE International Conference on Mechatronic and Embedded Systems and Applications (2011) 9. Dancer, B., Shenkiryk, M. and Day, S., Leak Detection Survey for a Large Diameter Transmission Main: City of Calgary, Proceedings of the 2009 Pipelines Conference, ASCE, San Diego, CA (2009) 10. James, A., The U.S. Wastes 7 Billion Gallons of Drinking Water a Day: Can Innovation Help Solve the Problem? Climate Progress, Think Progress (2011)

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

11. Khalifa A., Chatzigeorgiou D., YoucefToumi K., Khulief Y. and Ben-Mansour R., Quantifying acoustic and pressure sensing for inpipe leak detection, ASME International Mechanical Engineering Congress & Exposition (2010)

DOMAIN EQUATIONS The domain diffusion equation for the soundenergy density w = w(x, t) is given by :

10. Acknowledgements The author is grateful to the University of Louisiana at Lafayette and to the advisor Dr. Seibi, who supported the simulation work and provided continuous support to this research study.

where the local energy flux vector J (SI unit: J/m2/s = W/m2) is defined in the usual way, as:

The equation is simplified as:

11. Appendix GEOMETRY The total diffusion coefficient is Dt = D = λc/3 (SI unit: m2/s), λ is the mean free path (SI unit: m), c is the speed of sound (SI unit: m/s), and ma is the volumetric absorption coefficient (SI unit: 1/m). MESH

Figure A1. Solid Model of Pipe/Ball System

ACOUSTIC DIFFUSION MODEL The acoustic diffusion model is based on the assumption that the volumes studied contain scatters that uniformly scatter the sound field and that the sound field is diffuse (large number of reflections). Using the diffusion of light in a scattering environment as an analogy one can express a diffusion equation for the soundenergy density w = w(x,t) (SI unit: J/m3). The diffusion equation describes the energy flow from high to low energy regions.

Figure A2. Finite Element Model of Pipe/Ball System

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

Pore-Level Bénard–Marangoni Convection in Microgravity Peyman Mohammadmoradi, and Apostolos Kantzas* Chemical and Petroleum Engineering Department, University of Calgary *Corresponding author: 2500 University Dr. NW Calgary, Alberta, Canada T2N 1N4, [email protected]

Abstract: Pore-level displacement of heavy-oil during thermal operations such as SAGD and CSS is a complex multi-scale phenomenon. As gravity drainage is the main depletion mechanism within the intergranular pore space, the surface tension-related phenomena are dominant in intra-granular micro-pores. In a high capillarity pore element, the effect of Marangoni flow is amplified because of extremely reduced viscous and buoyancy-driven flows. Due to the low Capillary and Bond numbers, the fluid interface development is controlled by the geometry and wettability of the solid phase. In this work the COMSOL Multiphysics platform is used to simulate oil mobilization in an oil-wet single-pore geometry as a result of Bénard– Marangoni convection. In particular, the reduction of the trapped oil saturation due to the thermally induced interfacial tension gradient fluxes is monitored. Results indicate that, during the thermal processes, the fluids configuration in high capillarity zones is dictated by the capillary pressure and is deformed by the surface tension gradient across the interface.

are some of the challenges in dynamic microscale simulations. Assuming the capillary force as the only determining factor in microporosity spaces, pore-network, and direct pore morphology-based (App.1) approaches are proposed to eliminate numerical complexities and simulate two- and three-phase immiscible scenarios, Fig.1.

Keywords: Bénard–Marangoni, Pore-scale, Residual oil, Capillary Pressure, Micro-porosity.

Figure 1. Volume fraction distribution in a 3D medium applying pore morphology-based technique (blue and gold represent water and oil)

1. Introduction Recent advances in production technology from unconventional reservoirs and the evolving high-resolution imaging tools are the main parameters attracting much attention to the porelevel experimental and simulation studies. There is a myriad of experimental and numerical simulation studies conducted using pore-level images in a wide range of rock configurations and pore size distributions. In this regard, porescale simulation of thermal processes is not yet fully accomplished. Using CFD-based approaches, by coupling Navier-Stokes and energy equations, one can model immiscible thermal displacement scenarios. However, there are still numerical issues applying Level-Set (LS) or Volume of Fluid (VOF) volume tracking methodologies. Spurious velocities, incapability to conserve mass, interface curvature miscalculations, and high computational times

However, when it comes to thermal processes, surface tension and consequently capillary pressure are variables affecting the interface position and consequently the residual saturations. Here, the effect of thermally induced interfacial tension gradient fluxes, Bénard– Marangoni convection, on the amount of residual oil saturation is investigated in a microgravity environment. Bénard–Marangoni or thermocapillary convection happens when surface tension varies locally as a function of the temperature.

2. Assumptions The following hypotheses are considered in this study:  There are three phases: oil, hot fluid (water or steam) and solid.  The solid phase is strongly oil-wet.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

      

Fluids are compressible and immiscible. Heat transfer happens in both solid and fluid phases. Navier-Stokes and energy equations are solved, simultaneously. The process is non-isothermal. Surface tension and viscosity are temperature-dependent variables. Phase change is not taken into the account. Gravity effect is negligible; microgravity.

3. Use of COMSOL Multiphysics® Software COMSOL Multiphysics is used to solve the governing energy and flow equations. The LS approach is utilized to track the interface development through the pore space.

4. Governing Equations The case study includes two-phase immiscible invasion of hot fluid into an oil-wet geometry to drain oil. The schematic of the case is demonstrated in Fig. 2. The medium is unconsolidated sand and as a result of the nonwetting phase injection, oil remains trapped in high capillarity zones.

Figure 2. Demonstration of the case study (black and gold represent oil and solid phases)

Two studies are conducted simultaneously; heat and two-phase flow. The governing equations for laminar incompressible two-phase flow using LS advection equation are momentum, continuity and LS functions, Eq.1-3. 𝛿𝑢

𝜌 𝛿𝑡 + 𝜌(𝑢. ∇)𝑢 = ∇. [−𝑝 + 𝜌(∇𝑢 + (∇𝑢)𝑇 ] + 𝜌𝑔 + 𝐹𝑠𝑡 + 𝐹

(1)

∇. 𝑢 = 0 𝛿𝜙 𝛿𝑡

(2)

+ 𝑢. ∇𝜙 = 𝛾∇. (𝜀𝑙𝑠 ∇𝜙 − 𝜙(1 −

∇𝜙 𝜙) |∇𝜙|

(3)

The heat transfer in fluids is simulated using time-dependent energy equation, Eq.4-5: 𝜌𝐶𝑝

𝛿𝑇 𝛿𝑡

+ 𝜌𝐶𝑝 𝑢. ∇𝑇 + ∇. 𝑞 = 𝑄 + 𝑞𝑜 + 𝑄𝑝 + 𝑄𝑣𝑑

𝑞 = −𝑘∇𝑇

(4) (5)

𝑄𝑝 and 𝑄𝑣𝑑 refer to the work due to pressure changes and viscous dissipation term, respectively. The heat transfer in solids is only due to conduction and is described by Fourier’s law defining the conductive heat flux proportional to the temperature gradient, Eq.6. 𝛿𝑇

𝜌𝐶𝑝 𝛿𝑡 = ∇. 𝑞 + 𝑄

(6)

5. Case Definition The test case to simulate the thermocapillary effect is defined as introducing hot fluid (steam) to the vicinity of oil trapped in a high capillarity zone. This case would represent a late stage of SAGD when the displacing front has moved past the sand grains. The results are post-processed to determine the effect of temperature imposed phenomena on the residual saturation. 5.1 Solvers Three modules including:  Laminar Two-Phase Flow,  Heat Transfer in Fluids,  Heat Transfer in Solids, are utilized to conduct the numerical simulations including fluid flow in pore space and heat transfer in both solid and pore spaces. 5.2 Boundary conditions  The pressure difference is constant and equal to 30 Pa.  The hot fluid temperature is variable between 300 to 400 K  Oil can be produced via films or bulk flow.  The solid phase is oil-wet. 5.3 Initial conditions  The temperature of the whole system is 300K.  The medium is initially filled by oil.  Density, heat capacity and thermal conductivity of oil are assumed equal to 930 kg/m3, 2000 J/kg.K and 1.1 W/m.K, respectively.  Hot fluid properties are same as water/steam properties.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

5.4 Fluid properties The functionality of viscosity and surface tension of the oil phase are depicted in Fig.3, as two typical curves of bitumen.

constant surface tension, a similar result is achievable using direct quasi-static pore morphology-based technique.

Figure 3. IFT and viscosity of the oil phase

The hot invading phase is considered water/steam and its properties are predefined in COMSOL library.

6. Results Four different scenarios are deployed to distinguish the role of viscosity and surface tension and demonstrate how fluid distribution in high capillarity zones changes during immiscible thermal processes. The simulation scenarios are listed as follows:  Isothermal process  Non-isothermal processes  Variable viscosity  Variable surface tension and viscosity  Variable hot fluid temperature 6.1 Isothermal process The first scenario is to assume a constant temperature during the injection process. A constant pressure difference, 30 Pa, is applied under isothermal and microgravity conditions at 300 K. As the simulation starts, the non-wetting phase gradually invades the medium as a result of pressure difference between two opposite faces. Due to the high viscosity of the oil, the invasion process is very slow. As the interface goes to narrower portions, the interface reaches the equilibrium condition and stops. Fig.4 represents interface development in six steps. Fig.5 shows the final volume fraction distribution at equilibrium condition. Due to the

Figure 4. 2D demonstration of volume fraction distribution during the interface development (left to right, top to bottom) At equilibrium condition, 19.23% of oil remains trapped in the tight portion of the pore, forming a pendular ring, Fig.5.

Figure 5. 2D demonstration of volume fraction distribution at equilibrium condition (Sor=19.23%) In section 6.3, it is demonstrated that temperature gradients can cause further flows and relocate the interface position.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

6.2 Variable viscosity Here, the displacement scenario is a thermal process and the viscosity of the oil is temperature-dependent as shown in Fig.3, while the surface tension remains constant during the process. Hot fluid is injected at 380 K. As it is depicted in Fig.6, the same amount of residual saturation remains immobile. However, due to the thermally imposed lower viscosity of the oil, the equilibrium is reached in a shorter time.

6.4 Variable hot fluid temperature The hot fluid temperature is considered as a variable to perform a sensitivity analysis on the effect of temperate on trapping amount of the wetting phase. The residual oil decreasing trend versus invading phase temperature is plotted in Fig.8.

Figure 6. 2D demonstration of volume fraction distribution at equilibrium condition (Sor=19.61%)

It is concluded that the final saturation distribution in such a single pore is controlled by capillary force and the viscosity difference does not play a determining role. However, in multimacropore systems, the viscosity variations lead to different flow efficiency. 6.3 Variable surface tension and viscosity Here, both viscosity and surface tension are temperature-dependent parameters. The initial and hot fluid temperatures are 300 K and 380 K, respectively. Due to viscosity reduction by heat transfer, in comparison to the isothermal case, the process is faster. More importantly, the temperature distribution induces a surface tension gradient through the oil phase which causes additional streams and consequently, as it is shown in Fig.7, the amount of residual oil is lower.

Figure 8. Residual saturation versus hot fluid temperature

According to the results, it is illustrated that in high capillarity zones the thermally-induced surface tension variations can mobilize the interface and affect the residuals. Therefore, in quasi-static simulations, the surface tension should be adjusted to the invading phase temperature assuming negligible temperature gradient in a microstructure.

7. Conclusions A multidisciplinary study was conducted to investigate the effect of temperature on porelevel capillary dominant displacements. In particular, the effect of thermally induced thermocapillary convection on the amount of residual oil saturation was studied using an oilwet single pore geometry. Results demonstrate that in thermal-based EOR operations, the temperature rise can profoundly change the surface tension in micro-pores and decrease the residual saturation. As the viscous and gravity forces are the main production mechanisms in macro-pores, the surface tension gradient is one of the important phenomena in micro-pores affecting fluid-fluid interface equilibriums.

Figure 7. 2D demonstration of volume fraction distribution at equilibrium condition (Sor=10.2%)

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

8. Nomenclature         

𝑪𝒑 : heat capacity at constant pressure (J/kg/K) g: gravity acceleration (m/s2) k: thermal conductivity (W/m/K) q: heat flux (W/m2) Q: heat source (W/m3) T: temperature field (K) 𝑻𝒔 : steam temperature (K) u: velocity field (m/s) ρ: density (kg/m3)

9. References 1.

2.

3.

4.

5.

6.

7.

Doi, T., & Koster, J. N.. Thermocapillary convection in two immiscible liquid layers with free surface. Physics of Fluids A, Volume 5, 1914-1927, DOI:/10.1063/ 1.858817 (1993). Mohammadmoradi, P. and Kantzas, A., Toward Direct Pore-Scale Modeling of Three-Phase Displacements. Advances in Water Resources, under review (2016). Mohammadmoradi, P. and Kantzas, A., Petrophysical characterization of porous media starting from micro-tomographic images. Advances in Water Resources, Volume 94, 200-216. (2016) Mohammadmoradi, P., and Kantzas, A. Pore-Scale Permeability Calculation Using CFD and DSMC Techniques, Journal of Petroleum Science and Engineering, ISSN 0920-4105, doi:10.1016/j.petrol.2016.07.010(2016) Mohammadmoradi, P., and Kantzas, A., Pore Scale Investigation of Wettability Effect on Waterflood Performance. Society of Petroleum Engineers. doi:10.2118/181309-MS(2016) Mohammadmoradi, P., Behrang, A., and Kantzas, A., Effective Thermal and Electrical Conductivity of Two-Phase Saturated Porous Media. Society of Petroleum Engineers. doi:10.2118/180740MS(2016) Rao, A. R., and Biswal, P. C.. BénardMarangoni convection in two thin immiscible liquid layers with a uniform magnetic field. Acta mechanica, Volume 151, 61-73. (2001)

8.

Scriven, L. E., and C. V. Sternling. "The Marangoni effects. Nature 187, 186 – 188, doi:10.1038/187186a0 (1960) 9. Shevtsova, V.M., Indeikina, A.E., Ryazantsev, Y.S., Thermocapillary motion of a two layered liquid with nonlinear dependence of the surface tension on the temperature, Proceedings of the International Symposium on Hydromechanics and Heat/Mass Transfer in Microgravity, Perm, Russia, Gordon and Breach Science Publ., pp. 157–162. (1992) 10. Wang, C.H., Sen, M., Vasseur, P., Analytical investigation of BénardMarangoni convection heat transfer in a shallow cavity filled with two immiscible fluids, Applied Scientific Research, Volume 48, pp. 5–53. (1991)

10. Acknowledgements The authors gratefully acknowledge the financial support of the FUR program from NSERC, AITF/i-CORE, and the sponsoring companies: Athabasca Oil Corporation, Laricina Energy Ltd., Devon Canada, Foundation CMG, Husky Energy, Brion Energy, Canadian Natural, Maersk Oil, Suncor Energy, and Schulich School of Engineering (University of Calgary).

11. Appendix The pore morphology-based simulation of capillary dominant scenarios could be summarized in four main steps as follows:  Thresholding gray-scale raw images.  Extraction of pore size distribution using expanding sphere concept.  Applying geometrical rules and perform drainage and imbibition processes considering by-passing and snap-off as trapping mechanisms. The accessibility and geometry of each pore determines its saturation state.  Post-processing partially saturated media by extracting volume fraction propagation maps.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

1

Glass Transition of ABS in 3D Printing Miftahur Rahman*1, N. R. Schott2 and Lakshmi Kanta Sadhu3 1 Professor, North South University, Dhaka, Bangladesh. 2 Professor Emeritus, Department of Plastics Engineering, UMASS Lowell, USA. 3 IRays Teknology Ltd., Dhaka, Bangladesh. *Corresponding author: [email protected], [email protected] 

Abstract— In a commercial 3D printer head, plastic ribbon passes through a hot nozzle of an extruder to dispense liquid plastic droplets to construct the model. In this paper a 2D axisymmetric model of a 3D head is considered to study the secondary transition change from below the glass temperature to above the glass-transition temperature of ABS (Acrylonitrile Butadiene Styrene) using COMSOL Multiphysics Software. To achieve accurate results, the melt flow fields in combination with the heat transfer and change in tensile modulus are considered. The model includes the secondary transition, both in terms of latent heat and in terms of other thermodynamic and physical variables. The simulation is based on the assumption that there is no volume change during solidification of ABS. It is also assumed that the velocity of the melted ABS is constant and uniform throughout the modeling domain. The transition from the glassy to the molten state of ABS is modeled using the heat capacity model. A narrow secondary transition is observed during the glass-transition temperature. Index Terms—ABS (Acrylonitrile Butadiene Styrene), Computer Aided Design (CAD), Computer-Aided Manufacturing (CAM), Computer Numeric Control (CNC), Fused Deposition Modeling (FDM), Polylactic Acid (PLA), Selective Laser Sintering (SLS), Stereolithography (SLA), Thermoplastic Polyurethane (TPU).

N

I. INTRODUCTION

OT all 3D printers make use of the extruder head for dispensing melted plastic. There are several other 3D printing techniques that are available, differing mainly in the way the layers are built to create the final object. Selective laser sintering (SLS) and fused deposition modeling (FDM) are the most common technologies that make use of melting or softening of material to construct the layers. In a stereolithography (SLA) technique a photo-reactive resin is cured one layer at a time with a UV laser. Fused deposition modeling (FDM), a method of rapid prototyping, makes use of a nozzle ejecting molten thermoplastic and is used to model the object with a controlled movable vertical platform. The FDM technology works using a thermoplastic filament or metal wire which is unwound from a coil and supplying material to an extrusion nozzle which can turn the flow on and

off. The nozzle is heated to melt the material and can be moved in both the horizontal and vertical directions by a numerically controlled mechanism similar to conventional Computer Numeric Control (CNC) machine directly controlled by a computer-aided manufacturing (CAM) software package. The model is produced by extruding melted material to form layers as the material hardens immediately after extrusion from the nozzle. This technology makes use of two thermoplastic materials namely: ABS and PLA (Polylactic acid). Other thermoplastic materials that can also be extruded are high-impact polystyrene (HIPS), thermoplastic polyurethane (TPU), aliphatic polyamides (nylon) and recently also PEEK. Paste-like materials such as ceramics and chocolate can also be extruded using the fused filament process and a paste extruder. Hard plastics like ABS, PMMA (Polymethyl methacrylate) have end use applications well below their glass transition temperatures that is in their glassy state. Their Tg values are well above room temperature, both at around 105 °C (221 °F) [1]. Rubber elastomers like polyisoprene and polyisobutylene are used above their Tg, that is, in the rubbery state, where they are soft and flexible [2]. In polymers the glass transition temperature, Tg, is often expressed as the temperature at which the Gibbs free energy is such that the activation energy allows the cooperative movement of molecular chain segments to slide past each other when a force is applied. When the glass temperature has been reached, the stiffness stays the same until the temperature exceeds Tg, and the material turns rubbery. This region is called the rubber plateau.

II. 3D PRINTER EXTRUDER HEAD In a conventional extruder a high-volume thermoplastic plastic is melted and forced out through a nozzle to cool into various shapes and profiles such as pipes, tubes, ribbons, films, coatings, wire insulation, or sheet depending on the various kinds of die as shown in Fig. 1. During the extrusion process, plastic material (pellets, granules, flakes or powders) are fed from a hopper into the heated barrel of the extruder. The plastic pellets or powders are gradually plasticated by the heated barrel and pushed through the nozzle by the turning screws and by heaters. The molten polymer is then forced into

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

2 a die, which shapes the polymer into a desired shape that hardens during cooling. Fig. 1 [3] shows a schematic diagram of an extruder with a feeding hopper through which the plastic pellets or resins enter into a plasticating screw in which pellets get turned into a viscous liquid which is dragged forward by the rotating screw. The melt reaches the end of the cylinder where it is forced through a screen pack and the die to form the desired shape.

Fig. 1 Schematic diagram of a conventional extruder [3].

The technology of a conventional extruder has been adapted to a commercial 3D printer as shown in Fig. 2 [4]. In Fig. 2, instead of a hopper, a 3D printer uses a continuous plastic filament which is fed from a large spool or reel and pushed through a set of gears to the heated extruder barrel. Molten plastic is forced through the nozzle and is deposited on the working table. The 3D extruder head moves, under computer control, to define the printed shape. Similar to a CNC machine, the 3D printer head is controlled by a computer along with electronic controllers, stepper or servo motors to move and manipulate the mechanical axes of the 3D printer. Usually the head moves in layers, moving in two dimensions to deposit one horizontal plane at a time, before moving slightly upwards to begin a new slice.

molten ABS flows continuously through a narrow nozzle. Also, to optimize the casting process in terms of casting rate and cooling, the thermal and fluid dynamic aspects of the process are also considered. To obtain accurate results, the melt flow fields in combination with the heat transfer and glass (secondary) transition are considered. The model includes the transition from glassy solid to rubbery melt, both in terms of sensible heat along with other physical properties. The model assumes a steady state and continuous process. The heat transport is described by the equation: ρC P u.T  .(kT)  Q (1) where k, Cp and Q denote thermal conductivity, specific heat, and heating power per unit volume (heat source term), respectively. The velocity, u, is the fixed pulling speed of the plastic ribbon in both the rubbery liquid and the glassy states. As the melt cools down in the air, it becomes a glassy rigid liquid. During the glass transition, a significant amount of heat is released. The total amount of heat released per unit mass during the glass transition is given by the change in enthalpy, ΔH. The heat capacity, Cp, changes considerably during the glass transition. The difference in specific heat before and after the transition can be approximated by ΔH ΔCP  . (2) T As opposed to neat polymer, a polymer with fillers or reinforcements generally undergoes a broad temperature transition zone, over several degrees Kelvin in which a mixture of both solid and molten polymer co-exist in a “mushy” zone. To account for the sensible heat related to the glass transition, the apparent heat capacity method is used through the Heat Transfer with a Phase Change domain condition. The objective of the analysis is to make ΔT, the half-width of the transition interval small, such that the solidification front location is well defined. Furthermore, models considered the laminar flow by describing the fluid velocity, u, and the pressure, p, according to the equations u  ρu.u  t    2μ  .  pI  μ(u  u T    κ .u I   F (3)  3    ρ  .( u )  0 (4) t where ρ is the density (in this case constant), μ is the viscosity, and κ is the dilatational viscosity (here assumed to be zero). Here, the role of the source term, F, is to dampen the velocity at the modulus-change interface so that it becomes that of the rigid phase after the transition from the rubbery. The source term has been taken from equation [5] ρ

Fig. 2 Schematic diagram showing how plastic ribbon is fed from a large coil through a set of moving gears to a heated barrel of an extruder of a 3D printer head [4].

III. NUMERICAL MODELING A 2D axisymmetric heat transfer model is considered in the numerical modeling of the 3D printer head. In this simulation

1 - α 2 A

(5) mush  u - u cast  α3  ε where α can be seen as the volume fraction of the liquid phase; Amush and ε represent arbitrary constants (Amush should F

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

3 be large and ε small to produce a proper damping); and ucast is the velocity of the extruded plastic.

IV. MATERIAL PROPERTIES OF COMMONLY USED THERMOPLASTICS The properties of plastics and processes are influenced by the thermal characteristics such as melt temperature (Tm), glass transition temperature (Tg), thermal conductivity (k), and heat capacity (CP). Table I [6] shows thermal properties of some of the commonly used thermoplastics along with steel. The Tm occurs at a relatively sharp point for crystalline materials. Amorphous materials basically do not have a Tm; they are in a melt state but the viscosity is too high for flow. As the glassy rigid fluid is heated its viscosity decreases exponentially to enable flow. In reality there is no single flow temperature, but rather a range, which is often taken as the peak of a differential scanning calorimetry (DSC) curve. The glass-transition temperature (Tg) is the point below which plastics behaves as glass – it is very strong and rigid, but brittle. Above this temperature it is neither as strong or rigid as glass, but neither is it brittle. The amorphous TPs have a more defined Tg. TABLE I Thermal Properties of TPs [4] Plastics

Density

PP (C) HDPE (C) PTFE (C) PA (C) PET (C) ABS (A) PS (A) PMMA (A) PC (A) PVC (A) Steel

 (kg / m3 )

Melt Temp., Tm (K)

Thermal Conductivity k (W /(m.K ))

Heat Capacity C P ( J /( kg.K ))

900 960 2200 1130 1350 1050 1050 1200

441 407 603 533 523 378 373 368

0.1-0.22 0.45-0.52 0.25 0.24-0.28 0.15-0.4 0.17 0.1-0.13 0.17-0.19

1700-1900 1900 1000 1700 1200-1350 1080-1400 1200 1400-1500

1200 1350 7900

539 472 3023

0.19-0.22 0.12-0.25 43

1200 1000-1500 466

V. COMSOL SIMULATION PARAMETERS The simulation of a 3D extruder head is a highly nonlinear problem. As a result an iterative approach is taken for finding the solution. The location of the transition between the molten rubbery and the glassy state is a strong function of the casting velocity and the cooling rate in the air. A fine mesh is considered across the rubbery to glassy front to resolve the change in material properties. However, it is difficult to know where this solidification front is located. By starting with a gradual transition between rubbery and glassy, it is possible to find a solution even on a relatively coarse mesh. This solution can be used as the starting point for the next step in the solution procedure, which uses a sharper transition from the rubbery liquid to the glassy solid. This is done using the continuation method. Given a monotonic list of values to solve for, the continuation method uses the solution of the last case as the starting condition for the next. Once a solution is found

for the smallest desired ΔT, the adaptive mesh refinement algorithm is used to refine the mesh to put more elements around the transition region. This finer mesh is then used to find a solution with an even sharper transition. This can be repeated as needed to get better and better resolution of the location of the solidification front. In this paper, the parameter ΔT is first ramped down from 300 K to 100 K, then the adaptive mesh refinement is used such that a finer mesh is used around the solidification front. The resultant solution and mesh are then used as starting points for a second study, where the parameter ΔT is further ramped down from 20 K to 17 K. The double-dogleg solver is used to find the solution to this highly nonlinear problem. Although it takes more time, this solver converges better in cases when material properties vary strongly with respect to the solution. ABS parameters that are used in the COMSOL simulation are given in Table-II, and Table-III. TABLE-II THERMAL PARAMETERS OF ABS USED IN THE COMSOL SIMULATION. Description Processing temperature Temperature transition zone half width Heat of glass transition Heat capacity at constant pressure, glassy state Heat capacity at constant pressure, rubbery state Ambient temperature Melt inlet temperature Casting speed Heat transfer coefficient, ABS Heat transfer coefficient, air Surface emissivity

Data Used 378[K] 75[K] 207[kJ/kg] 1200 [J/(kg.K] 1797.6 300 [K] 378 [K] 1.0 [mm/s] 2000 [W/(m^2.K] 10 [W/(m^2.K)] 0.95

Table-III Thermal parameters of ABS used in the COMSOL simulation. Glassy ABS

Value

Unit

Thermal conductivity

0.3

W/(m.K)

Density

1050

kg/m^3

Ratio of specific heats

1

1

Rubbery ABS

Value

Unit

Thermal conductivity

0.2

W/(m.K)

Density

1050

kg/m^3

Ratio of specific heats

1

1

VI. RESULTS AND ANALYSIS The plot in Fig. 3 shows the streamlines close to the nozzle resulting in a vortex. This eddy flow could create problems with a non-uniform surface quality in a real process. Process engineers can thus use the model to avoid these problems and find an optimal nozzle design. To help determine how to

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

4 optimize process cooling, Fig. 4 plots the conductive heat flux. It shows that the conductive heat flux is very large in the nozzle area. This is a consequence of the heat released during the secondary phase transition, which is cooled by air. An interesting phenomenon of the process is the occurrence of the peak of the conductive heat flux appearing in the center of the flow at the transition zone. Furthermore, by plotting the conductive heat flux at the outer boundary for the process as in the lower plot in Fig. 4, one can see that a majority of the cooling process occurs just outside the nozzle. More interestingly, the heat flux varies along the nozzle wall length. This information can help in optimizing the cooling of the model (that is, the cooling rate and choice of cooling method). One can solve the model using a built-in adaptive meshing technique. This is necessary because the transition zone—that is, the region where the glass transition change occurs— requires a fine mesh size. Fig. 5 is a 3D plot of the velocity magnitude obtained by a revolution of the 2D axisymmetric data set. In Fig. 5, molten ABS passes through the narrow hot nozzle and solidifies as an extruded ABS. Fig. 6 shows the surface temperature distribution with the heat flux direction as indicated by arrows. Fig. 7 shown the fraction of liquid rubbery state along the centerline for all values of ΔT. For smaller values of ΔT, the transition is steeper. Fig. 8 shows the sensible heat as a function of arc length for various values of ΔT. As ΔT gets smaller, in particular ΔT = 10 K, is not entirely smooth. As a result, to model with reduced ΔT increments one requires to increase the mesh resolution.

VII. CONCLUSION Modeling and simulation of 3D printer head using COMSOL software tools may predict the die design that gives the best properties for the 3D model. Die design may control the heat flux to give the fastest cure rate to make the model strong in least amount of time. The effects of fillers and their influence on cure rate and end use properties can be predicted.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

5

REFERENCES [1] [2] [3] [4] [5] [6]

PVC Handbook. Hanser Verlag. ISBN 1-56990-379-4. ISO 11357-2: Plastics – Differential scanning calorimetry (DSC) – Part 2: Determination of glass transition temperature (1999). http://www.pitfallsinmolding.com/extrusion1.html https://en.wikipedia.org/wiki/Fused_filament_fabrication V.R. Voller and C. Prakash, “A fixed grid numerical modeling methodology for convection—diffusion mushy region phase-change problems,” Int.J.Heat Mass Transfer, vol. 30, pp. 1709–1719, 1987. Donald V. Rosato, David P. DiMattia and Dominick V. Rosato, “Designing with Plastics & Composites, A Handbook”, Van Nostrand Reinhold, New York, 1991.

Fig. 7 The fraction of liquid rubbery state along the centerline for all values of T. For smaller values of T, the glass transition is steeper.

Fig. 8 As T gets smaller, in particular T  10 K, the curve is not entirely smooth. As a result, to model with reduced T increments one must to increase the mesh resolution.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

Study Thermohydraulic a Fixed Bed for Core of Nuclear Reactor J. Almachi*1, J. Montenegro1, A. Portilla1 1 Escuela Politécnica Nacional, Quito, Pichincha, Ecuador *[email protected], [email protected]

Abstract: The fixed beds have the advantage of large heat transfer area, for they are used in some designs of innovative nuclear reactors as the reactor FBNR. Inside of study of fixed beds is important to define the following parameters: flow minimum and velocity profile of cooling fluid, temperature profile of cooling fluid and the fuel elements. In order to determine the parameters before described the Software COMSOL Multiphysics® was used, the geometry was established in 3D for core of nuclear reactor, and for termohydraulic study modules used were heat transfer, conjugated heat transfer and flow. In this study was considered to nuclear fuel spheres as heat sources, the cooling fluid used in the simulation was water and the type of study was stationary. The study of termohydraulic parameters is important to establish safe operating conditions to verify that the thermal limitations of the materials involved in the reactor core are not exceeded. Keywords: Fixed Bed, Heat Transfer, FBNR, Nuclear Reactor.

1.1 Description of Nuclear Reactor FBNR The fixed bed nuclear reactor, FBNR bases its operation in pressurized water reactors (PWR). Has a small size (2m in diameter and 6 m high) compared to PWRs, generates an output of 70 MW, is modular in design, it does not require refueling site, since it can operate long periods with fuel loaded in the factory building and has an inherent safety system. All these features allow the reactor FBNR be part of the IAEA in the INPRO program, because it meets the requirements set for nuclear reactors of fourth generation, security and economic profitability [6][7][8]. The FBNR reactor is a research proposal of PhD. Farhang Sefidvash, fundamentals and technology of this reactor is public domain and is currently in the study phase. The graphic representation of the FBNR can be seen in the Figure 1 [9].

1. Introduction Faced with the growing demand of energy of world population, nuclear energy captures attention because is secure, competitive and does not produce greenhouse gases, there more covered the 15 % of demand the electricity worldwide in 2015[1][2]. Since 2000 until today so as to potentiate the inclusion of nuclear energy in the energy matrix, are being developed new concepts of nuclear reactors, known as nuclear reactors of fourth generation, through which seeks to improve the efficiency of heat transfer in comparison with reactors of previous generations, therefore the importance of thermohydraulics study of fixed beds has like advantage their larger surface area of heat transfer in the reactor core [3][4][5]. Under this concept are being developed the reactor FBNR (Fixed Bed Nuclear Reactor).

Figure 1. Schematic Design of FBNR

In the design of FBNR are identified 3 modules, with a specific and independent function which are: the core, the steam generator and the fuel chamber.

1.2 Core of Nuclear Reactor FBNR According to the design proposed by Sefidvash F., the reactor core is constituted by a cylinder 1.71 m in diameter and 2 m high, the fuel spheres are inside forming a fixed bed, through which pressurized water circulates, the water acts as cooling fluid, which also ensures that heat removing the materials from which the fuel elements do not melt, the refrigerant ascends through the core and is recirculated to operate in a closed circuit [10]. The fuel elements are based on High Temperature Gas Cooled Reactor (HTGR) technology, in which the fuel elements are spherically shaped. Fuel spheres are made with CERMET, this material is characterized by the high temperature resistance. The spherical CERMET fuel consists of coated UO2 kernels embedded in a zirconium, or zirconium hydride, matrix which is then overcoated with a protective outer fuel-free layer [11]. The Table 1. shows the values of the physical and thermal properties CERMET material used in the present study [12]. Table 1. CERMET properties Properties Thermal Conductivity Specific Heat Density

Temperature 300°C 28.2 W/m°K 431 J/kg.°K 6940 kg/m3

 Select module conjugate heat transfer within this module select laminar flow section  Select stationary study 2.1.2 Model Creation Model creation is done in the geometry module by performing the following steps:  Definition of material, choose the material (water) from the library of Software for the cooling fluid and create a new material for the spheres (CERMET) entering the thermophysical properties described in Table 1.  In the CAD module, create the cylindrical container of the bed and create the spheres placing them in their respective coordinates according to the required structure BCC and FCC 2.1.3 Multiphysics study  Select one of the circular faces of the cylinder and in the section laminar flow create a condition of entry that corresponds to the workflow and temperature input, output condition corresponding to the output  Select the spheres and in the section of heat transfer in solid, creating a heat source, which should indicate the overall transfer rate which corresponds to the value of the total power 2.1.4

Meshing

 In the module meshing create a normal mesh  In the study section place the value of 200 iterations and an error range of 20 2.1.5 Results

2. Use of Software

COMSOL

Multiphysics®

In this section is described as COMSOL Software modules that were used, the required steps to create the model for core of the reactor FBNR, the initial conditions and established considerations for simulation 2.1 Modules of COMSOL Software used in the simulation 2.1.1 Entrance to simulation interface  Select the wizard modeling  Select 3D Model

 In the results module select the data set  Create lines and planes in the CAD model to later extract data  Select show graphics  Import data from lines and planes created in TXT format. 2.2 Considerations for the simulation The considerations for the simulation were:  Steady State Study  Types of fixed bed structures: BCC and FCC

temperature of the fluid is important because this is then studied for the design of the heat exchanger. The energy that wins the cooling fluid of the spheres is represented in the below Figure 4.

2000 1500

Temperature [°K]

500 0.5

1.0

1.5

2.0

0.0

0.5

1.0

1.5

2.0

Length [m]

Temperature distribution, Flow 0.8[m/s]

Temperature distribution, Flow 1[m/s] 2500

Length [m]

2000 1500

Temperature [°K]

1500

500

1000

T min=563.01 °K, T max=2100.97 °K

1000

2000

T min=563.01 °K, T max=2113.18 °K

500

Temperature [°K]

2500

0.0

Table 2. Initial conditions for simulation of core of FBNR Initial conditions Value Inlet temperature 290 °C Heat generation 134.4 W per sphere Inlet Pressure 16 MPa

T min=563.01 °K, T max=2127.56°K

1000

2500 2000 1500

Temperature [°K]

T min=563.01 °K, T max=2173.89 °K

500

The initial conditions established by Sefidvash F. are shown in Table 2 [13].

Temperature distribution, Flow 0.6[m/s] 2500

Temperature distribution, Flow 0.4[m/s]

2.2 Initial Conditions

1000

 Thermo physical properties of the materials involved remain constant  There is no phase change  There uniform heat generation in the particles that make up the fixed bed

0.0

0.5

3. Experimental Results

1.0

1.5

2.0

0.0

0.5

1.0

Length[m]

1.5

2.0

Length [m]

Figure 3. Temperature of spheres in BCC configuration

3.1 BCC Configuration The nuclear fuel in the form of spheres can support a maximum temperature of 3113 ° K. This fundamental parameter in the simulation because the critical temperature points are located in the central part of the fixed bed, as shown in Figure 2.

700

Temperature Profile Fluent, BCC Arrangement Fluid Speed 0.6 m/s

0.8 m/s

1 m/s

680

0.4 m/s

640 620

T=631.55 °K T=614.41 °K T=604.09 °K

560

580

600

Temperature [°K]

660

T=665.74 °K

0.0

0.5

1.0

1.5

2.0

Length [m]

Figure 2. BCC configuration

3.1.1 Temperature of spheres The critical area analysis is in the middle of the bed is therefore important to know what happens to each area along the bed length, the Figure 3 shows the status of the spheres at different flow. 3.1.2 Temperature of fluid Spheres transferred heat energy to the fluid, this enters with a temperature of 563 ° K, the outlet

Figure 4. Temperature of fluid in BCC configuration

3.1.3 Pressure BCC arrangement has greater porosity because in this type of bed more open spaces, this parameter helps determine the pressure drop that the reactor core may have, the pressure curves are shown in the following Figure 5. 3.2 FCC Configuration FCC conditions mentioned in accordance BCC are valid in the FCC arrangement. The following

Figures 6, 7, 8, 9 helps visualize is happening in the reactor core. 800

Temperature Profile Fluent, FCC Settlement

0.6 m/s

0.8 m/s

1 m/s

700

750

T=750.94 °K

T=688.04 °K

T=656.48 °K

650

Temperature [°K]

T=637.44 °K

600

15950000

Fluid Velocity 0.4 m/s 0.6 m/s 0.8 m/s 1 m/s

0.0

0.5

1.0

1.5

2.0

Length [m]

Figure 8. Temperature of fluid in FCC configuration 0.5

1.0

1.5

2.0

Different Pressures profile Flows, Arrangement FCC

16000000

length [m]

Fluid Velocity 0.4 m/s 0.6 m/s 0.8 m/s 1 m/s

14500000

Pressure [Pa]

Figure 5. Pressure in BCC configuration

15500000

0.0

15000000

Pressure [Pa]

15850000 15750000

Fluid Speed 0.4 m/s

Different Pressures Profile Flows, Arrangement BCC

0.0

0.5

1.0

1.5

2.0

Length [m]

Figure 6. FCC configuration

2500

Temperature distribution, Flow 0.6[m/s]

2000 1500

Temperature [°K]

1500

500

1000

0.5

1.0

1.5

2.0

0.0

0.5

1.0

1.5

Length [m]

Temperature distribution, Flow 0.8[m/s]

Temperature distribution, Flow 1[m/s] 2500

Length [m]

1500

2000

T min=563.01 °K, T max=2140,3 °K

500

500

1000

1500

Temperature [°K]

2000

T min=563.01 °K, T max=2164 °K

2.0

1000

2500

0.0

Temperature [°K]

4. Discussion

T min=563.01 °K, T max=2203.3 °K

1000

2000

T min=563.01 °K, T max=2281.9 °K

500

Temperature [°K]

2500

Temperature distribution, Flow 0.4[m/s]

Figure 9. Pressure in FCC configuration

0.0

0.5

1.0 Length [m]

1.5

2.0

0.0

0.5

1.0 Length [m]

Figure 7. Temperature of spheres in FCC configuration

1.5

2.0

Figure 10 shows that there are points where the fluid stagnates and as a result those points higher concentration generate heat this happens for both FCC and BCC arrangement. Comparing Figure 10 with Figure 2 and 3 must points having lower speed higher temperature acquired, this is because the convection coefficient at these points is small. This phenomenon helps determine until temperature can reach areas in all flows presented the interior temperature of the area does not exceed the maximum working. As shown in Figures 5 and 9 the pressure drop is very small and increases with the fluid velocity. The ideal design parameter is 0.6 m / s because

pre hydraulic designs consent to be used a pump with features of 16000000 Pa and 0.6 m / s. With these parameters fluid temperatures and areas under safe working parameters, shown in Figures 3 , 4, 7, 8 decreased fluid velocity will cause the spheres acquire larger temperatures implying increased monitoring reactor core.

As seen in Figure 11 the arrangement BCC has lower pressure drop and higher heat concentration which is beneficial for the cooling fluid. Figures 3, 4, 7, 8 shows that the arrangement BCC is more stable and easier to create a traceability control each of the areas along the bed, it benefits the periodic maintenance and analysis and material control nuclear which is operational. Create a BCC arrangement is less costly than an FCC arrangement.

5. Conclusions

Figure 10. Graphic of velocity

 BCC arrangement is more stable than the FCC arrangement has greater benefits to the cooling fluid, Traceability control of each of the areas of nuclear material increases compared the FCC.  BCC configuration having higher porosity has a lower pressure drop compared to the FCC structure, this causes reduction of energy loss in the operation of recirculation pumps predesigned.

Figure 11. Pressure drop and Temperature for BCC and FCC configuration

Pressure parameters are initial conditions of the simulation for determining the temperature profile both fluid and spheres.  In all cases both under FCC and BCC nuclear material does not exceed the limit operating temperature as long as the fluid velocity is between 0.4 and 1 m / s, if the fluid velocity decreases can raise the temperature of the spheres and creating an unstable and dangerous system.  This type of analysis by simulation with COMSOL allows evaluate, correct and improve designs without creating highly costly prototypes, simulation time in each of the cases had an average of 5 hours, the interactive interface that COMSOL presents excellent since you can see the convergence of the results in real time and make the necessary corrections to improve results  The data obtained in this paper will used to set the required parameters in the study of neutronic in nuclear material present in the Fixed Bed Nuclear Reactor. Such studies can be performed with other modules of COMSOL multiphysic offers, remaining well as future work to study this type of physical phenomenon.  The stability of treatment and data acquisition COMSOL is robust, it allows linking several physical phenomena in the same time interval thus reducing costs and post processing simulation data. Whether the graphic presentation and interactive design makes this software unique when simulate geometries with different metaphysical robust interactions.

6. References 1. IEA, Worl Energy Outlook 2015, 4-6 (2015) 2. Perera J, Impulso a la innovación. Boletín IAEA 2004, 45-46 (2004) 3. Diamonds P, Abundancia el futuro es mejor de lo que piensas, 123-126. Antoni Bosh Editor S.A (2013) 4. Nuclear Energy Agency. The Nuclear Energy Today. 54. (2012)

5. Achenbach E, Heat and Flow Characteristics of packed beds, Experimental thermal and fluid science, 17-25 (1995) 6. IAEA, Status report 72- Fixed Bed Nuclear Reactor (FBNR), 1-3 (2011) 7. Gonzáles, L. Reactores Nucleares de Cuarta Generación,. 5, (2010) 8. Sefidvash F, Thermal hydraulics of the fixed bed nuclear reactor concept, Kerntechnik, 12-18, (20069 9. Sefidvash F, Non-Proliferation Resistance and Physical Protection of FBNR Nuclear Reactor, IAEA Technical Meeting, 5-7 (2010) 10. Sefidvash F. Conceptual design of the fixed bed nuclear reactor (FBNR) Concept, IAEA Report, 123-125, (2005) 11. Vatulin A, Stetsky Y et all, Developed of Cermet fuel, 1-18, (2005) 12. Senor D, Painter C, et all, A new Innovative Spherical Cermet Nuclear fuel Element to Achieve an Ultra Long Core Life for use in Grid Appropiate LWRs, 3.14-3.15 (2007) 13. Sefidvash F, Water Cooled FBNR Nuclear Reactor, IAEA-CN- 164, 1-10, (2008)

7. Acknowledgements The authors are grateful to Prof. Angel Portilla Head of Heat Transfer Laboratory of the National Polytechnic School, Ing. William Venegas Teacher of the Finite Element of the National Polytechnic School and the PhD. Farahng Sefidvash “Prometeo” Professor of the National Polytechnic School, for providing the necessary facilities for the preparation of the paper.

Modeling Research Reactor Fuel Plate Hotspots with COMSOL’s Thin Layer and Thermal Contact Features Michael J. Richards*1, Arthur E. Ruggles1, James D. Freels2 1 The University of Tennessee, 2Oak Ridge National Laboratory *Corresponding author: 5029 Jade Pasture Ln, Knoxville, TN 37918, [email protected]

Abstract: One challenge with converting nuclear research reactors to low enrichment fuels is the thermal-fluid performance of those fuels. Local areas of high temperature, known as hotspots, limit reactor performance and thus require accurate modeling. A simplified fuel plate model is developed to compare traditional FEA techniques with COMSOL’s thin layer (TL) and thermal contact (TC) features for modeling regions of low thermal conductivity and high energy generation that produce hotspots. Temperatures, conductive heat fluxes, and energy balances are reported. TL and TC offer similar performance with < 0.3% error in energy balance for all cases. Both experience oscillations in boundary, normal heat flux at discontinuities in conductivity that decrease with grid refinement and increased element order. TC offers options that make it more appropriate for hotspot models. Keywords: hotspot, thermal contact, thin layer

1. Introduction Nuclear research reactors offer unique capabilities that provide opportunities for material irradiation and neutron scattering not available in power reactors. The use of high enriched uranium (HEU) in some research reactors creates a nuclear material proliferation vulnerability. To help prevent proliferation of nuclear weapons while allowing the continued use of these reactors, efforts are underway to convert from HEU to low enriched uranium (LEU) fuel. Among the challenges with this conversion is the thermal-fluid performance of the fuel. A typical fuel plate consists of fuel sandwiched between aluminum cladding. During the manufacturing process two defects of interest may occur. First, the fuel may not be distributed evenly and local areas of higher enrichment, known as fuel segregations, may result. These generate more energy than the surrounding fuel.

Second, the fuel plate may contain non-bonds, regions where the cladding does not bond to the fuel. These reduce the thermal conductivity across the fuel-cladding interface which creates a hotspot on the opposite side of the plate. The worst case scenario is when these two defects are coincident. Hotspots caused by these defects limit platefueled research reactor performance. Accurate modeling of hotspots is required for qualification of new LEU fuels and also improves the safety analysis and operations for current HEU fuels. A portion of this effort to improve hotspot models is to accurately model non-bonds, which is the focus of this paper.

2. Numerical Model For this study, a simplified 2D model of a fuel plate was developed, shown in Figure 1. The cladding is 38.1 mm [1.5 in] long and 1.27 mm [0.05 in] thick. The fuel is centered within the cladding and is 25.4 mm [1 in] long and 0.76 mm [0.03 in] thick. A number of simplifying assumptions were made in the development of the model. The first is that the defect-free portion of the fuel experiences uniform energy generation. In reality, the configuration of the fuel plate within the reactor will cause variations in energy generation across the fueled region. A second assumption is that a constant convective heat transfer coefficient exists with a constant coolant temperature adjacent to the fuel plate. In reality, the fuel plate is submerged in a distributed turbulent flow of coolant that warms as it flows across the plate. The built-in material properties for aluminum were used for both cladding and fuel. Heat generation rates and the convective heat transfer rate were set to 1.63e10 W/m3 and 85.2 W/(m2 K) respectively. The fluid temperature was set to 322 K [120°F] for the convective boundary condition. The left and right edges of the model were insulated.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

Figure 1 Simplified fuel plate model including a fuel segregation (green), altered boundary (red), and a nonbond/adiabatic disc (black) [Note that the vertical scale is significantly enlarged].

The meshes were generated using the semiconductor settings with free quadrilaterals, and the built-in mesh-density sizes.

thermal contact, that may be used to model a non-bond with reduced mesh requirements. 3.1 Thin Layer

3. Use of COMSOL Multiphysics® Software Nuclear fuel plates present a conjugate heat transfer challenge. Energy produced within the fuel plate is conducted through the fuel plate and into the flowing coolant. The primary equations describing these processes are an energy balance on both the fuel plate and the coolant, and the continuity equation and Navier-Stokes equations on the coolant. The simplifying assumptions of a constant convective heat transfer coefficient and constant bulk fluid temperature, mentioned above, reduce the equations describing the fluid to a boundary condition in the fuel plate energy balance. The focus of this effort was on modeling non-bonds. Historically, the non-bond has been modeled as an adiabatic disc between the fuel and the cladding [1-4]. This results in a conservative prediction of the peak cladding surface temperature. By more accurately modeling the non-bond, some of this conservativeness may be reduced while still maintaining appropriate safety margins. While the non-bond can be, and has been, modeled using traditional FEA techniques, this requires the generation of a very fine mesh in the vicinity of the non-bond. COMSOL Multiphysics offers two features, thin layer and

For the thin layer and thermal contact models, the bottom of the fuel portion, the red line in Figure 1, was used as the altered boundary. For these features to be of value, the temperature fields and heat fluxes should approach what would be expected without a modified boundary. To make this happen, the heat flux through the altered boundary had to be as similar to the unaltered boundary as possible. The thin layer is mathematically modeled as 𝑇𝑢 − 𝑇𝑑 −𝒏𝑑 ⋅ 𝒒𝑑 = −𝑘𝑠 𝑑𝑠 where the d and u subscripts indicate the “up” and “down” side of the boundary, ks is the thermal conductivity of the layer, and ds is the thickness of the boundary. The up and down side of the boundary represent the same location and operate at co-incident nodes between two domains. When the boundary is between two materials with different thermal conductivities the geometric mean of the two conductivities is used. Within COMSOL Multiphysics the options available are to use the material properties of the two surfaces to determine the conductivity or to provide a user-defined thermal conductivity. For this model, the built-in material properties were used, and a thin layer thickness of 1 μm [39 μin] was used.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

3.2 Thermal Contact COMSOL Multiphysics offers two models for thermal contact conductance: CMY for plastic contact and Mikic for elastic contact. This study used the Mikic model because the non-bond contact is elastic. The model is represented mathematically by −𝒏𝑑 ⋅ 𝒒𝑑 = −ℎ𝑐 (𝑇𝑢 − 𝑇𝑑 ) + 𝑟𝑄𝑏

(1)

there were complete thermal contact, this coefficient must equal 1. Uncertainty in the values for 𝑚𝑎𝑠𝑝 and 𝜎𝑎𝑠𝑝 typical for a non-bond, suggested they be examined in a later sensitivity analysis. For this study, values of 1 μm [39 μin] for the asperity 𝑝 height, a slope of 0.2 and the corresponding ′ of 𝐸 0.4950 (to bring 𝐶𝑡𝑐 to 1) were used.

4. Method

with

0.94

𝑚𝑎𝑠𝑝 √2𝑝 ( ) . 𝜎𝑎𝑠𝑝 𝑚𝑎𝑠𝑝 𝐸 ′ where masp is the average slope of the surface topology, σasp is the average asperity height, or surface roughness, E’ is the reduced modulus of the two materials, kc is the geometric mean of the thermal conductivity of the two materials, and p is the contact pressure between the two surfaces. Qb is for heat generation within the boundary. As explained above, in order to reduce the impact of the altered boundary on the temperature field, the heat transfer through this thermal contact boundary should equal that through an unaltered boundary. The following steps were taken to achieve this result. By rearranging the conductance relationship, hc becomes ℎ𝑐 = 1.54𝑘𝑐

0.94

𝑘𝑐 √2𝑝 ℎ𝑐 = [1.54 ( ′ ) ] . 𝐸 𝜎𝑎𝑠𝑝 By substituting this into (1), and by assuming heat transfer can only occur normal to the contact and that there is no heat generation within the boundary, this can be rewritten as 0.06 𝑚𝑎𝑠𝑝

−𝑞𝑡𝑐 = [1.54

0.06 𝑚𝑎𝑠𝑝

√2𝑝 ( ′ ) 𝐸

0.94

]

𝑘𝑐 (𝑇 − 𝑇𝑑 ). 𝜎𝑎𝑠𝑝 𝑢

By letting 𝐶𝑡𝑐 = 1.54

0.06 𝑚𝑎𝑠𝑝

This becomes

√2𝑝 ( ′ ) 𝐸

0.94

𝑇𝑢 − 𝑇𝑑 . 𝜎𝑎𝑠𝑝 Recognizing that 𝜎𝑎𝑠𝑝 , the mean asperity height of the surfaces in contact, is the un-modeled length separating the up and down temperature, it can be seen that this is essentially Fourier’s law with an added coefficient, Ctc. This coefficient can be thought of as the modification to heat transfer due to the constrictions imposed by the contacting surfaces. In order to make this heat transfer equal to what would take place if −𝑞𝑡𝑐 = 𝐶𝑡𝑐 𝑘𝑐

The purpose in examining the thin layer and thermal contact features was to evaluate their ability to simulate a non-bond between the fuel and cladding. To this end, an adiabatic disk was introduced as a model for a non-bond in the thin layer boundary. This mimics how non-bonds have been modeled in the past. The adiabatic disc was produced by using a function to provide the user-defined thermal conductivity within the thin layer. The center 1.59 mm [0.0625 in], a length corresponding to the maximum acceptance criteria for non-bond inspections, of the function was set to 0 W/(m K) while the remainder of the length was set to the thermal conductivity of the material. A 0.254 mm [0.01 inch] transition, an arbitrarily chosen length, was included to soften the otherwise discontinuous thermal conductivities. The transition maintained a continuous second derivative. The non-bond was similarly simulated in the thermal contact boundary: a function was developed where 1.59 mm [0.0625 in] at the center of the boundary used the realistic values of 3.1 MPa [450 psi] for the contact pressure and 70 GPa for the reduced modulus, while the remainder of the boundary used the same 0.4950 𝑝 ratio for ′ identified above. Again, a 0.254 mm 𝐸 [0.01 in] transition region was used between these two values. This function of pressure over elastic modulus was then supplied in the “pressure” field and 1.0 was supplied in the “elastic modulus” field of the thermal contact menu. To simulate the non-bond directly in the base model, a geometric entity with the same length and a 1 μm thickness was inserted at the center of the plate immediately below the fuel. The non-bond was assigned aluminum properties, as well, though the thermal conductivity set to 0.001 W/(m K) to model extremely poor conductivity.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

A fuel segregation, which introduced higher heat flux as well as a region with significantly lower thermal conductivity, was modeled with the addition of a geometric feature. The built in aluminum properties were again used, except that the thermal conductivity was lowered to 25% of the value for aluminum. The energy generation rate for the fuel segregations was set to 10.3 times the fuel generation rate. At 0.417 mm [.0164 in] long and the thickness of the fuel, the fuel segregation was located at the center of the model and is shown in green in Figure 1. A number of cases were run (1) to observe how the inclusion of TL and TC boundaries without non-bonds or fuel segregations change the results of the base model, (2) to compare the results of models with non-bonds created with TL, TC and by conventional FEA methods, and (3) to observe the effects of the inclusion of fuel segregations both independent of and coincident with non-bonds.

Energy balance errors were evaluated by subtracting the energy leaving through the boundaries of each domain from the energy generated within those boundaries, and dividing that by the energy generated. The base model showed no error in energy balance. Both TC and TL showed very similar errors to one another (identical to >4 digits of accuracy), shown in Table 1. Within this table, Fine, Finer, and Extra Fine refer to the built-in sizes for mesh generation.

5. Results

The normal heat flux along the modified (TL or TC) boundary was very similar to that seen from the base model, except near the edges, where the modified boundaries oscillated and terminated at 0 while the base remained smooth and did not drop to 0, Figure 3.

The analysis examined the temperature distribution, the energy balances, and behavior of the heat fluxes within each model. The visual appearance of the results was very similar, and in some cases indistinguishable between the three modeling options (base, thin layer, and thermal contact).

Table 1: Errors in energy balance for both TL and TC at different grid refinements for the base case without fuel-defects or non-bonds. Grid Fine

Finer

Extra Fine

Fuel

0.0759%

0.0217%

0.0107%

Whole Model

0.0000%

0.0000%

0.0000%

Domain

5.1 Models without Fuel Defects All three models without the presence of fuel defects produced visually indistinguishable temperature fields, Figure 2. By subtracting the results of the base model from the results of the models with the thin layer or thermal contact boundaries, a maximum difference of 1e-2 K was observed.

Figure 3 Comparison of the boundary, normal heat flux along the modified boundary for the base case without fuel segregations or non-bonds.

5.2 Models with Simulated Non-bonds

Figure 2 Temperature field of base case without fuel segregations or non-bonds.

Introduction of a non-bond produces a localized hotspot at the surface of the fuel plate opposite from the non-bond. The intensity of the hotspot depends on the conductivity through the nonbond. Figure 4 shows the temperature increase caused by an adiabatic non-bond. The figure was generated by subtracting the results of the defect-

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

free base model from the TL with an adiabatic disc model.

altered boundary, Figure 7, are not as severe as with a non-bond. This plot shows just to the left of center of the fuel plate where the discontinuity in conductive heat flux, due to the discontinuity in thermal conductivity between the fuel segregation and the fuel, is captured.

Figure 4 Temperature increase caused by an adiabatic non-bond.

Table 2 shows the maximum temperatures predicted and the maximum increase in temperature over the defect-free base. Table 2: Maximum temperature (K) predicted and increase in temperature over the defect-free base caused by an adiabatic (wad) or a conductive (wnb) non-bond. Base wad

TL wad

TC wnb

max

393

393

390

increase

14.4

15

11.6

Introduction of a non-bond, either adiabatic or conductive, causes only minor changes ( a Here A is a given constant, ˆı3 is a constant unit vector parallel to the axis of symmetry, and r is

(rv)T ]a = curl v ⇥ a

(ru)T = W .

(1.4)

In view of (1.3) W must satisfy W (a) = [ru (ru)T ](a) = curl u ⇥ a = [(Aˆı3 ) ⇥ r] ⇥ a. It follows from the expansion of the vector triple product and the definition of the tensor product* that W = A(r ⌦ˆı3

ˆı3 ⌦ r) .

(1.5)

2. Variational principle for the flow inside Hill’s spherical vortex In Ri consider the problem of minimizing the expresson F defined by ZZZ {(div u)2 +(1/4)kru (ru)T W ]k2 } dV F := Ri

(2.1) over all smooth functions r 7! u. Here dV is the volume of a typical small part of Ri . Equation (2.1) employs the notation kSk2 := S • S for the squarenorm of a tensor, S. The inner product, A • B, of tensors A and B is defined by the rule A • B := * For any two given vectors, a and b, their tensor product, denoted a ⌦ b, is a linear vector-to-vector operator whose action upon a third vector, v, is defined by the rule (a ⌦ b)v = a(b • v).

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

COMSOL Conference Boston 2016

page 2

tr (AT B), in which tr ( . ) denotes the trace operator and ( . )T denotes the transpose. Note that F is nonnegative by construction. Moreover, F attains the value zero when u satisfies equations (1.1) and (1.4). By requiring the first variation, F , of F to be zero for arbitrary variations, u, of u within Ri and on @Ri one my derive a the di↵erential equation for u in Ri (a.k.a. the Euler-Lagrange equation) for the variational problem and a natural (a.k.a. fluxsource) boundary condition for u on @Ri . To be specific, if one takes the first variation of (2.1), applies some identites*, and introduces the abbreviation := 2(div u)I + ru

(ru)T

W

(2.2)

one may arrange the result in the form F

ZZZ

•

r( u) dV = 0 .

(2.3)

Ri

If one sets F to zero (as is appropriate when F attains a minimum) and interprets the first variation, u, of u as what the COMSOL documentation denotes by v and calls a test function then equation (2.3) becomes what that documentation calls the weak form. By application of some more identities† one may write (2.3) in the equivalent form F+

ZZ

@Ri

u •[

(ˆ n)]dA +

ZZZ

surface. If F is a minimum then F = 0 subject to arbitrary variations u on @Ri and in Ri . Equation (2.4) then implies that div

ˆ=m, u•n

Ri

* Three such identities are: (i), the divergence of a vector field equals the trace of its gradient; (ii) the trace of any tensor equals the inner product of that tensor with the identity, I; and (iii), the inner product of two tensors equals zero whenever one factor is symmetric and the other is skew. † One such identity [cf. Gurtin (Ref. 3), equation (4.2)5 on p30] is the di↵erentiation formula div [S T (v)] = S • rv + v • div S (in which v and S are generic vector and tensor fields); another is the divergence theorem; and third is the definition of the transpose, according to which S T (a) • b = a • S(b).

(ˆ n) = 0 ,

(2.5)1,2

which constitute the Euler-Lagrange equation and the natural boundary condition, respectively. The system consisting of (2.2) and (2.5)1 is suitable for input to the COMSOL general form PDE physics interface and (2.5)2 is the corresponding null flux boundary conditon. The solution for u of the minimization problem as posed thus far is not unique. To see why note that replacement of u by u + r' will leave the expression under the integral sign in (2.1) unchanged provided rr' (rr')T = O—which is always the case—and div (r') = 0. There are, of course, infinitely many solutions of div (r') = 0 in Ri so if there is one vector field u in Ri for which F is zero there must be infintely many of them. To remove this ambiguity one may recall that the most general solution of the problem consisting of the equation div (r') = 0 subject to the boundary ˆ = 0 in a simply connected domain condition r' • n (a.k.a. the homogeneous Neumann problem) is a constant. But the gradient of a constant is the zero vector so r' must reduce to 0 in that case. Now a ˆ in classical boundary condition that specifies r' • n ˆ potential theory corresponds one that specifies u • n in the problem of minimizing F , i.e. to a boundary condition of the form

u • div dV = 0 .

(2.4) Here dA is the area of a typical small part of @Ri and n ˆ is the outward unit normal vector on that

=0 ,

(2.6)

in which r 7! m is subject to the constraint that the total volumetric outflow across @Ri be zero (as incompressiblity requires) but is otherwise arbitrary. If the observer for whom u is the fluid velocity is at rest relative to the boundary sphere then m in (2.6) is identically zero.

3. Variational principle for the flow outside Hill’s spherical vortex In Re consider the problem of minimizing the expresson F defined by ZZZ ⇥ ⇤ (div u)2 + (1/4)kru (ru)T ]k2 dV F := Re

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

(3.1)

COMSOL Conference Boston 2016

page 3

over all smooth functions r 7! u. The fact that Re is unbounded poses both practical and mathematical challenges. One way of addressing these challenges is precede the treatment of the minimization problem for F by a change of independent variable r ! q, which maps points in the original domain, Re —which is not bounded—to points in proxy domain, Q, which is. To this end, let a > 0 be a given constant length and let r and q denote the magnitudes of the vectors r and q, respectively. One may then define a change of variable r ! q by the following equations: r/r =

q/q

,

a2 = rq .

(3.2)1,2

One may visualize r and q as position arrows springing from a common point, namely the center of the boundary sphere r = a or q = a. If the minus sign in (3.2)1 were replaced by a plus sign then the resulting change variable r ! q would represent Kelvin inversion, as introduced by Kelvin in 1845 (Ref. 1) and employed in classical treatises on potential theory since then [e.g. Kellogg (Ref. 4, Chapter IX, §2]. I will postpone the discussion of my reason for introducing the minus sign in (3.2)1 until after I have discussed an orthogonal curvilinear coordinate system that arises naturally in this problem. I introduced the constant unit vector, ˆı3 , in (1.2)1 above. Now let ˆı1 and ˆı2 be any two constant unit vectors chosen such that {ˆı1 ,ˆı2 ,ˆı3 } forms a right-handed orthogonal system. One may expand a generic vector v into P3components with respect to {ˆı1 ,ˆı2 ,ˆı3 }, viz. v = i=1 viˆıi . Here and elsewhere, italic symbols with numerical subscripts denote the scalar components with respect to {ˆı1 ,ˆı2 ,ˆı3 } of a vector denoted by the corresponding letter, without subscripts, in bold face. According to this convention once {ˆı1 ,ˆı2 ,ˆı3 } is fixed the list (q1 .q2 , q3 ) determines q, and the system (3.2)1,2 then determines r, so (q1 , q2 , q3 ) 7! r is now a known function and (q1 , q2 , q3 ) constitutes a set of curvilinear coordinates. One may show that this curvilinear coordinate system is orthogonal, i.e. (@r/@qi ) •(@r/@qj ) = 0 for i 6= j

(3.3)

It will be convenient at this point to recall some textbook results that apply to all orthogonal curvilinear coordinate systems [i.e ones which satisfies

(3.3)]. Given any orthogonal curvilinear coordinate system, one may associate each coordinate, qi , with a corresponding scale factor, hi , defined by hi := k@r/@qi k ,

i 2 {1, 2, 3} .

(3.4)

From (3.3) and (3.4) one may construct a system of unit vectors ˆ ei := (1/hi )(@r/@qi ) i 2 {1, 2, 3} belonging to any orthogonal curvilinear coordinate e2 , ˆ e3 } one may system. Having the system {ˆ e1 , ˆ with expand a generic vector v into components P3 e ei , in respect to that system, i.e. v = i=1 vi ˆ which the scalar components (v1e , v2e , v3e ) are not to be confused with the scalars (v1 , v2 , v3 ) in the expanP3 sion v = i=1 viˆıi . Textbooks that treat orthogonal curvilinear coordinates [e.g. Phillips (Ref. 5, pp 88–90) and Hildebrand (Ref. 6, pp 306–311)] include derivations of expansions of the divergence and curl of a generic vector, v, and of the gradient of a generic scalar, ', with respect to orthogonal curvilinear coordinates. These authors, and others, assume in their derivations that the system e2 , ˆ e3 } is right-handed. It so happens that the {ˆ e1 , ˆ omission of the minus sign in the definition (3.2)1 above would result in a left-handed system. This fact drove my decision to employ an inversion rule motivated by Kelvin’s but not Kelvin Inversion proper. Returning, now, to the textbook results one may write them in the following forms: div v =

✓ ◆ 3 X @ (h1 h1 h3 vie ) 1 , h1 h2 h3 i=1 @qi hi

curl v = and

3 X 3 X ˆ ei ⇥ ˆ ej @(hj vje ) , hi hj @qi i=1 j=1 3 X ˆ ei @' r' = . h @qi i=1 i

(3.5)

(3.6)

(3.7)

One may employ (1.3) to deduce the corresponding representation of rv (rv)T with respect to orthogonal curviliner coordinates as follows. If a is any vector equation (3.6) implies that curl v ⇥ a =

3 X 3 X [(ˆ ei ⇥ ˆ ej ) ⇥ a] @(hj vje ) i=1 j=1

hi hj

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

@qi

, (3.8)

COMSOL Conference Boston 2016

page 4

But (ˆ ei ⇥ ˆ ej ) ⇥ a = (ˆ ei • a)ˆ ej (ˆ ej • a)ˆ ei = (ˆ ej ⌦ ei ⌦ ˆ ej )a, so (3.8) is equivalent to ˆ ei ˆ

so (3.12)–(3.14) are equivalent to

X 3 X 3 (ˆ ej ⌦ ˆ ei ˆ ei ⌦ ˆ ej ) @(hj vje ) (a)=curl v ⇥ a. hi hj @qi i=1 j=1

(3.9) Since the right members of (1.3) and (3.9) are equal their left members must be equal as well for all a. The operators that act upon a must then be equal, i.e. rv

h1 = h2 = h3 := h = (a/q)2 .

(3.11)

Equations (3.5), (3.10), and (3.7) then take the simpler forms div v =

rv

rv

(rv)T =

3 1 X @(h2 vie ) , h3 i=1 @qi

3 3 1 XX T (rv) = 2 (ˆ ej ⌦ ˆ ei h i=1 j=1 3 X ˆ ei @' , r' = h @qi i=1

(3.13) (3.14)

respectively. Let Q denote the linear vector-to-vector operei , i 2 ator—i.e. the tensor—that takes ˆıi to ˆ {1, 2, 3}. Then ˆ ei = Q(ˆıi ). Since Q takes unit vectors to unit vectors, it must have the feature that the magnitude of its input must equal the magnitude of its output. Such an tensor is called orthogonal and has the feature that its transpose equals its inverse. In view of the identities ˆ ei = Q(ˆıi ) we have ei • v = Q(ˆıi ) • v = ˆıi • QT (v) , vie = ˆ

(3.15)

@[ˆıj • QT (hv)] , @qi

(3.17)

3

r' =

@' 1X Q(ˆıi ) , h i=1 @qi

(3.18)

respectively. If a and b are any two vectors and S is any tensor then two rules of tensor algebra state that S(a ⌦ b) = S(a) ⌦ b and (a ⌦ b)S = a ⌦ S T (b) [See e.g. Gurtin (Ref. 3, p9, Exercises 6 a,b)]. From these general rules one may deduce, for example, that Q(ˆıj ) ⌦ Q(ˆıi ) = Q(ˆıj ⌦ ˆıi )QT , Thus (3.17) is equivalent to rv

X 3 X 3 1 (ˆıj ⌦ˆıi (rv) = 2 Q h i=1 j=1 T

ˆıi ⌦ˆıj )

(3.12)

@(hvje ) ˆ ei ⌦ ˆ ej ) , @qi

(3.16)

3 3  1 XX Q(ˆıj ) ⌦ Q(ˆıi ) h2 i=1 j=1

Q(ˆıi ) ⌦ Q(ˆıj )

3 X 3 X (ˆ ej ⌦ ˆ ei ˆ ei ⌦ ˆ ej ) @(hj vje ) . (rv)T = hi hj @qi i=1 j=1

(3.10) The discussion in the paragraphs containing (3.4) through (3.10) applies to all orthogonal curvliner coordinates systems. I will now revert to the discussion of the particular orthogonal curvilinear coordinate system defined by (3.2)1,2 and the system {ˆı1 ,ˆı2 ,ˆı3 }. Here, one finds that

3 1 X @[ˆıi • QT (h2 v)] , h3 i=1 @qi

div v =

@[ˆıj • QT (hv)] T Q . (3.19) @qi

One may express these results more compactly by introducing some abbreviations. To this end, let 3 X @vi i=1

@qi

=

3 X @(ˆıi • v) i=1

@qi

:= divq v ,

3 X @' ˆıi := rq ' , @qi i=1

(3.20)

(3.21)

3 X 3 3 3 X @vi X X @(ˆıi • v) (ˆıi ⌦ˆıj ) = (ˆıi ⌦ˆıj ) := rq v . @qj i=1 j=1 @qj i=1 j=1

(3.22) Then equations (3.16), (3.19), (3.18) takes the more compact forms div v =

1 divq [QT (h2 v)] , h3

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

(3.23)

COMSOL Conference Boston 2016 rv =

page 5

(rv)T

1 Q{rq [QT (hv)] h2 r' =

rq [QT (hv)]T }QT ,

(3.24)

1 Q(rq ') , h

(3.25)

respectively. Returning, now, to the variational principle (3.1) one may transform the right member from an integral over the physical exterior domain, Re , to one over the proxy domain, Q, by application of the transformation rules (3.23) and (3.24) in the case when the generic vector, v, is the fluid velocity, u. In carrying out this transformation the volume element transforms according to the rule dV = h3 dVq , in which dVq is the volume of a typical small part of Q. The simplicity of this rule is due to both the orthogonality of the coordinates and the equality of three scale factors, as stated in (3.11). Equation (3.1) thus becomes F :=

ZZZ 

1 Q{rq [QT (hu)] rq [QT (hu)]T }QT + 4 h2

2

#

h3 dVq , (3.26)

which one may simplify in three ways, namely: (i) by cancellation of common powers of h in the numerator and denominator under the integral sign; (ii), by introducing the change of variable u ! U defined by QT (hu) := U; and (iii), by appealing to an algebraic rule according to which kQSQT k2 = kSk2 for any tensor, S, and orthogonal tensor, Q. The resulting simplified form of (3.26) is F :=

1 {divq (hU)}2 h3

Q

+

1 [rq (U) 4h

rq (U)T ]

Q

1⇥ + rq U h

(rq U)T

⇤

•

rq ( U) dVq = 0 . (3.28)

Three identities are useful here. The first asserts that r('v) = 'rv + v ⌦ r' for a generic scalar ' and vector v [cf. Gurtin (Ref. 3, equation (4.2)1 , p30] and the second and third assert that I •(a ⌦ b) = tr (a ⌦ b) = a • b [cf. Gurtin (Ref. 3, pp 5– 6)]. With the aid of these identities one may deduce that I • rq (h U) = hI • rq ( U) + U • rq h .

1 {divq [QT (h2 u)]}2 h6

Q

ZZZ 

equation for U in Q (a.k.a. the Euler-Lagrange equation) for the variational problem and a natural (a.k.a. flux-source) boundary condition for u on @Q. To be specific, if one takes the first variation of (3.27) and applies some identites as described in the footnote prior to equation (2.2) one may arrange the result in the form ZZZ  2 div (hU)I • rq (h U) F h3

2

dVq . (3.27)

At this point the reasoning follows logic similar to that described in the text containing (2.1)– (2.5) above. By requiring the first variation, F , of F to be zero for arbitrary variations, U, of U within Q and on @Q one my derive a the di↵erential

(3.29)

If one substitutes (3.29) into (3.28) and introduces the abbreviations ⇤ 2 1⇥ := 2 div (hU)I + rq U (rq U)T (3.30) h h 2 div (hU)rq h h3 that equation becomes ZZZ ⇥ ⇤ • r ( U) + f • U dV = 0 F q q f :=

(3.31)

(3.32)

Q

If one sets F to zero (as is appropriate when F attains a minimum) and interprets the first variation, U, of U as what the COMSOL documentation denotes by v and calls a test function then equation (3.32) becomes what that documentation calls the weak form. By application of the identities described in the footnote preceding (2.4) above one may write (3.32) in the equivalent form ZZ U •[ (ˆ nq )]dAq F+ @Q

+

ZZZ

U •(divq

f ) dVq = 0 . (3.33)

Q

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

COMSOL Conference Boston 2016

page 6

Here dAq is the area of a typical small part of @Q and n ˆq is the outward unit normal vector on that boundary surface. If F is a minimum then F = 0 subject to arbitrary variations U on @Q and in Q. Equation (3.33) then implies that divq

=f

,

(ˆ nq ) = 0 ,

(3.34)1,2

which constitute the Euler-Lagrange equation and the natural boundary condition, respectively. The system consisting of (3.30), (3.31) and (3.34)1 is not quite suitable for input to the COMSOL general form PDE physics interface owing to the presence of the derivative of a product, namely divq (hU). If one substitutes h = a2 /q 2 [cf. (3.11) above] and evaluates the indicated derivatives one finds that (3.31) and (3.32) are equvialent to ◆ ✓ 2 ⇤ 4 q2 ⇥ 2q • div U q U I + rq U (rq U)T = q a2 a2 a2 (3.35) and ◆ ✓ 2 4 q • (3.36) q divq U q U . f= a2 q q The system consisting of (3.34)1,2 , (3.35), and (3.36) is now suitable for input to the COMSOL general form PDE physics interface.

4. Conditions on normal and slip velocities; Results As noted in the discussion of the flow in Ri [see the paragraph containing (2.6) above] the EulerLagrange equation and natural boundary condition are not sufficient, by themselves, to ensure uniqueness of u, or, in the present context, U. One can determine U uniquely by augmenting the natural boundary condition (3.34)2 , with a condition analogous to (2.6), namely a specification of the ˆq on @Q. If the observer for whom value of U • n u is the fluid velocity in Re is at rest relative to the boundary sphere, @Re , then the impermeableˆ = 0. sphere boundary condition takes the form u • n Nontrival solutions are possible if the fluid at a remote distance from the sphere is in uniform motion, with, say, a downward fluid velocity w1ˆı3 with w1 < 0. in that case U = QT (hu) = QT [(a2 /q 2 )w1ˆı3 ], whose rightmost member goes to

infinity as q ! 0. Such a behavior for the unknown U is unsatisfactory for computation. A alternative approach is to frame the impermeable-sphere boundary condition on @Re under the assumption that the observer for whom u is the velocity in Re is at rest relative to the remote undisturbed fluid. In such a reference frame the bounding sphere, @Re , translates upward with velocity w1ˆı3 and the impermeable sphere boundary condition beˆ = ( w1ˆı3 ) • n ˆ on @Re . comes u • n • ˆ at a typical point on @Re to To relate u n • ˆq at its image on @Q let be a scalar-valued U n function of position with the feature that both @Re and @Q are surfaces of constant , say = 0. For definiteness let be negative valued in both Re and Q. Then the outward unit normal vectors on @Re and @Q satisfy n ˆ = (r /kr k)

=0

, n ˆq = (rq

/krq

k) =0 , (4.1)1,2

respectively. If one re-expresses r in the right member of (4.1)1 by means of (3.25), cancels the factors 1/h in the numerator and denominator, and notes that kQ(a)k = kak for every vector a and orthogonal tensor, Q, equation (4.1)1 becomes ⇥ n ˆ = Q(rq

)/krq

⇤ k

=0

.

(4.2)

If one substitutes (4.1)2 into the right member of (4.2) that equation becomes n ˆ = Q(ˆ nq ) .

(4.3)

But QT (hu) = U, so u = (1/h)Q(U). If one evaluates this last equation on the boundary sphere q = a and notes that h = 1 there one gets u = Q(U)

(4.4)

Now the inner product of the left members of (4.3) and (4.4) must equal the inner product of their right members, so ˆ = Q(U) • Q(ˆ nq ) = U • n ˆq , u•n

(4.5)

in which the last equality follows from the orthogonality of Q. In view of (4.5) the impermeable-sphere

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

COMSOL Conference Boston 2016

page 7

boundary condition u • n ˆ = ( w1ˆı3 ) • n ˆ on @Re is equivalent to ˆq = ( w1ˆı3 ) • n ˆ. U•n

(4.6)

The system consisting of the di↵erential equation (3.34)1 subject to the natural boundary condition (3.34)2 and the impermeable-sphere condition (4.6) is now sufficient to ensure a unique solution for U. Having U in Q, one could compute u in Re from u = (1/h)Q(U), which would represent the velocity relative to the remote undisturbed fluid. If, alternativley one wants u to represent the velocity as seen by an observer who moves with the sphere then one may take u = (1/h)Q(U) + w1ˆı3 . The vorticity constant, A, for the flow in Ri and the free-stream velocity, w1ˆı3 cannot be specified independently without causing a discontinuity in the tagential velocity, or slip, across the boundary sphere. To avoid this slip one specifies only one of the two parameters A and w1 and computes the other. The present COMSOL model implements no-slip by requiring that the average around the equator of the ˆı3 -component computed from the usolutions in Ri and Re agree. The present COMSOL model employs three components: the first contains a General Form PDF interface to solve for u in Ri ; the second contains a General Form PDF interface to solve for U in Q; and the third calculates u in a subrgion of Re —namely the region a < r < 2a—from U in the corresponding subregion of Q. This third component enables graphical illustration of u in the immediate neighborhood of the spherical vortex. To accomplish this calculation from the identity u = (1/h)Q(U) one requires an explicit representation of Q, namely Q = 2(q/q) ⌦ (q/q) I. These three components employ General Extrusion Model Coupling Operators to exchage information between components. The first component contains a node which enforces no-slip by means of a Global ODE and DAE node. There are two Study Steps: the first carries out the computations in components 2 and 3; and the second carries them out in component 1. Fig. 5.1 nearby illustrates the results.

5. Conclusion COMSOL’s General Form PDE interface enables computation of the three-dimensional velocity field

Figure 4.1 Computed velocity field in Hill’s spherical vortex in a cut plane containing the axis of symmetry: shading represents fluid speed; solid lines represent streamlines; and arrows represent velocity vectors.

of a fluid flow with nonzero vorticity in a bounded subdomain of a larger domain that lacks any exterior boundary.

6. References 1. Kelvin, Lord Extrait d’une lettre de M. William Thomson `a M. Liouville. Journal de Math´ematique Pures et Appliqu´ees, 10, 1845. In Reprint of Papers on Electrostatics and Magnetism by Sir William Thomson, second edition, Cambridge, pages 144–146, 1884. 2. Hill, M.J.M. On a Spherical Vortex. Philosophical Transactions of the Royal Society of London (A), 185, pp 213–245, 1894. 3. Gurtin, Morton B. An Introduction to Continuum Mechanics. Academic Press, 1981. 4. Kellogg, O.D. Foundations of Potential Theory, Springer, 1929. 5. Phillips, H.B. Vector Analysis, Wiley, 1933. 6. Hildebrand, F.B. Advanced Calculus for Applications, Second edition, Prentice-Hall, 1976.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

The Use of Finite Element Analysis in the Design of Oil-Water Separators Mohamed E. Wanas1*, Yehia M. S. Elshazly1, Dina A. ElGayar1 1 Chemical Engineering Department, Faculty of Engineering, University of Alexandria; Alexandria; Egypt *

[email protected], [email protected]

Abstract: Oil-Water separators play an important role in several industries as well as in waste water treatment. However, no basic principles have been set to guarantee the designed separators would work according to the desired efficiency due to the effect various factors. Both time and money could be saved by simulating the separation process using Comsol software without the need of building a prototype and testing it. This study aims to use Comsol software to simulate the separation process and test its accuracy by comparing the modeling and experimental results. Upon confirmation of modeling accuracy, the study also aims to test the effect of different operational and design parameters on the efficiency of the separation process. Different conclusions drawn from this study could be very helpful in the design and operation of the separators. Keywords: Computational Fluid Dynamics, Two phase flow, Comsol, oil- water Separators.

1. Introduction Separation processes play an important role in several industries such as in oil and gas refineries. Gravity settling separators are used widely in oil and gas industry to separate water from oil produced from wells. These separators depend mainly on the difference in specific gravity between the two fluids[1, 2]. As the oil wells get older the oil-water separators' importance becomes even greater in order to overcome the oil fields aging phenomenon by which the water content in the oil produced from wells increases [3]. This increase causes corrosion problems to the piping and equipment used in the production and reduces the heating value of the produced oil.[4] Another important use of gravity settling separators is in the treatment of wastewater produced from different industries or from ships bilge tanks which contain mixtures of different fluids along with water that could have harmful

effects on the environment; hence, it has to be separated prior to the disposal to comply with environmental laws and regulations.[5-7] Thus, the design of separators has achieved great attention as many factors interfere in the separation process. The design of such separators have only depended on previous experiences or study prototypes to achieve the desired separation efficiency.[8-10] In the past few years, Computational fluid dynamics (CFD) and finite element analysis have been used as a replacement to the previous methods. CFD has many advantages such as cost, time and effort reduction over traditional methods. A major advantage is also that not all separator prototypes are easy to build or perform experimentally; however, these experiments could be easily simulated or applied using CFD modeling.[11-17] In this study, Comsol was used to build horizontal oil-water separators and study the effect of different variables and designs on the separation process efficiency. Similar separator prototypes were built and tested at similar settings to those in the Comsol models in order to test their accuracy compared to the experimental results. Moreover, a tracer (Methylene blue dye) was injected into the prototype using a syringe from the inlet nozzles and the color was monitored to study the flow pattern.

2. Laboratory Experiments The laboratory experiments used a tank equipped with a rotating stirrer to mix transformer oil (who's viscosity was determined using a viscometer and found to be 0.0603Pa.s and a Pycnometer to measure the Specific gravity which was found to be 0.870) and tap water. A 0.5hp centrifugal pump was then used to pump the fluid into the horizontal separator. Valves controlled the flow which was measured using a rotating vane flow meter. (Fig.1) The separator is a plexi glass cylindrical horizontal vessel of 45cm length and a diameter

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

of 15cm [18]. It is equipped with half elliptical ends on both sides. The separator had one inlet nozzle at the center of one of its ends and two outlet nozzles: one at the top and one at the bottom, near the other elliptical end. The inlet and outlet nozzles had an inner diameter of 1.5cm. Samples were collected at the inlet and the two outlets for measurement (Fig.2). The samples were left to settle in graduated measuring cylinders to measure their oil in water volumetric fractions.

Reynolds number calculation confirmed the flow was turbulent, moreover, the presence of several obstacles in the system which would increase the level of turbulence. Thus, the Turbulent Two-Phase Flow, Level Set interface was selected to track the separation between the two liquids.

3.1 Separator Geometry The separator model was built in similar dimensions to those used in the laboratory prototype. Vanes were added to some models at different location to study their effect on the process.

Figure 3. 3D separator model with inlet diverter.

Figure 1. Process Flow Diagram for the prepared laboratory experiment

Figure 2. Laboratory setup for the experiment before placing into operation

Although 3D simulation would generate more precise results for cylindrical and spherical designs than 2D simulations. Moreover, the 2D model does not describe the actual cylindrical shape accurately. However, the results were close to those of the experiments. Therefore, there was no need for the excessive computer. Only 2D simulation was used in the study.

Figure 4. The simplified 2D separator model.

3. Modeling Experiments

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

3.2 Governing Equations The chosen Turbulent Two-Phase Flow, Level Set interface depends mainly on some basic concepts and equations: The conservation of mass and momentum. The software then solves Navier-Stokes equations in order to the track the interface between the two non mixing fluids.[19]

where u is the fluid velocity, p is the fluid pressure, ρ is the fluid density, μ is the fluid dynamic viscosity, I is the identity matrix and F is the external forces acting on the fluid. The equation is built on the equalization between the different forces acting on the fluid: the inertial forces on one side and the sum of the pressure forces, viscous forces, and the external forces applied to the fluid on the other side. The turbulence effects are modeled using the standard two-equation k-ε model with realizability constraints.[20, 21] Time dependent simulation was used to track the interface until the steady state conditions of operation is reached.

c.

Separator internal design i. Presence of baffles The effect of the presence of baffles (as inlet diverters or bottom baffles) on the separation efficiency has been studied. Experimental results were compared with modeling results to ensure the relationship is observed correctly.

a) Empty separator with no internal baffles Inlet diverter

b) Separator with an inlet diverter Bottom baffle

4. Variables studied The relation between variables and the separation efficiency was monitored by comparing the bottom outlet sample compositions. Different variables studied were: a. Inlet Composition Fluids of different inlet mixture compositions 10%, 20% and 30% by vol. and the out coming fluids from the bottom outlet nozzles were analyzed to study the effect of the inlet composition on the separation process. b. Inlet Velocity Mixtures of a constant composition of 30% oil in water by Vol. were used at various inlet velocities in a separator with no internal baffles to study the effect of the inlet velocity on the separation process. The inlet velocities ranged between 0.2 to 1 m/s in laboratory experiments and between 0.5m/s and 3m/s in Comsol simulations. The velocity in laboratory experiments only reached 1m/s as this was the maximum velocity that could be reached by the used pump's power.

c) Separator with bottom baffle placed 10cm ahead of bottom outlet nozzle. Figure 5. The three different separator designs.

ii. Baffle location inside the separator The location of the bottom baffle ahead of the bottom outlet nozzle was then studied to test its effect on the composition of the bottom outlet.

5.Identification of the flow pattern Another step taken in order to confirm the accuracy of Comsol multiphysics simulations has been done by injecting a tracer (Methylene blue dye) through the inlet nozzle of different separator designs and observing how the colored fluid flows then inside the separator. Comsol was then used to monitor the fluid flow pattern in similar separator models.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

6.Mesh design

Modeling (%)

Practical (%)

Different mesh sizes were tested. The results showed that the size of the mesh is irrelevant. Thus, the Coarse mesh size was used for the rest of the simulations.[22]

Reduction in oil (%)

100 80 60 40 20 0 0

0.5

1

1.5

2

2.5

3

3.5

Inlet velocity (m/s)

Figure 6. 2D discretization of the separator to a coarse mesh size.

Figure 7. Percent reduction in oil results for experiments using an empty separator.

7. Results 7.1 Effect of different variables

7.1.2 Effect of separator internal design

7.1.1 Effect of Inlet velocity The results from laboratory experiments were close to those from the Comsol simulations. Fig.7 also shows clearly that as the inlet velocity increases the % reduction in oil decreases which also indicates a drop in the separation process efficiency.

Fig. 8 indicates that experimental tests separators gave close results as those obtained from their corresponding models. Moreover, the figure also shows that both the inlet diverter and bottom baffle had a positive influence on the separation efficiency over the other separator with no baffles inside. Although, the effect of both baffles seemed very much alike from Fig.8 which could be also confirmed from Comsol simulations in Fig.9. However, the location had another influence, the bottom baffle resulted in greater amounts of water to escape from the upper outlet along with the oil than when using the inlet diverter. Practical No vane

Practical with Inlet diverter

Practical 10cm Bottom Baffle

No vane separator Model

Model with inlet diverter 100

Reduction in oil (%)

The first step taken to study the effect of different variables was developing an equation that would help monitoring the separation process efficiency and especially for different inlet composition where it is not logical to just compare the compositions of the bottom outlet samples. The percentage reduction in oil in the bottom outlet was considered as an indication for the separation efficiency. It was calculated as follows:

Model 10cm Bottom Baffle

80 60 40 20 0 0

0.5

1

1.5

2

2.5

3

3.5

Inlet velocity (m/s)

Figure 8. Percent reduction in oil results for separators of different internal designs.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

Practical No vane separator

Practical 10cm Bottom Baffle

Practical 20cm Bottom Baffle

Modeling No vane separator

Modeling 10 cm Bottom Baffle

Modeling 20 cm Bottom Baffle

100

Reduction in oil (%)

80

60

40

20

0 0

0.5

1

1.5

2

Inlet velocity (m/s)

2.5

3

3.5

Figure 10. Percent reduction in oil results for separators with different bottom baffle locations.

7.1.4 Effect of mixture inlet composition Fig.11 showed that the higher mixture compositions resulted in greater oil composition reductions than the lower mixture compositions. Practical 10% mix Practical 30% mix Modeling 20% mix

Practical 20% mix Modeling 10% mix Modeling 30% mix

100

7.1.3 Effect of bottom baffle location inside the separator Two bottom baffle locations were tested with reference to the bottom outlet nozzle. Locations of the bottom baffles were 10cm and 20cm ahead of the bottom outlet. Their results were compared to those from the no baffle separator. At low velocities the results for both of the baffle locations were close, however, as the velocity increases the one placed closer to the bottom outlet (10cm distance) showed much better separation efficiency than the other one. The effect of the bottom baffle that was distant from the bottom baffle nearly vanished at high inlet velocities.

Reduction in oil (%)

Figure 9. Simulations of different separator designs. 80

60

40

20

0 0

0.5

1

1.5

2

Inlet velocity (m/s)

2.5

3

3.5

Figure 11. Percent reduction in oil results for different inlet mixture compositions.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

7.2 Identification of the Flow Pattern 7.2.1 Separator with no internal vanes Fig. 12 and Fig.13 indicate that the fluid enters with high momentum therefore, it does not spread for a certain distance inside the separator. After the momentum reduces gradually, the blue dye and the color in the model (Fig.13) starts spreading till filling the separator. Finally, the fluid reaches the outlet nozzles where the color by then the color reaches all the ends and fills the entire separator volume. Figure 14. Laboratory experiments using a tracer to identify the flow pattern for a separator with an inlet diverter.

Figure 15. Identification of flow pattern by Comsol simulations for a separator with an inlet diverter. Figure 12. Laboratory experiments using a tracer for identification of flow pattern in an empty no vane separator.

Figure 13. Identification of the flow pattern by Comsol simulations for an empty no baffle separator.

7.2.2 Separator with an inlet diverter Fig. 14 and Fig.15 indicate that the fluid enters with high momentum but starts spreading rapidly to the sides due to the collision with the inlet diverter. The interference of the diverter, decreases the velocity and forms eddies at the center. This causes the color to much longer time to reach the outlet nozzles and fill the separator's entire volume.

8. Conclusions After comparing the Comsol modeling results with the experimental results for different separator designs. A great resemblance was seen between the modeling and the experimental results. Referring to these results we could conclude that: a) Comsol CFD simulations could be used with great confidence to replace the construction of a prototype. b) The Introduction of a vane ( as an inlet diverter or bottom baffle) plays an important role in enhancing the separation efficiency. The introduction of a vane could also offer a cheaper solution that would increase the separation efficiency than other possible solutions such as extending the separator size. c) The location of the bottom baffle inside the separator should be studied with reference to the inlet mixture flow rate to achieve the desired separation efficiency. d) Mixtures with low oil content need longer residence times to reach similar percentages reduction in compositions as those with higher oil contents.

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston

9. References 1. Gas Processors Suppliers Association Engineering Data Book, Twelfth Edition ed., (2004). 2. J.L. Humphrey, G.E. Keller, Separation Process Technology, McGraw-Hill, (1997). 3. H.K. Abdel-Aal, M. Aggour, M.A. Fahim, Petroleum and Gas Field Processing, Taylor & Francis, (2003). 4. H. Devold, Oil and Gas Production Handbook: An Introduction to Oil and Gas Production, (2013). 5. J.H. Vandermeulen, S.E. Hrudey, Oil in Freshwater: Chemistry, Biology, Countermeasure Technology: Proceedings of the Symposium of Oil Pollution in Freshwater, Edmonton, Alberta, Canada, Elsevier Science, (2013). 6. M. Stewart, K. Arnold, Produced Water Treatment Field Manual, Elsevier Science, (2011). 7. C.A.S. Hall, K.A. Klitgaard, Energy and the Wealth of Nations: Understanding the Biophysical Economy, Springer, (2011). 8. E.G. Sinaiski, E.J. Lapiga, Separation of Multiphase, Multicomponent Systems, Wiley, (2007). 9. M. Stewart, K. Arnold, Gas-Liquid And Liquid-Liquid Separators, Elsevier Science, (2008). 10. C.E. Brennen, Fundamentals of Multiphase Flow, Cambridge University Press, (2005). 11. J. Tu, G.H. Yeoh, C. Liu, Computational Fluid Dynamics: A Practical Approach, Elsevier Science, (2012). 12. J. Wendt, J.D. Anderson, G. Degrez, E. Dick, R. Grundmann, Computational Fluid Dynamics: An Introduction, Springer Berlin Heidelberg, (2013). 13. J. Blazek, Computational Fluid Dynamics: Principles and Applications, Elsevier Science, (2015). 14. G. Venkatesan, N. Kulasekharan, S. Iniyan, Numerical analysis of curved vane demisters in estimating water droplet separation efficiency, Desalination, 339, (2014) 40-53. 15. G. Venkatesan, N. Kulasekharan, S. Iniyan, Influence of turbulence models on the performance prediction of flow through curved vane demisters, Desalination, 329, (2013) 19-28. 16. P.W. James, B.J. Azzopardi, Y. Wang, J.P. Hughes, The role of Drainage Channels In the

Performance Of Wave-Plate mist eliminators, Chemical Engineering Research and Design, 81, (2003) 639-648. 17. C. Galletti, E. Brunazzi, L. Tognotti, A numerical model for gas flow and droplet motion in wave-plate mist eliminators with drainage channels, Chemical engineering science, 63, (2008) 5639-5652. 18. F.S. Manning, R.E. Thompson, R.E. Thompson, Oilfield Processing of Petroleum: Crude oil, PennWell Books, (1995). 19. H. Lomax, T.H. Pulliam, D.W. Zingg, Fundamentals of Computational Fluid Dynamics, Springer Berlin Heidelberg, (2013). 20. R.W. Pryor, Multiphysics Modeling Using COMSOL?: A First Principles Approach, Jones & Bartlett Learning, (2011). 21. Comsol, Comsol 4.3 CFD Module User's guide, DOI, (2012). 22. W.B.J. Zimmerman, Multiphysics Modelling with Finite Element Methods, London, (2006).

Excerpt from the Proceedings of the 2016 COMSOL Conference in Boston