IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 53 LARGE EDDY SIMULATION FOR TRANSESTERIFICATION OF WASTE VEGETABLE OIL A. K.M. MOHIUDDIN AND NABEEL A. ADEYEMI Department of Mechanical Engineering, Faculty of Engineering, International islamic University Malaysia, Jalan Gombak, 53100, Kuala Lumpur, Malaysia. mohiuddin@iium.edu.my (Received 15 March 2015; accepted 26 March 2015; published on line 29 May 2015) ABSTRACT:This paper provides detailed information involved in the numerical simulation of transesterification of waste vegetable oil (WVO). The main objective of this work is to perform a mixing study based on large eddy simulation particle image Velocimetry (LES-PIV) which is resolved in the turbulent scale. Reynolds stress model (RSM) was subsequently used to validate the result using a multiple reference frame (MRF) approach for the impeller-vessel geometry. Experimental FAME yield and liquid velocities were found to be dependent on stirrer speeds, impeller bottom distance and bulk flow pattern. Thermodynamic properties of the reaction components were incorporated as user defined function (UDF) for the mixing models. FAME yield was predicted in terms of species concentration and compared fairly well with experimental condition for 1-L and 2-L STR, where yield from the numerical model varied by about 18 and 23 % for 1-L and 2-L STR respectively. The characteristic time scales were used to show the relevant mixing scale to describe the process. ABSTRAK: Kertas kerja ini memberikan informasi terperinci dalam pentransesteran simulasi berangka sisa minyak sayuran (WVO). Matlamat utama penyelidikan ini adalah untuk menjalankan kajian bauran berdasarkan kepada Velosimetri imej partikel simulasi pusaran besar (LES-PIV) yang terlerai pada skala turbulen. Model stress Reynolds (RSM) kemudiannya digunakan bagi mengesahkan keputusan menggunakan kaedah kerangka rujukan berbilang (MRF) untuk geometri bekas-pendesak. Keputusan eksperimental FAME dan halaju bendalir didapati bergantung kepada kelajuan pengacau, jarak dasar pendesak dan corak aliran pukal. Ciri-ciri termodinamik komponen reaksi digunakan sebagai fungsi mentakrif pengguna (UDF) untuk model bauran. Hasil FAME dijangkakan dalam terma konsentrasi spesis dan dibandingkan dengan adilnya dengan kondisi eksperimental untuk STR bagi 1-L dan 2-L, dimana hasil daripada model berangka didapati berbeza di antara 18 dan 23% untuk STR bagi 1-L dan 2-L. Untuk menggambarkan proses, skala berciri masa digunakan bagi menunjukkan kaitan hubungan skala bauran. KEYWORDS: transesterification; CFD; waste vegetable oil; STR; impellers 1. INTRODUCTION The main objective of this paper is to use computational fluid dynamics (CFD) modeling approach to integrate reaction measurements and computational modeling in a multi-scale framework. The CFD model of the WCO transesterification is proposed to be used in the development process condition with few experimental data regarding transesterification of WVO. IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 54 The production of biodiesel in STR takes place in the turbulent range. Understanding the mixing regime for process improvement is quite important. Literatures on biodiesel production mostly report radial impellers to achieve mixing during oil transesterification [1]. At a stoichiometric molar ratio of 3:1 (alcohol to oil), the conversion of oil to biodiesel takes place in the presence of a homogenous catalyst where the reaction is understood to be either the formation of two intermediates before forming the final product [2] or the formation of only one intermediate subsequent to formation of the final product [3]. Formations of undesirable products are often mitigated by increasing rate of reaction with excess alcohol [4]. However, studies on effects of impellers are few. Different mixing impellers produce certain different effects in multi-component reactions [5]. Where it is necessary to promote yield of a single component without unnecessary side products, mixing velocities is selected to produce the right mixing-scale and a mass transfer rate. Flow varies significantly both in magnitude and direction within the STR and difficult to experimentally quantify. Where the flow in STR have been quantified and compared reasonably with a numerical model [6], the added complexity due to reaction gives rise to uncertainties which are not addressed in kinetic reaction models. Dynamic modeling couples the physical phenomena (reaction and flow) using numerical methods to study flows and reactions. It provides a broader insight into the equipment designs and compares performance of the reactors under different operating conditions [7]. In principle, numerical solutions of flows, in confined space, with appropriate boundary and initial conditions are used to study mixing processes and their solutions are applicable over a wide range of conditions that are rather expensive to obtain experimentally. 2. WVO TRANSESTERIFICATION PARAMETERIZATION Transesterification was carried out in a 2-L stirred tank reactor (STR) (Fig. 1). The STR dimensions are listed Table 1. The NaOH catalyst (1-1.5 ± 0.02 % oil wt) was dissolved in methanol (140 ml). The oil-to-alcohol stoichiometric ratio was kept at 1:6 and the excess methanol was used to drive the reaction forward and prevent auto-catalysis. The reaction conditions with regards to the molar ratio of oil to alcohol and catalyst amount was according to the optimized value [8] used for fryer grease. All mixing was carried out for 60 min with a variable speed motor, automatic timer and temperature controller. Products samples were rinsed with hot water (90°C) and dried for 6 hr. Fig. 1 2-L laboratory STR. IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 55 Table 1: Description of tank an impeller geometry with dimensions Tank Height, H (mm) 200 Tank diameter, T (mm) 130 Mixed flow impeller, Impeller Diameter, D (mm) Impeller blade height (mm) Impeller blade length (mm) No of blades Blade angle to the shaft (downward/axial flow) 55 D/4 D/1.5 3 50 Impeller bottom distance (IBC), C T/4 T/5 T/6.5 2.1 Evaluation of Free Fatty Acid (FFA) in WCO For the evaluation of the FFA in WCO sample, 50 ml isopropyl alcohol was prepared and added to 1 g of WCO sample in the presence of 1 g of phenolphthalein. The solution was titrated against 0.1 M KOH until color change was noticed in the solution. This was done in triplicates and the average FFA was calculated by Eq. (1). (1) 2.2 Taguchi Method Design of Experiment The experiment was based on the Taguchi Orthogonal Array (OA) L9 method [9]. Three variables (temperature, speed and IBC) at 3 levels were used to design the experiment (Table 2) in MINITAB 14. The Taguchi method based on OA reduces variance for the experiments by optimizing setting of control parameters. To identify the effect of temperature, impeller speed, and IBC on yield, the Signal to Noise ratio (S/N) of the Taguchi method (which are log functions of desired output) was used for data analysis and prediction of optimum parameters. Table 2: Taguchi design variables Variables Levels 1 2 3 Temperature, T (°C) 60 65 70 Impeller speed, S (rpm) 600 650 700 Impeller bottom clearance, IBC (mm) 20 25 30 The S/N is defined as (2) where, n is the number of observations in the sub-sample and y is the data observations in the subset. The S/N is used to calculate the performance statistics values and contribution ratios of each parameter. Thus the combination of design of experiments with optimization of control parameters to obtain the best results is achieved in the Taguchi method. OA provides a set of well balanced (minimum) experiments. Without the Taguchi method, 16 experimental runs would have been needed to investigate the effect of each of the three 25.6 0.0864 NaOH volume % difference in FFA sample weight      2 1 1 10 log Bigger yS N n                IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 56 parameters instead of the 9 used. With the main effect analyses, possible combination of optimum parameters is predicted. These were further analyzed to determine the interactive effect of the three parameters. 2.3 Standards and Sample Preparation for T, DG, MG and Glycerol Analysis GC-MS analysis of biodiesel according to ASTM D6584 was used to verify the G, MG, DG, T, and total glycerin content in WCO and FAME sample. Experimental samples were analyzed on a DB 23 Agilent 6890 GC system with a split/splitless inlet and a column of 60 m × 0.25 mm id × 0.15 μm [10]. Five GC calibration standards were prepared by mixing aliquots of the individual stock standards in proportions specified by the ASTM method. After mixing, 100 μl of the derivatization agent, N-Methyl-N- (trimethylsilyl) trifluoroacetamide (MSTFA) was added to each calibration standard. After 20 min, 8 ml of reagent grade n-heptane was added to each calibration standard. These final reaction mixtures were directly injected into the gas chromatograph. 2.4 Sample Preparation for Fatty Acid Analysis WCO samples were analyzed by weighing 100 mg sample in a 20 ml test tube (with screw cap) and dissolved in 10 mL of hexane. 100 μl of 2 N KOH in methanol (11.2 g in 100 ml) was added and vortex for 30 seconds. The supernatant was then transferred into a 2 ml auto sampler vial for GC analysis. 2.5 GC-MS Analysis of T, DG, MG and Glycerol FAME yield was obtained by Eq.3 during the reactions at time interval of 5, 10, 15, 30, 45 and 60 min. T, DG, MG and Glycerol content in the FAME samples were estimated using the relative retention times from the GC reading [10]. The retention time of the first internal standard, 1,2,4-butanetriol was used to identify and quantify free glycerin (G) while that of the second internal standard, tricaprin was used to identify and quantify the MGs, DGs, and Ts in from the monoolein, diolein, and triolein calibration functions respectively. Total glycerin (TG) was obtained using Eq.4. (3) where, MFAME = mass of ester, MWFAME = molecular mass of ester, MWoil = molecular weight of oil and Moil = mass of oil (4) 2.6 Analyses of Biodiesel Properties The ignition quality tester (IQT) was used to measure the cetane number. Kinematic viscosity and density were measured by using a KVS 702 tensiometer according to ASTM D445; pour point according to ASTM D97 by using a cloud and pour point analyzer (ISL CPP 97-2); flash point by using a Pensky-Martens flash point tester (Petrotest PMA 4) based on ASTM D93; moisture content was determined by the Karl Fisher method; and oxidative stability by Metrohm 743 Rancimat according to EN 14115 method and glyceride content was measured according to the method described above. All tests were carried out at the Malaysian Palm Oil Board, Kajang, Selangor, Malaysia. 2.7 Single Phase Two Dimensional (2-D) Particle Image Velocimetry In order to model the reaction, the hydrodynamic parameter such as velocity of impeller, turbulent kinetic energy and dissipation rate are needed. These were later ( %) 100 3 FAME oil oil FAME M yield mol M MW MW     TG G 0.255M 0.146D 0.103T    IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 57 compared with numerical results obtained from the CFD model and estimated from flow measurement of two mixing impellers using PIV technique. 2-D velocity data of the liquid flow filed in a geometrically similar cylindrical acrylic flat-bottomed tank (Fig. 2a) was set up as (i) unbaffled and (ii) baffled (fitted with 4 vertical baffles equally spaced radially at 90) tank were used along with a Dantec PIV equipment at high spatial resolution, which allowed the estimation of spatial velocity gradients and ε. The laser light sheet thickness was approximately 2 mm (0.15% of tank diameter) and aligned parallel to the z axis, 30 mm from the shaft center, just at the tip of the impeller. When the laser sheet is focused at the cente of the vessel, reflections from the impeller/shaft surface cast a shadow over the half of the PIV in the images as a result of the reflective steel surfaces. The camera was focused on a 150 × 120 mm2 area which covered the whole impeller stream region and nearly to the tank wall in the horizontal. For PIV measurement the set up parameters are listed in Table 3. The vessel is filled with water ( =1000 kgm3, = 10-6 m2/s). (a) (b) Fig. 2 (a) Cylindrical acrylic flat-bottomed tank (b) Mixed flow impeller. Table 3: PIV recording parameters for flow in STR Flow geometry : parallel to light sheet and plate Max in-plane velocity : Umax≈ 1.2m/s Field of view : 150 × 120mm2 Observation distance : Z0 ≈ 0.6m Recording method : dual frame/single exposure Recording medium : full frame interline transfer CCD Illumination : Nd:YAG laser 120mJ/pulse Pulse delay : Δt = 263μs Seeding material : psp (dp ≈ 20 μm), =1030 kg/m3 2.8 Flow Measurement The PIV readings were performed on the plane x-y and y-z to the impeller shaft. To define the flow conditions in the 2 L batch STR, cylindrical acrylic tank (3 mm thick) was placed in a square tank with water to observe distortions in the PIV.   Field of IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 58 Fig. 3 Plane of PIV measurement (a) z-x plane (b) z-y plane. Measurements were taken at impeller speeds of 600, 650 and 700 rpm. In the turbulent flow regime, the instantaneous velocity of the fluid at any point is impossible to predict, hence a statistical approache is used to calculate mean and root mean square (rms) velocities provided the flow is steady on average (constant flow rate). For each set of 2-D measurements, 600 pairs of pictures were taken at a frequency of 4 Hz and the measurements were repeated a number of times to establish reproducibility. Basically, a pair of images is divided into smaller regions (interrogation windows). The cross- correlation between these image sub-regions measures the optic flow (displacement or velocity of the objects) within the image pair. Two consecutive images were cross- correlated to give the velocity field with a total of 39 × 31 velocity vectors for the entire plane. The time separation between two exposures was estimated using Eq.5 (5) where S is scale factor, N is interrogation area, d is pixels pitch, and V is velocity (estimated). In our case N=64 pix, d = 6.7 μm where the time used in the experiment should be less than the calculated t. The maximum desired particle displacement, Dmax, between pulses was ensured to be one quarter of the interrogation length, LIV, where maximum separation was then given as a function of the magnification of the image, M and the maximum velocity, Utip, within the system ௠௔௫ܦ = ଴.ଶହ௅಺ೇெ ௎೟೔೛ (6) where M is image magnification in μm/ pixel. Interrogation area (IA) of 64 × 64 pixels was used for the first pass and 32 × 32 pixels were adopted for the second pass at 50% overlap between IA, and the final resolution of the velocity measurements was 32 × 32 pixels. Post-processing of the PIV data for an interrogation area of 32 × 32 and 64 × 64 pixels were compared to determine the best size. Mean velocity using 64 × 64 was seen to be closer to the theoretical value compared to results using 32 × 32 as shown in Fig. 4. In order to obtain good quality results, the vectors obtained from the PIV were processed by removing the outer rows of vectors. Spurious vectors around the image edges occur due to high pixel sensitivity in that region. Also vectors which have a velocity greater than 1.5 times the tip velocity were taken to be spurious.    t S 0.25N 2 d / Vmax Vmin   IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 59 Fig. 4 Static PIV calibration. Impeller rotation is clockwise as viewed from above and origin of the coordinate system used is the base of the center of the vessel with and being the mean and root-mean-square (RMS) velocities in the axial (z, directed upwards), radial (r, directed outward from the shaft) and tangential ( , directed out of the plane) directions respectively. The tip velocity (Vtip) and Reynolds number (Re) are estimated for the impeller according to Eq.7 and Eq.8 and data from Table 1. (7) (8) and average volume dissipation rate in STR, is (9) and (10) where power, P is the impeller power needed for the mixing, m is the number of blades and T is the torque. Considering that eddies present in the flow must be accounted for, spatial resolution (Δ) must be similar to the Kolmogorov length scale ( ), where the mean Kolmogorov length scale is estimated as: (11) In previous studies, spatial resolution has also been shown to strongly affect the value of dissipation [11]. Spatial resolution is required to be similar to when studying turbulence in order to take into account the smallest turbulent scales in PIV measurement [12]. In order to capture the turbulence scale, the PIV resolution (spatial resolution) must be similar to the Kolmogorov length scale. To resolve the seeding particle relaxation time , Eq.12 is a convenient measure for the tendency of particles to attain velocity equilibrium with the fluid. This Eq.13 is compared with Eq.12 such that time due to u ,u ,uz r  , ,z ru u u    D NtipV    2Re /N D   P V    P 2 NmT   (1/ 4) 3v           p IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 60 diameter of the largest particle is smaller than time associate with the Kolmogorov time scale, Eq.9 is computed from the overall mean dissipation rate at the rotation speed of the impeller [12]. (12) (13) (14) (15) (16) (17) (18) Using the analogy of LES, a large eddy PIV method is used for estimation of the dissipation rate [13]. Here, the biggest eddies in high Reynolds number turbulent flows absorb energy from the mean flow, and transport it to smaller eddies. The large eddies are highly anisotropic, and strongly flow-dependent whereas the smallest eddies involved in TKE dissipation tend to be flow-independent. This is why the spatial resolution of PIV exceeds the smallest eddy sizes if dissipation rate is to be measured. Hence the LES analogy can resolve both velocity field (analogous to solving the NS equations), and model the unresolved scales using an interrogation cell size that is naturally the spatial filter. Going by the LES analogy, an inertial sub-range exists between the integral and Kolmogorov scales. It is assumed that for dynamic-equilibrium that the flux of TKE through the inertial sub-range is equal to the dissipation rate. To estimate this, only length scales within the inertial sub-range are required and it is not necessary to resolve the velocity field down to the Kolmogorov scale. For the mixed flow impeller, the mass flow rate, Q (m3 /s) through the impeller is taken as mass flow out of the impeller both at the side and the bottom. Radial mass flow rate, Q (m3 /s), was calculated by (19) over a radial distance equal to the radius of the impeller and integrated at the height, Z, of the impeller. Pumping number and power number are estimated using Eq.20 and Eq.21 respectively. (20) and p   2 p p p 18 d v     1/ 2 v         Q 3 Q N = ND b c b bN w DFl a b D T                0 3 2 IR z QR v rdr N ND    0 3 2 RR z QS v rdr N ND      i ii i p a r   2 1 Q 2 . z r rz r U dz  2 1 32 / z Q z N rU ND  IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 61 (21) To characterise the performance of impellers, non-dimensional power number of Np (Eq.21) is calculated. Impeller speed 600, 650 and 700 rpm corresponds to Re = 3.16 × 104, 3.43 × 104, and 3.69 × 104 with μ=0.789e-3 Ns/m2 at 25C for water. 2.9 Pre-processing of Computational Domain In the pre-processing stage, the tank and impeller were drawn in a CAD environment and imported into GAMBIT as a para-solid file for the meshing process. This ensured that the profile of the vessel is preserved when exported into Gambit. Even though other options like STP, IGES or STL file format can be used to transfer the geometry into neutral format, loss of geometry is often experienced. Edges that needed to be reconnected in GAMBIT when IGES format is used require time consuming refinement for the grid process. Although GAMBIT is not user friendly, as a pre-processor tool for geometry modeling, generated mesh is transferable to FLUENT without loss of features. 2.10 Meshing The computational grid was divided into three parts; vessel, rotation part and impeller geometry. A hybrid mesh over the full three dimensional domain was adopted. The multiple reference frame (MRF) approach was selected for handling the stationary and the rotating walls in a fully-predictive way (Fig. 5 and Fig. 6). The vertical size of the rotating reference frame was 60 mm in diameter and 55 mm in height. The solution convergence was monitored as the residual of the variable approaches to 10−7. Mesh resolution was obtained with 3 mesh sizing using GAMBIT. The full vessel was modeled as the result of the shape of the impeller where blade overlaps into the other sector of the next blade. Table 4 shows the mesh employed in this study. Unstructured hexahedral cells ensured the edge of the geometry were not lost during geometry translation. Skewness of all the meshes was below 0.81 as required for ANSYS FLUENT. (a) (b) Fig. 5 (a) Plane cut (b) 3-D unbaffled STR with mixed-flow impeller using the MRF. 3 5 2 p NN N D     IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 62 (a) (b) Fig. 6 (a) Plane cut (b) 3-D baffled STR with mixed-flow impeller using MRF. Table 4: Total cell number in baffled and unbaffled vessel for mixed flow impeller with multiple reference frame (MRF) baffle vessel unbaffled vessel MRF tank MRF Tank Mesh 1 64437 98353 88043 69591 Mesh 2 187601 98353 82973 214985 Mesh 3 110396 407917 85912 321827 Mesh 4 187601 645504 761130 10171 Mesh 5 560013 1225747 934194 90345 The MRF boundary was located at r/R = 0.46. The vessel wall was specified as stationary in the moving zone. For the MRF approach, no averaging was carried out while exchanging the information between the rotating and stationary zone in a fully-predictive way. The distance between the impeller blades and baffles were assumed to be sufficient so that the flow around the impeller is unaffected by the rest of the tank using the MRF. Wall functions were used to specify wall boundary conditions. The top surface of the liquid pool was assumed to be flat and was modeled as zero normal velocity with zero- shear. The tank wall, shaft and impeller had a no-slip boundary condition. The torque from the numerical simulation is plotted against the different grid numbers (Fig. 7) as a preliminary validation of the mesh size. From the grid dependency study, calculated torque from the numerical simulation was varied by 41% for mesh 1 (162790 cells), 23% for mesh 3 (518313 cells), 21% for mesh 4 (833105 cells) and 24% for mesh 5 (1785760 cells) from the experimental torque estimated as an average value of 0.01 Nm. Fig. 7 Optimization check using five different mesh sizes. M IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 63 The difference between the torque estimation of 800,000 cells and the experimental value was considered to be within an allowable estimate and thus was used for this study. Further increment in cell numbers did not improve calculated torque. 2.11 Convergence Monitoring Scaled residuals were monitored to ascertain convergence until the values of residuals were reduced to three to five orders of magnitude. The momentum, lift and drag were also monitored thru the impeller surface. The integral of static pressure on the faces of the impeller blade was also included. 2.12 CFD Mixing Model The mixing model presented here follows the work of Baldgya [6] using a closure modeling approach of transport equations. The average reaction rate is calculated based on a model of the local mixture structure accounting for concentration fluctuations on different scales [14]. Vicum et al. [15] employed a multi-scale model, where a mean and fluctuating concentration of scalar tracer captured the different scales [16]. To consider the macro-mixing, the concentration of a scalar trace f from a mixture fracture is taken as (22) where Dm is the molecular diffusivity and DT the turbulent diffusivity. To consider the micro-mixing, the variance of the mean of a tracer concentration and local concentration of f is decomposed to . The model of the decomposed variance represents the inertial-convective, viscous convective and viscous diffusive mechanism. The time scale for the changes in variance due to the inertial- convective is thus defined as . Similarly dissipation of and production of is and . This way the model accounts for mixing at all the scale. A probability density function , approximated by a beta function is assumed for the mixture fraction. Putting in terms of and f gives; (23) and are local concentration of the WCO and FAME. The concentration of FAME is thus solved by a mass balance approach (24) DR j m T j j j f f p f u D D t x x x                    r   2 2 s f f   2 2 2 1 2 3, ,   2 1 / s  2 2 2 3 2 2 / E  2 3 / G  ( )f ( )f 2i 1 2 2 2 0 ( ) ( ) ( )r k CaCb k Ca f Cc f f df     ( )aC f ( )cC f 2j m T j j j Cc Cc p Cc u D D r t x x x                     r IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 64 However, the function describes the mixture fraction and can be used to determine , and f simultaneously [15]. For further details see Baldyga and Bourne [6] and Baldyga et al. [16]. To solve conservation equations for chemical species, the conservation equations take the form of Eq.16. The mass and transport equations have been studied with various approaches. In this work Reynolds averaging approach is adopted, where the velocity and concentration of the transport equation are decomposed into a fluctuating and average concentration and velocity components: (25) (26) (27) Eq.24 and Eq.25 inserted into Eq.23 yields (28) The terms and give an infinite set of moments because they must be resolved for each set of the species involved in the reaction. Ordinarily, these terms could be ignored but the reactions are at the sub-grid level and these are the terms that explain the phenomenon at this level. First and second [17] order closure models have been used to close the equation. Second order models yield higher order unclosed terms when a relatively larger number of equations have to be solved. When the species are perfectly mixed, the fluctuation term will approach to zero. However, the Magnussen model, developed for combustion is used to incorporate the effect of the reaction rate on the fluctuation and average component implemented in the FLUENT solver as: (29) Two mixing rates from the reactants and products that depend on the local turbulence kinetic energy and dissipation rate are computed as (30) (31) 0 , , a b c b c a b bulk c bulk C C C C C C C C       ( )aC f ( )cC f ', jC U  ,j jC U 2 12 j j j i a b i i C C C U D k C C t x x          j j jC C C  j j jU U U      2 ' ' ' ' 12 j j j i j j a b a b i i i C C C U D U C k C C C C t x x x              ' ' j jU C ' ' a bC C ', _ 1 1 exp j k f k k n nN N a k i k i k i k j i k i j j j E R M A T C K M C RT                              1 R M i k i k i k i k R m R M A M          2 , ' ' ' R P M i k i k i N j k j j m R M AB M            IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 65 The approach accounts for the high turbulence mixing. Since this study is turbulence oriented, the yield in the STR is determined by the degree of turbulence, power requirement, pumping capacity of the impeller and recirculation time through the STRs. Hence, the accuracy of the simulation is evaluated by comparing simulated results with available field measurements. A variety of turbulence models, convection schemes as well as different grid resolutions were employed and the models performance assessed. The major challenge was in the coupling of the reaction model, which contained unresolved component. 2.13 CFD Process parameter The density for the WCO and FAME alcohol is estimated as 848.8 – 907.9 kg/m3 and 857 – 885 kg/m3 respectively, which has been verified experimentally. According to Kay’s rule [18], thermo-physical properties of FAME are estimated from their methyl ester composition. Using the linear regression of the molecular weight of FAME to specific heat [19], the specific heat for the FAME and WCO can be obtained as; (32) Specific heat for T, DG was calculated as (33) where is the ideal gas specific heat capacity from the Joback group contribution method and obtained as. (34) are parameters obtained based on contribution of methyl, alkyl, ester groups present in methyl ester. Better estimation has been reported with modified Rowlinson correlation with Constantinou and Gani acentric factor correlation, [19]. R is the universal gas constant. Specific heat for MG and Glycerol is provided by the expression in Eq.35 and Eq.36 respectively [20]; (35) (36) The viscosity of FAME as a function of molecular weight and number of bonds is expressed as; (37) 21 06( ) 0.0009( ) 1.9291pC E M M     0 1/3 1 1 1.45 .45 / (1 ) 0.25 [17.11 25.2(1 ) 1.742(1 )] p p r r r r C C T T T R T             0 pC 0 2 3 37.93 0.210 0.000391 2.06 7 p j j j j j j j j C n a n b T n c T n d e T                                        , , and a b c d     2955.54 1703.1 4.1768pC T T   78, 468 480.714pC T  ( ) 12.503 2.496. ( ) 0.178i iLn In M N     IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 66 where M is the molecular weight and N is the number of double bond [21]. With the strong dependency of viscosity on temperature [20], dynamic viscosity for T and MG is expressed as Eq.38and Eq.39, respectively; (38) (39) Glycerol viscosity is taken as 3.5-5 mm2/s according to EN14214 standard. The thermal conductivity, which is linearly dependent on temperature, is expressed as Eq.40 and Eq.41 for MG and glycerol respectively (40) (41) 2.14 Time Scaling The time-scale of reaction are defined in Eq.42, Eq.43 and Eq.44 for, macro, micro and meso-mixing respectively; (42) where is the pumping volume around the impeller. This is obtained as the integral of the area over the impeller using with a three-dimensional CFD model. is a constant value. Eq.48 accounts for the diffusivity at the feed point of mixing. However, the feed point was not included in this case study as the NaOH/alcohol fraction was assumed to be part of the mixture. (43) The micro-scale mixing is reported to be driven by viscous-convective deformation [6], which has a time characteristic expressed as Eq.44. (44) The molecular diffusion that produces micro-mixing at this level where the eddies, r is larger the Kolmogorov mixing is accounted for by Eq.45. (45) v is the kinematic viscosity, Dm is the turbulent diffusivity and 218584exp 136.48 18.61ln 2596500T T T         1731.2 exp 11.67 T        0.22919 0.00019T   0.258 0.0001134T   C feed V Q   3 1 impcQ C Nd 1C feed D T Q uD   1/ 2 1 17.24E v E           1 17050 0.030G E Sc           m v Sc D  IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 67 The combination of all the above time scales is almost impossible to be included in a single model. The multi-scale time-scale model was used to identify the relevant mechanism to be considered in modeling a bimolecular reaction [14]. It was concluded that a time-scale based on micro-mixing did not affect the yield prediction; hence it was safely ignored. Although, that study was a parallel reaction involving bimolecular species with feed point, our current study is limited to the single step transesterification without consideration of a feed point. Thus, time scale due to feed addition time has is neglected. The comparison of the time scale forms the basis of characterization and subsequent CFD modeling of the reaction with focuses on the micro and macro-mixing effect. The characteristic time for the reaction is taken as Eq. (46), with mixing time-scales calculated as average values. (46) The present CFD mode employs the cell-centered finite-volume method that allows use of computational elements with arbitrary polyhedral shape method [22]. The convective terms were solved by second order accurate upwind scheme and the diffusion terms as second order accurate central-differencing scheme. The velocity-pressure coupling was achieved using the SIMPLEC type segregated algorithm for instructed, non- uniform mesh [23]. Simulation of the full 3-D geometry of the STR was implemented in ANSYS Fluent. For the cases conducted, convergence was set at 0.5  10-5 for the mass residual. These calculations required 42-76 hours of CPU time using a 3 GHz Intel processor and 4GB RAM. 3. DISCUSSION AND ANALYSIS 3.1 PIV: Mixed Flow Impeller Figure 8a and 8b show the PIV result for the mixed axial flow in unbaffled and baffled vessel. The maximum radial velocity, Ur/Utip is 0.1 and 0.38 in unbaffled and baffled vessel respectively. The flow feature in the unbaffled vessel presents a streamline away from the impeller face with a distinct up-pumping flow field. In the baffled vessel, the mixed flow feature was however more obvious with the recirculation loop/ vortexes appearing near the wall at the baffle. Axial flow was also observed to be parallel to the shaft, towards the center of the tank. The velocity due to the impeller at the bottom of the tank did not change even with increase in impeller distance (not shown). Consequently, it is assumed that corresponding energy will not vary significantly with the change in impeller distance from the bottom of the vessel. It was also noticed that axial and radial flow were combined with this impeller. This flow is similar to impellers with strong primary and secondary flow [24]. In this PIV study, this feature was noticed to mask velocity disparity in the unbaffled tank (Fig. 8a). Due to the combined axial and radial flow for the baffled vessel, near the baffle and shaft respectively (Fig. 8b), the result was referenced against the data for the axial impeller. It showed consistent similarity to flow generated by an axial impeller [24]. 1 2 1 ( )Ri ik c c    IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 68 (a) (b) Fig. 8 PIV Velocity magnitude at impeller speed of 600 rpm in (a) unbaffled and (b) baffled vessel. 3.2 LES analogy Using LES analogy the dissipation rate was determined from the TKE balance [13]. The LES analogy accounts for flow and resolves turbulence at the smallest scales using a sub-grids scale model. The spatial distributions of ε (0.4N3D2) from this method are shown. The single phase flow generated by the mixed flow impeller is an STR. Fig. 9 Dissipation rate using TKE from pseudo-isotropic assumption, assuming 2-D PIV ensemble-averaged measurements at 600 rpm. Fig. 10 Field of view for mixed flow impeller. 5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 0.5 1 1.5 2 2.5 3 FOV IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 69 In order to determine whether the magnitude of the estimated dissipation rate were reliable, the power dissipated in the field of view (FOV) was compared to the total power input by the impeller. Power no. for the axial impeller is substituted for the power no. of mixed flow impeller [25]. Power input, was calculated to be 0.4983 W. The power dissipated in the field of view, simplified for 2-D application was obtained by double integration of the area over the field of view; (47) The power dissipated in this field of view is calculated as 0.5344 kW. was 0.04. The ratio of the two parameters is generally agreed to be < 1 [26]. However, when the dissipation is known to be non-homogenous, variation could be significant. For comparison purposes, the dimensional analysis and LES estimation methods yielded ratios of 0.04 and 0.05 respectively. Based on the average dissipation rate from the power input (ε) (Table 5) the Kolmogorov length scale was found to be 31.88 μm at 600 rpm which corresponds to ~35 PIV resolution. This should have been ~2 in order to capture 90 % of the dissipation rate [27]. However, this requirement for research is difficult to resolve, probably due to the high solidity ratio and high magnitude in the discharge stream. The ε estimated here was subsequently used for the CFD simulation for the mixed flow. 3.3 Experimental Validation To validate the experimental result, power drawn for the impeller was calculated from current consumed during stirring [28]. Similarly, the power was used to estimate the torque required to rotate the impeller. The impeller power number, compared reasonably with the theoretical (Table 5) between 600 to 700 rpm (i.e. 10 – 11.7 rps). Table 5: Hydrodynamic properties of mixed flow impeller N (rps) D (m) ReL x 10 3 Power (kW) ε (W/kg3) η (μm) Tk (μs) 10 0.055 3.03 - 0.34 0.44 0.96 31.88 5.16 Mixed 10.8 0.055 3.27 - 0.32 0.46 1.22 30.09 4.10 flow 11.7 0.055 3.57 - 0.34 0.44 1.55 28.33 3.23 No technical information was available for the mixed-flow impeller used. Thus the Reynolds mixing number was read from a power curve. This approach was similarly used to arrive at the Reynolds number, Npc for the axial impeller [29]. The power obtained from the electric current measured during unloads and loading at 100-700 rpm is shown in Fig. 11, corresponding to a specific power (P/V) of 82.07 W/m3 for the mixed flow impeller respectively. The net volumetric power being higher in the mixed flow impeller increases the amount of material transported through the vessel. Power consumption increased steadily with speed for both impellers until 400 rpm. At speed higher than 400 rpm (i.e. turbulent region), this change becomes less significant [30]. INP 2 2 1 1 2 z r FOV z r P rdrdz    FOVP /FOV INP P /FOV INP P /FOV INP P pN T pN c pN T pN c IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 70 Fig. 11 Power drawn measurements for mixed impeller between 100-700 rpm in baffled and unbaffled vessels. 4. NUMERICAL SIMULATION The numerical simulation carried out in 3-D using the standard к-ε, sst (к-ω) and RSM turbulent models is presented. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm for velocity-pressure coupling and second-order upwind discretization scheme was adopted in obtaining the mean velocity at mass residuals 0.510–5. This strategy did not produce satisfactory result by the models for the baffled vessel. Sahu et al. [31] had earlier reported that simulation at large grid size could not converge easily using the SIMPLE algorithm. However, the problem was subsequently solved by initiating the model as a lamina flow using the first-order scheme and SIMPLEC algorithm to obtain convergence. This was the only algorithm that led to convergence for the baffled vessel. Similarly, none of the models used in this work were successfully solved with the SIMPLE algorithm using unstructured mesh. Several authors [32-34] had noted that discretization scheme had no effect on predicted mean velocity in unbaffled tanks. However, the turbulence prediction is seriously affected. 4.1 Comparison of PIV-CFD Velocity Data The comparison of the values of the velocity components by the к-ε, sst (к-ω) and RSM and PIV data over the field view at x/R = 0.3 (bottom) and 0.46 (top) in an r-z plane are shown in Figs. 12–14 for the mixed flow impeller in baffled vessel. For the baffled vessel, normalized mean velocity was better predicted by the RSM between R = -0.02 to 0.02 m at x/R = 0.38. The highest error between PIV data and RSM occurred near the baffled zone with difference of about 60 %. It was observed that a higher deviation occur on the opposite side of the vessel. Meanwhile, the к-ε and sst (к-ω) model grossly under-predicted the mean velocity at x/R = 0.3 between R = -0.06 m and vessel center. At x/R = 0.46, maximum difference was 102 % and 82.5 % for the к-ε and sst (к- ω) model respectively at R = -0.06 m to the center of the vessel. In Fig. 15, the normalized axial velocity as a function of the radial position showed that the region occupied the impeller had the highest difference for the к-ε and sst (к-ω) models. At x/R = 0.46, the difference observed between the PIV and RSM model over the axial region at x/R = 0.46 was lower compared to the к-ε and sst (к-ω) model. Under-prediction of the axial velocity similar to that for the mean was recorded at x/R = 0.38 for the к-ε and sst (к- ω). 100 200 300 400 500 600 700 154 156 158 160 162 164 166 168 170 172 Speed, rpm P o w er (W at ts ) mixed flow-baffled mixed flow-unbaffled IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 71 Fig. 12 Normalized mean velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for baffled vessel using mixed flow impeller. Fig. 13 Normalized axial velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for baffled vessel using mixed flow impeller. -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0 0.5 1 1.5 Vessel radius, R(m) N or m al is ed m ea n ve l (U /U tip ) top, x/R = 0.46 rsm rke sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0 0.5 1 Vessel radius, R(m) N or m al is ed m ea n ve l (U /U tip ) bottom, x/R = 0.38 rsm rke sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -1 0 1 Vessel radius, R(m)N or m al is ed a xi al v el , U z/U tip x/R = 0.46 rsm rke sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -2 -1 0 1 Vessel radius, R(m)N or m al is ed a xi al v el , U z/U tip x/R = 0.3 rsm rke sst PIV IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 72 Fig. 14 Normalized radial velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for baffled vessel using mixed flow impeller. The normalized radial velocity at x/R = 0.46 was also observed to be better predicted between the baffle and impeller tip with difference value of 17-72 % for the RSM and up to 200 % for the к-ε and sst (к-ω) models. At x/R = 0.3, the RSM performed better over the radial distance measured (Fig. 17). Flow away from the impeller, towards the tank surface is generally poorer for the к-ε model as demonstrated so far. However, the area closer to the baffles is well represented by all the models at distances closer to the liquid surface. Even though the RSM gave a better overall results, values close to the mid of the vessel did not tally with experimental result. By this analysis, the RSM model produced a more realistic outcome. The sst (к-ω) have been reported to give reasonable results where other models are known to fail severely [31]. Similarly, the mean axial and radial velocity component in the unbaffled vessel (Figs. 15 – 17) revealed better prediction at x/R = 0.3 with the RSM when compared with the PIV data. Except where a very poor result was obtained for the radial component at x/R = 0.46 (Fig. 17), all other region had errors less than 35-50%. Again, the к-ε model gave greater error (under 120 %). Overall, the RSM gave reasonable velocity values except at the impeller region, which has the rotating reference frame. The simulation could probably be further improved by using a finer grid around the MRF. Further refinement of the grid in the MRF did not improve the results. This is kept as one of the challenge faced during this study. Furthermore, the κ-ε model seems to be an inappropriate model especially with baffled vessel. The magnitude of the difference in velocity component in the unbaffled vessel varied similarly to that of the baffled vessel as explained above for the axial and mean velocity. For the radial velocity component, the κ-ε model had the least difference in magnitude compared to the RSM and sst (κ-ω) models. Earlier, the radial component was not successfully resolved at locations below and above the impeller in the baffled vessel (Fig. 6). However, the radial velocity in the unbaffled vessel (Fig. 5) span of the plane measured. -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -1 0 1 Vessel radius, R(m) N or m al is ed ra di al v el (U r/U tip ) x/R = 0.46 rsm rke sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -0.5 0 0.5 1 Vessel radius, R(m) N or m al is ed ra di al v el (U r/U tip ) x/R = 0.3 rsm rke sst PIV IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 73 Fig. 15 Normalized mean velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for unbaffled vessel using mixed flow impeller. Fig. 16 Normalized axial velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for unbaffled vessel using mixed flow impeller. Fig. 17 Normalized radial velocity at x/R = 0.3 and 0.46 for PIV and CFD turbulent model at 600 rpm for unbaffled vessel using mixed flow impeller. -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0 0.5 1 Vessel radius, R(m) N or m al is ed m ea n ve l (U /U tip ) x/R = 0.46 rke rsm sst PIV -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0 0.5 1 1.5 Vessel radius, R(m) N or m al is ed m ea n ve l (U /U tip ) x/R = 0.3 rke rsm sst PIV -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 -0.4 -0.2 0 0.2 0.4 Vessel radius, R(m) N or m al is ed a xi al v el (U z/U tip ) x/R = 0.46 rsm rke sst PIV -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 -0.4 -0.2 0 0.2 0.4 Vessel radius, R(m) N or m al is ed a xi al v el (U z/U tip ) x/R = 0.3 rsm rke sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -0.4 -0.2 0 0.2 0.4 Vessel radius, R(m) N or m al is ed ra di al v el (U r/U tip ) x/R = 0.46 rke rsm sst PIV -0.06 -0.04 -0.02 0 0.02 0.04 0.06 -0.05 0 0.05 0.1 0.15 Vessel radius, R(m) N or m al is ed ra di al v el (U r/U tip ) x/R = 0.3 rke rsm sst PIV IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 74 Table 6: Magnitude of difference between CFD and PIV result over vessel radius, R = -0.06 m to 0.065 m for baffled vessel Mean velocity R = -0.06 m to 0 m R = 0 m to 0.06 m κ-ε 7-107.1 17.6-70.6 x/R = 0.38 sst (κ-ω) 5.5-210 33-63.5 RSM 3 - 60.1 0-51 κ-ε 17.2 - 102.5 4 - 200 x/R = 0.46 sst (κ-ω) 10.8 – 82.5 30 – 110 RSM 12.4 - 127 1.2 – 95.1 Axial velocity κ-ε 18-226 ~ x/R = 0.38 sst (κ-ω) 1.1-37.1 ~ RSM 100-112 ~ κ-ε 60-110 ~ x/R = 0.46 sst (κ-ω) 1.7-101.1 ~ RSM 1-23 ~ Radial velocity κ-ε 36-144 3.8-146 x/R = 0.38 sst (κ-ω) 10-125 50-93 RSM 1-17 17-72 κ-ε 91.5-181.1 99-200 x/R = 0.46 sst (κ-ω) 2.1-94.9 1-200 RSM 6-76 5-181 Table 7: Magnitude of difference between CFD and PIV result over vessel radius, R = -0.06 m to 0.065 m for unbaffled vessel Mean velocity κ-ε 11-97 - x/R = 0.38 (sst) (κ-ω) 46-131 73-99 RSM 0-25 - κ-ε 46-200 98-200 x/R = 0.46 sst (κ-ω) 25-1000 65-600 RSM 16-250 61-700 Axial velocity κ-ε 385-100 57-309 x/R = 0.38 sst (κ-ω) 156-2000 10-350 RSM 115-135 13-350 κ-ε 79-140 ~ x/R = 0.46 sst (κ-ω) 56-475 ~ RSM 8-990 ~ radial velocity κ-ε 34-87 36-99 x/R = 0.38 sst (κ-ω) 65-155 ~ RSM 12-97 27-169 κ-ε ~ ~ x/R = 0.46 sst (κ-ω) ~ ~ RSM ~ ~ IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 75 Tables 6 and 7 summarize the magnitude of difference between mean, radial and axial velocity in CFD models and PIV in STR (vessel radius, R = - 0.06 m to 0.065 m) for baffled and unbaffled vessel. In summary, the RSM resolves the mean and axial component better than the κ-ε and sst (κ-ω) models away from the wall and below the impeller until the tip of the impeller region R = -0.025 m. The significance is that the mean/ axial component of the velocity can be estimated accurately going by the comparison with PIV particularly below the impeller. From the estimated values which coincide with PIV data, it could also be assumed that the model has been verified. However, as the assumptions behind the models are time-average based, a model including the fluctuating component would improve the results considerably. 4.2 Comparison of the Mass Flow Rate, Torque, Flow Number and Power Number The comparison of the power number (Np), pumping number (Nq) and mass flow rate (Qr) is shown in Table 8. Since the flow field by the mixed flow impeller was taken as that produced by the axial impeller, the Nq value was expected to fall within the range of 0.68 – 0.86 [35]. The к-ε model did not predict the Nq adequately, with a value of 0.19. While the κ-ω (sst) and RSM gave 0.53 and 0.58 respectively. Np was obtained within the expected range of 0.9 – 1.62 with both κ-ω (sst) and RSM at 0.99 and 1.44 respectively. Np for the mixed flow impeller was over-estimated by the κ-ε model in both vessels. Table 8: Comparison of the mass flow rate, torque, flow number and power number by κ-ε, sst (κ-ω) and RSM models at 600 rpm axial (m3/s) radial (m3/s) Qr (m3/s) torque (N.m) Nq Np (CFD) κ-ε 0.45 -0.13 0.32 0.031 0.19 3.75 mixed flow sst (κ-ω) 0.75 0.13 0.89 0.015 0.53 1.84 (unbaffled) RSM 1.00 -0.036 0.97 0.008 0.58 1.00 κ-ε 0.27 -0.0064 0.26 0.0658 0.16 8.23 mixed flow sst (κ-ω) 1.00 -0.17 0.81 0.0079 0.49 0.99 (baffled) RSM 0.98 0.40 1.10 0.012 0.72 1.40 4.3 Time Scale of Mixing The actual experiment was conducted by adding methanol and NaOH into the STR containing the WCO. The reaction in the 2-L was extended to 60 min. The mixing time was taken as the concentration of the local species of interest instead of inert tracer due to difficulty in measurement; hence the species of interest were directly measured. The mixture model assumed that all the species were contained in the STR at the commencement of mix. This eliminates the need to consider the effect of meso-mixing in the time scale. The mixing-time scales are compared in Table 9. Table 9: Time-scales of the mixing WCO transesterification in batch STR Macro-mixing time Micro-mixing time Circulation time Engulfment Molecular diffusion N(rpm) (s) (s) (s) 600 0.401 0.603 1.705 1-L 650 0.370 0.550 1.421 700 0.343 0.492 1.137 C E G IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 76 The value of being lesser than and . As explained by Vicum et al. [15], when >> 1 molecular diffusion are negligible and micro-mixing is by diffusion especially when Sc are large in the range of thousand. Were this analysis based solely on comparison of and , this would have led to the conclusion that macro-mixing control the reaction as is smaller than and . This observation suggests that model may neglect micro-mixing, however in order to capture the dynamics of the process all the scales needs to be consider. 4.4 Simulation Reaction and Computational Details User-defined function (UDF) and user defined scalar (UDS) were written in C language for ANSYS Fluent to define the thermo-physical and chemical properties of WCO, FAME and the individual glycerides derived from Eqs. 38–47. Based on these considerations, the batch transesterification in STR was modeled using the RSM model. At first and independent flow field was calculated without inclusion of the energy or transportation of the species. This is only calculated once, to obtain the distributions of the average velocity, the average kinetic energy, κ, and the average rate of the energy dissipation, ε. Occasioned by the difficulty in obtaining convergence with multiple reaction steps the FAME – WCO was made the focus of this analysis. The CFD simulation at 600 and 60 C was carried out for 360 time steps. Each time step represented 5 sec of the physical period for a total simulation of 30 min. FAME yield in terms of mass fraction from the CFD simulation result were extracted for a total of 10 intervals. Two monitoring points were used to take record the concentration of FAME during the simulation and compared with experimental data. These two points were kept close below the impeller. Convergence of the solution was enabled by patching an initial concentration for the methanol, WCO and FAME to the fluid domain at an initial concentration of 31% and 16% of the triglyceride. As in the previous experiments, the FAME yield takes an exponential form (Eq.48): (48) Fig. 18 1-L Experimental and simulation result of FAME yield at 600 rpm, 60C in 1-L STR. C E G/E G  G C E C E G 2 2( ) ( ) 1 2 max 2 1 1 k t k tk exp k exp FAME X k k         0 5 10 15 20 25 30 0 20 40 60 80 100 120 140 time (min) FA M E c on c. (d m ol /d m 3 ) Exp 60°C, 600 rpm Eddy dissipation model finite rate model IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 77 A significant agreement between the model and experiment is shown in the 1-L and 2-L cases. Experimental yield data range from XFAME =0.56 in the 1-L STR to XFAME = 0.52 in the 2-L STR on the average. The conversion, XFAME in both cases were > 0.5, indicating that kinetic does not play a role in the mixing and confirming the indication of mixing effectiveness. Similarly, the finite rate model, which also employs the Arrhenius rate in place of the ε/κ ratio as well, under predicts the FAME in both unbaffled STR. The yield trend was similar showing that the exponential form of yield can be used to model transesterification. The results shown is limited to FAME using the mass fraction from the CFD simulation and compared with the experimental results in Figs. 18–19. Fig. 19 2-L Experimental result of FAME yield at 600 rpm, 60C in 2-L STR. 5. CONCLUSION A model for the transesterification of WCO in an STR is presented using a mixed flow impeller. In this work, the hydrodyamic parameter from a single-phase flow field based on the turbulence model was obtained in the rotating reference frame. The process parameters and proporties of WCO and FAME from experiemental work and published literature have been sourced for two mixing models based on (i) kinetic-controlled flow and (ii) turbulence-governed reaction have been compared. The time scale of the reaction showed the kinetics-independence of the reaction, while a reasonable estimate of FAME yield for the 1-L and 2-L reaction was given by the turbulence based model. It should be pointed out that the validation of the stirred flow is based on the RSM model. The 3-D model demonstrates the capabilty to obtain near-experimental result within the set limit, indicating that it is a fully reliable CFD model, not far away for transesterification of WCO. REFERENCES [1] Stamenkovic O, Lazic M, Todorovic Z, Veljkovic V, Skala D. (2007) The effect of agitation intensity on alkali-catalyzed methanolysis of sunflower oil. Bioresource Technology, 98(14):2688-2699. [2] Meher L, Dharmagadda V, Naik S. (2006) Optimization of alkali-catalyzed transesterification of Pongamia pinnata oil for production of biodiesel. Bioresource Technology, 97(12):1392-1397. [3] Freedman B, Butterfield R, Pryde E. (1986) Transesterification kinetics of soybean oil. J. Am. Oil Chemists' Soc., 63(10):1375-1380. [4] Wahlen B, Barney B, Seefeldt L. (2008) Synthesis of biodiesel from mixed feedstocks and longer chain alcohols using an acid-catalyzed method. Energy & Fuels, 22(6):4223-4228. 0 10 20 30 40 50 60 0 20 40 60 80 100 120 140 160 180 time(min) FA M E c on c. ( d m ol /d m 3 ) Exp 60°C, 600 rpm Eddy dissipation model finite rate model IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 78 [5] Yang Y, Yoon RH, Telionis DP, Weber A, Foreman D. (2009) Flow property measurements of stirred-tank flow across three Reynolds number decades. ASME Conference Proceedings, 329-336. [6] Baldyga J, Bourne J R. (1999) Turbulent mixing and chemical reactions. Wiley, New York. [7] Tu J, Yeoh GH, Liu C. (2008) Computational fluid dynamics: A practical approach. Butterworth Heinemann. [8] Issariyakul T, Kulkarni MG, Dalai AK, Bakhshi NN. (2007) Production of biodiesel from waste fryer grease using mixed methanol/ethanol system. Fuel Processing Technology, 88(5):429-436. [9] Rao RS, Kumar CG, Prakasham RS, Hobbs PJ. (2008) The Taguchi methodology as a statistical tool for biotechnological applications: A critical appraisal. Biotechnology J., 3(4):510-523. [10] McCurry JD, Wang CX. (2007) Analysis of glycerin and glycerides in biodiesel (B100) using ASTM D6584 and EN14105. Application-HPI/Petrochemicals/Polymers. [11] Baldi S, Yianneskis M. (2004) On the quantification of energy dissipation in the impeller stream of a stirred vessel from fluctuating velocity gradient measurements. Chemical Engineering Science, 59(13):2659-2671. [12] Delafosse A, Collignon ML, Crine M, Toye D. (2011) Estimation of the turbulent kinetic energy dissipation rate from 2D-PIV measurements in a vessel stirred by an axial Mixel TTP impeller. Chemical Engineering Science, 66(8):1728-1737. [13] Sheng, J., Meng, H., Fox, R. (2000) A large eddy PIV method for turbulence dissipation rate estimation. Chemical Engineering Science, 55(20):4423-4434. [14] Baldyga J, Bourne JR. (1990) The effect of micromixing on parallel reactions. Chemical Engineering Science, 45(4):907-916. [15] Vicum L, Ottiger S, Mazzotti M, Makowski L, Bałdyga J. (2004) Multi-scale modeling of a reactive mixing process in a semibatch stirred tank. Chemical Engineering Science, 59(8- 9):1767-1781. [16] Baldyga J, Henczka M, Makowski L. (2001) Effects of mixing on parallel chemical reactions in a continuous-flow stirred-tank reactor. Chemical Engineering Research and Design, 79(8):895-900. [17] Fox RO. (2003) Computational models for turbulent reacting flows. Cambridge Univ. Press. [18] Kay WB. (1936) Density of hydrocarbon gases and vapors. Industrial and Engineering Chemistry Research, 28:1014-1019. [19] Anand K, Sharma RP, Mehta PS. (2011) A comprehensive approach for estimating thermo- physical properties of biodiesel fuels. Applied Thermal Engineering, 31:235-242. [20] Wang WC, Natelson RH, Stikeleather LF, Roberts WL. (2012) CFD Simulation of transient satge of continuous countercurrent hydrolysis of canola oil. Computers and Chemical Engineering, 43:108-119. [21] Ramírez-Verduzco LF, Rodríguez-Rodríguez JE, Jaramillo-Jacob ADR. (2012) Predicting cetane number, kinematic viscosity, density and higher heating value of biodiesel from its fatty acid methyl ester composition. Fuel, 91(1):102-111. [22] Rhee SH, Joshi S. (2005) Computational validation for flow around a marine propeller using unstructured mesh based navier-stokes solver. JSME International Journal Series B, 48(3):562-570. [23] Versteeg HK, Malalasekera W. (2007) An introduction to computational fluid dynamics: The finite volume method. Prentice Hall. [24] Aubin J, Mavros P, Fletcher D, Bertrand J, Xuereb C. (2001) Effect of axial agitator configuration (up-pumping, down-pumping, reverse rotation) on flow patterns generated in stirred vessels. Chemical Engineering Research and Design, 79:845-856. [25] Aubin J, Le Sauze N, Bertrand J, Fletcher DF, Xuereb C. (2004) PIV measurements of flow in an aerated tank stirred by a down-and an up-pumping axial flow impeller. Experimental Thermal and Fluid Science, 28(5):447-456. [26] Unadkat H. (2009) Investigation of turbulence modulation in solid-liquid suspensions using FPIV and micromixing experiments. Loughborough University, London, Loughborough. IIUM Engineering Journal, Vol. 16, No. 1, 2015 Mohiuddin and Adeyemi 79 [27] Saarenrinne P, Piirto M, Eloranta H. (2001) Experiences of turbulence measurement with PIV. Measurement Science and Technology, 30(1):27-35. [28] Ascanio G, Castro B, Galindo E. (2004) Measurement of power consumption in stirred vessels - A review. Chemical Engineering Research and Design, 82(9):1282-1290. [29] Coroneo M, Montante G, Paglianti A, Magelli F. (2010) CFD prediction of fluid flow and mixing in stirred tanks: Numerical issues about the RANS simulations. Computers & Chemical Engineering, 100(4):277-299. [30] Mununga L, Hourigan K, Thompson M. (2003) Numerical study of the effect of blade size on pumping effectiveness of a paddle impeller in an unbaffled mixing vessel. Third International Conference on CFD in the Minerals and Process Industries CSIRO, Melbourne, Australia. [31] Sahu A, Kumar P, Joshi J. (1998) Simulation of flow in stirred vessel with axial flow impeller: Zonal modeling and optimization of parameters. Industrial & Engineering Chemistry Research, 37(6):2116-2130. [32] Aubin J, Fletcher DF, Xuereb C. (2004) Modeling turbulent flow in stirred tanks with CFD: the influence of the modeling approach, turbulence model and numerical scheme. Experimental Thermal and Fluid Science, 28(5):431-445. [33] Brucato A, Ciofalo M, Grisafi F, Tocco R. (2000) On the simulation of stirred tank reactors via computational fluid dynamics. Chemical Engineering Science, 55(2):291-302. [34] Deglon DA, Meyer CJ. (2006) CFD modelling of stirred tanks: Numerical considerations. Minerals Engineering, 19:1059-1068. [35] Chapple D, Kresta SM, Wall A, Afacan A. (2002) The effect of impeller and tank geometry on power number for a pitched blade turbine. Chemical Engineering Research and Design, 80(4):364-372.