Multiscale Modeling And Simulation Of Macromixing .
Articlepubs.acs.org/IECRCite This: Ind. Eng. Chem. Res. 2018, 57, 5433 5441Multiscale Modeling and Simulation of Macromixing, Micromixing,and Crystal Size Distribution in Radial Mixers/CrystallizersCezar A. da Rosa*,† and Richard D. Braatz‡†School of Chemistry and Food Science, Federal University of Rio Grande, FURG, Rio Grande 96201-900, BrazilMassachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United StatesDownloaded via MASSACHUSETTS INST OF TECHNOLOGY on February 5, 2020 at 04:03:09 (UTC).See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.‡ABSTRACT: Continuous-flow tubular crystallization in whichseed crystals are continuously generated is of interest due to itsenabling of tighter control of crystal properties. This article isthe most detailed simulation study on the design and operationof continuous-flow crystallizers using radial mixers, which havepotential for inducing rapid and intense turbulent mixing andhaving easy construction, high reliability, and low operatingcosts. A multiscale model is employed that couples computational fluid dynamics (CFD), micromixing modeling, energybalance, and population balance equation (PBE) using theopen-source CFD package OpenFOAM. The approach isdemonstrated for the methanol water antisolvent crystallization of lovastatin. A new crystallizer design with multiple radial inlets is proposed and shown to deliver improved mixingcompared to one radial inlet. The effects of varying operating conditions on micromixing and crystal size distribution areanalyzed. A systematic approach is provided for the design of continuous-flow tubular crystallizers with radial mixing.1. INTRODUCTIONCrystallization is widely used in the chemical and pharmaceuticalindustries to perform solid liquid separation and purification toproduce high-value materials such as pharmaceuticals, catalysts,and pigments.1 Several crystallizer designs and crystallizationtechniques have been applied in order to obtain products withhigh purity and desired crystal size distribution (CSD).2,3Among various methods of crystallization used, especially bythe pharmaceutical industry, the mixing of a liquid solutioncontaining the desired solute with a miscible antisolvent toreduce the solubility has the advantage of inducing crystallizationof thermally sensitive pharmaceuticals without large temperaturevariations.1,4 Since the solubility of the solute in the antisolvent isvery low, supersaturation is quickly induced, creating a drivingforce for crystallization. Since this method requires rapid andsufficient mixing of the antisolvent with the solute dissolved insolvent, the design and optimization of such crystallizers play animportant role in achieving crystallization with effective controlof the CSD.5Many different antisolvent mixer designs have been exploredto obtain high supersaturation in order to generate consistentcrystal nuclei that are subsequently grown to a desired size.6 10Over the past decade, state-of-the-art mixers/crystallizers such asimpinging jet and coaxial nozzles have gained more attention.5,11 14 The radial mixer is an alternative that is largelyunknown in crystallization applications but is widely used as athermal mixer in nuclear power plants, chemical plants, andcombustors.15 17 Its success in inducing intense turbulentmixing in other applications suggests that the radial mixer haspotential in antisolvent crystallization due to its easy 2018 American Chemical Societyconstruction, reliability, and low operating costs. However, theuse of the radial mixer in antisolvent crystallization applications,which is referred to here as a radial crystallizer, is limited by a lackof information regarding design and operation for that purpose.Although numerous experimental studies have been carriedout to gain a better understanding of the operation of antisolventcrystallizers,18 24 the number of possible designs and operatingconditions that can be investigated is large. As such, performingbench-scale experiments over the variety and range ofpossibilities can be time-consuming and costly. The applicationof mathematical modeling including transport phenomena tosuch complex systems as radial crystallizers can facilitate thesearch for more efficient processes, to improve the productioncapacity, reduce operating costs, and identify potential operational problems such as fouling on the pipe walls. In addition, thedetailed modeling and simulation of these processes enable theanalysis and a building of understanding of variables that aredifficult to measure experimentally, such as the spatial variationof extent of mixing and nucleation and growth rates.In this regard, Woo et al.13,25 applied an isothermal singlephase model with constant properties coupled with micromixingmodels, the Reynolds-averaged Navier Stokes equation, and aspatially varying population balance equation using a highresolution central difference discretization scheme to simulatethe behavior of batch and impinging jet antisolvent 5433January 23, 2018March 28, 2018April 2, 2018April 2, 2018DOI: 10.1021/acs.iecr.8b00359Ind. Eng. Chem. Res. 2018, 57, 5433 5441
ArticleIndustrial & Engineering Chemistry ResearchPirkle et al.5 extended the model and software of Woo et al.13,25to account for the nonisothermal operation of antisolventcrystallizers and studied the effect of different operatingconditions in the behavior of a coaxial nozzle crystallizer. Suchdetailed investigations have not been published for radialcrystallizers.The objective of this work is to investigate different radialantisolvent crystallizer designs and operating conditions and toimprove the performance of this type of crystallizer viamathematical modeling and numerical simulations. A singlephase model with variable properties coupled with the Fox26micromixing model, a population balance equation using a highresolution finite-volume method, an energy balance, and scalartransport equations was implemented in the open-source CFDpackage OpenFOAM. The methanol water antisolvent crystallization of lovastatin, using kinetics reported in the literature,27was chosen as the model system in the simulations. The influenceof different numbers of radial inlets and operating conditions inthe micromixing, crystal size distribution (CSD), and soluteconversion was investigated.in the CFD grid is divided into Ne different probability modes orenvironments, which correspond to a discretization of thepresumed composition PDF into a finite set of delta (δ)functions:Necontinuity equation:momentum conservation equation: (ρv ⃗) ·(ρvv⃗ ⃗) p ·(τ ) ρg ⃗ t(4)where fϕ is the joint PDF of all scalars, Ns is the total number ofscalars (species), pn is the probability of mode n or volumefraction of environment n, and ⟨ϕα⟩n is the mean composition ofscalar α corresponding to mode n. The weighted concentration isdefined as⟨s⟩n pn ⟨ϕ⟩n(5)The transport of probability and species in inhomogeneous flowsis modeled by p t ⟨vi⟩ i ⟨s⟩n t p p Dt G(p) Gs(p) xi xi xi ⟨vi⟩i (6) ⟨s⟩n ⟨s⟩n Dt xi xi xi Mn(p, ⟨s⟩1 , ., ⟨s⟩Ne ) M ns (p, ⟨s⟩1 , ., ⟨s⟩Ne ) pn S(⟨ϕ⟩n )(7)where G and Mn are the rates of change of p [p1, p2, ., pN] and⟨s⟩n due to micromixing, respectively; Gs and Mns are additionalmicromixing terms to eliminate the spurious dissipation rate inthe mixture-fraction-variance transport equation (see Fox26 fordetails); and S is the chemical source term. The conservation ofprobability requires thatN pn 1(8)n 1andNe Gn(p) 0(1)(9)n 1The mean compositions of the scalars are given byNe⟨ϕ⟩ (2)Ne pn ⟨ϕ⟩n ⟨s⟩nn 1(10)n 1and since the means remain unchanged by micromixing,standard k ε equations: μ (ρk) ·(ρkv ⃗) · μ t ·k Gk ρε Sk tσk Ne Mn(p, ⟨s⟩1 , ., ⟨s⟩N ) 0en 1 μ εε2(ρε) ·(ρεv ⃗) · μ t ·ε C1ε Gk C 2ερ Sε kk tσε μt ρCμα 1n 12. MODEL EQUATIONSThis article employs a multiscale mathematical modelingapproach that couples the dynamic Reynolds-averaged Navier Stokes equations with a multienvironment probabilitydensity (PDF) model26 that captures the micromixing in thesubgrid scale, a population balance equation (PBE) that modelsthe evolution of the crystal size distribution, and the energybalance equation to account for the heat transfer between thesolvent and antisolvent, as well as the heat of mixing andcrystallization.2.1. Conservation of Mass and Momentum Equations.The macromixing was modeled by the Reynolds-averagedNavier Stokes (RANS) model and the standard k ε turbulencemodel with enhanced wall treatment. In order to incorporate theeffect of density difference between the solution and antisolvent,an ideal mixture model was employed to calculate the mixturedensity at every computational grid cell. In general form, theequations are ρ ·(ρv ⃗) 0 tNs pn (x,t ) δ[ψα ⟨ϕα⟩n (x,t )]fϕ (ψ ;x,t ) (11)must be satisfied. In this article, a three-environment model waschosen to account for the micromixing effects. In this approach,the solution of solute/solvent is the environment 1, theantisolvent represents the environment 2, and the mixture ofenvironments 1 and 2 forms the environment 3. According toMarchisio et al.,28 30 the use of three environments is sufficientto capture the micromixing effects in nonpremixed flows withsatisfactory accuracy.Following Fox,26 the micromixing terms for the threeenvironment model are summarized in Table 1, where thevalues of ⟨φ⟩n ⟨s⟩n/pn denote the unweighted variables. Thek2ε(3)The symbols are defined in the Nomenclature list.2.2. Micromixing Model Equations. As in Marchisio etal.,28 30 Woo et al.,13,25 and Pirkle et al.,5 the micromixing effectswere considered by applying the finite-mode PDF modelproposed by Fox.26 In this approach, each computational cell5434DOI: 10.1021/acs.iecr.8b00359Ind. Eng. Chem. Res. 2018, 57, 5433 5441
ArticleIndustrial & Engineering Chemistry Research2.4. Conservation of Energy Equation. In order to applythe energy balance, the three environments are assumed to be inthermal equilibrium at the cell level. This assumption is based onthe time required to achieve thermal equilibrium in a turbulentflow at the cell level, which is, in a worst-case scenario, the sameorder of magnitude as the cell residence time. Also,compressibility effects are neglected since the fluids are in theliquid phase. Thus, the general form of the energy equation canbe written asTable 1. Micromixing TermsmodelvariableG, MnGs, Mnsp1p2p3⟨s⟩3 γp1(1 p1) γp2(1 p2)γ[p1(1 p1) p2(1 p2)]γ[p1(1 p1)⟨φ⟩1 p2(1 p2)⟨φ⟩2]γsp3γsp3 2γsp3 γsp3 (⟨φ⟩1 ⟨φ⟩2)γ γs εξp1(1 p1)(1 ⟨ξ⟩3)2 p2 (1 p2 )⟨ξ⟩32 ⟨ξ⟩3 ⟨ξ⟩32Dt(1 ⟨ξ⟩3)2 ⟨ξ⟩32 xi xi2⟨ξ′ ⟩ p1 (1 p1 ) 2p1 p3 ⟨ξ⟩3 p3 (1 p3 )⟨ξ⟩3 (ρE) ·[v (⃗ ρE p)] ·[keff T (τeff ·v ⃗)] S h 2pv2 ρ2E h scalar dissipation rate (εξ) was calculated according to Pirkle etal.,5 and mixture fractions in environments 1 and 2 are ⟨ξ⟩1 1and ⟨ξ⟩2 0, respectively.2.3. Population Balance Equation. In order to account forthe spatially inhomogeneous crystallization, a population balanceequation (PBE), f tN i(15)where keff is the effective conductivity and the source term (Sh)accounts for the heat of crystallization and heat of mixingbetween methanol and water in environment 3, S h S3( ΔHmix) S f ( ΔHcrys)w, j j [vjf ] [Gi(ri ,c ,T )f ] f Dt xj xj rij xj 3 B(f ,c ,T ) δ(ri ri0) h(f ,c ,T )inwhere S3 (Msolvent antisolvent in environment 3, ( jSf w,j) is the rate ofincrease in total crystal mass in environment 3, ΔHmix is the heatof mixing of methanol with water in mass basis, and ΔHcrys is theheat of crystallization of lovastatin from a methanol watermixture in mass basis. The values of ΔHmix depend on the masfraction of methanol in the mixture and are taken from Bertrandet al.32 The heat of crystallization ΔHcrys is derived from a van’tHoff relation used to fit the solubility data (see next section).2.5. Crystallization Kinetics of Lovastatin. Following thework of Pirkle et al.,5 the solubility and nucleation and growthrates for lovastatin are calculated from(12)31was used. The PBE is a continuity statement expressed in termsof the particle number density function f, which is a function ofexternal coordinates (X, Y, and Z in the Cartesian 3D case),internal coordinates ri which are the size dimensions of thecrystals, and time t.In the PBE (eq 12), the rates of growth Gi and nucleation B arefunctions of the vector of solution concentrations c and thetemperature T, δ is the Dirac delta function, and h describes thecreation and destruction of crystals due to aggregation,agglomeration, and breakage. For size-dependent growth, therates of growth Gi also varies with the ri.The PBE (eq 12), discretized along the crystal growth axisusing high-resolution finite volume method,25 was rewritten on amass basis and solved as a set of scalar transport equations in theCFD solver:rj 1/2fw, j ρc k v rj 1/2r 3f j dr ρc k vf j4c* (kg/kg of solvents) 0.001 exp(15.45763(1 1/θ)) (2.7455 10 4)W 3 (3.3716 10 2)W 2 asas 1.6704Was 33.089, for Was 45.67 (1.7884 10 2)W 1.7888, for W 45.67 asasθ 44 (r j 1/2) (rj 1/2) Tref 296 K(17) 15.8 B homogeneous at 23 C (#/(s ·m 3)) 6.97 1014 exp [ln S]2 f ] [vf i w,j Dt w,j xxx iiii 3 ρk Δr c v [(rj 1/2)4 (rj 1/2)4 ] Gj 1/2 f (f )j j2 r 4Δr Δr Gj 1/2 f j 1 (f )j 1 B ; Δc 0 j 0 2 r ρc k v Δr[(rj 1/2)4 (rj 1/2)4 ] Gj 1/2 f j 1 (f )j 1 2 r 4Δr Δr(f )j ; Δc 0 Gj 1/2 f j 2 r T,TrefB B homogeneous B heterogeneous(13) f t w, j(16) Mns ) is the rate of increase in the concentration of 0.994 B heterogeneous at 23 C (#/(s· m 3)) 2.18 108 exp [ln S]2 (18) 30G at 23 C (m/s) (8.33 1036.7)(2.46 10 ln S)(19)where Was is the weight percent of antisolvent (H2O), S c/c* isthe relative supersaturation, and c and c* are the solution andsaturated concentration, respectively, and the coefficient15.45763 in the temperature-dependence factor infers a heat ofcrystallization value of ΔHcrys 38 042.5 kJ/kmol. Thesolubility (eq 17) was fit to experimental data from three sourcescited by ref 5, whereas the nucleation and growth rate expressions(14)where f w,j is the cell-averaged crystal mass and has the units kg/m3, Δr rj 1/2 rj 1/2, ρc is the crystal density, kv is the crystalvolume shape factor, ( f r)j is the derivative approximated by theminmod limiter,32 and Δc is the supersaturation.5435DOI: 10.1021/acs.iecr.8b00359Ind. Eng. Chem. Res. 2018, 57, 5433 5441
ArticleIndustrial & Engineering Chemistry Researchsolutions were obtained by comparing the steady-state solutionfor different grid sizes.3.4. Operating Conditions Studied Here. Simulationswere performed with the solution of lovastatin/methanol(solute/solvent) fed through the main inlet at a temperature of305 K, and the antisolvent (pure water) was fed through theradial inlet at a temperature of 293 K. In all the simulations, thetotal mass flow rate (solution antisolvent) was kept constantand equal to 1.0 kg/s which corresponds to an approximateaverage residence time of 1.0 s.First, the effect of taking into account the heat of mixing andthe heat of crystallization in the temperature and consequently inthe CSD and solute conversion was analyzed. After that, differentradial inlet configurations and radial inlet velocities were studied.In this first set of simulations the methanol/water mass flow ratiowas set to 1.0. Further, the influence of methanol/water inletmass flow ratio on the mixing properties and crystallization wasanalyzed for the best radial inlet configuration.(eqs 18 and 19 are from an early experimental study onantisolvent crystallization.243. NUMERICAL SOLUTION PROCEDURE3.1. Computational Domain. The effect of the number ofradial inlets was investigated by creating five different 3Dcomputational domains, with XY z 0 plane of symmetry, withone, two, three, and four radial inlets, as well as a 360 radial inlet,as illustrated in Figure 1 for four radial inlets. In all of the4. RESULTS AND DISCUSSION4.1. Effect of Heat of Mixing and Heat of Crystallizationon Temperature and CSD. In order to analyze the effect ofheat of mixing and heat of crystallization separately, threesimulations were performed: (a) without considering the heats ofmixing and crystallization; (b) only considering the heat ofcrystallization; (c) considering both heat of mixing and heat ofcrystallization. The spatial temperature fields and mass-weightedoutlet CSDs for these simulations are shown in Figures 2 and 3,respectively.Figure 1. Illustration of the computational domains used in thesimulations.domains, the diameter and length of the main pipe were 0.0363and 1 m, respectively. The radial inlets were positioned at an axialposition of X 0.1 m. In order to keep the same specific kineticenergy and total mass flow rate at the radial inlets, the area of theinlets was calculated according to the number of inletsconsidered.3.2. Mesh. The numerical solutions were performed on 3Dcomputational meshes. GAMBIT 2.13 software was used to setup the computational grid, and the fluentMeshToFoam tool,which is available in OpenFOAM, was used to convert the grid toOpenFOAM standards. Triangular and rectangular cell faceswere used, when needed, to improve the mesh quality. Theaverage grid spacing between nodes was set to 1 mm.3.3. Model Implementation and Numerical Solution.The model equations were implemented on OpenFOAM 2.3 viaobject-oriented C programing language. A set of dictionarieswas used to input the transport, PBE, and finite-mode PDFproperties and variables. The population balance equation wasdiscretized into 30 bins for the longest growth axis, with δr 8μm. The 30 semidiscretized PBE equations, resulting from PBEgrowth axis discretization, were implemented in OpenFOAMcode using the PtrList T C template which constructs anarray of classes or templates of type T. The merged PISOSIMPLE (PIMPLE) algorithm was applied to run thesimulations. This algorithm combines the SIMPLE algorithmand then uses pressure implicit with splitting the operators(PISO) algorithm to rectify the second pressure correction andcorrect both velocities and pressure explicitly.33 The schemesimplemented for both convection divergence and diffusion(Laplacian) terms were the bounded second-order linear upwindand the unbounded second-order linear limited differencingschemes, respectively. Transient simulations were run until thesolutions achieved the steady state. Grid-independent numericalFigure 2. Temperature contour plot: (a) no heat of mixing andcrystallization; (b) only heat of crystallization; (c) heat of mixing andcrystallization.The inclusion of the heat of crystallization term in thesimulations, as expected for the system studied here, does notaffect the temperature significantly (Figure 2b), with the averageoutlet temperature being 298.6 K compared to 298.1 K when theterm is not included (Figure 2a). The heat of mixing, on the otherhand, strongly affects the spatial distribution of temperature inthe crystallizer (cf. Figure 2a,b and Figure 2c), with an averageoutlet temperature of 306.2 K. The heat of mixing also notablyaffects the CSD, as shown in Figure 3, with a narrower CSD andmuch smaller mean crystal size compared to simulations in whichthe heat of mixing is not taken into account. The soluteconversion into crystals is affected by the heat of mixing, reducingits value from 81.5% for simulation (a) to 70.4% for simulation(c). The higher temperatures caused by the heat of mixing resultin higher solubility (eq 17) and lower supersaturation and growthand nucleation rates (eqs 18 and 19), which together result inboth lower solute conversion and lower mean crysta
Multiscale Modeling and Simulation of Macromixing, Micromixing, . capacity, reduce operating costs, and identify potential opera- . the evolution of the crystal size distribution, and the energy balance equation to account for the heat transfer between the
methods [29,30], the equation-free multiscale methods [31,32], the triple-decker atomistic-mesoscopic-continuum method [23], and the internal-flow multiscale method [33,34]. A nice overview of multiscale flow simulations using particles is presented in [35]. In this paper, we apply a hybrid multiscale method that couples atomistic details ob-
1 Simulation Modeling 1 2 Generating Randomness in Simulation 17 3 Spreadsheet Simulation 63 4 Introduction to Simulation in Arena 97 5 Basic Process Modeling 163 6 Modeling Randomness in Simulation 233 7 Analyzing Simulation Output 299 8 Modeling Queuing and Inventory Systems 393 9 Entity Movement and Material-Handling Constructs 489
The multiscale model described in this paper combines continuum submodels and a discrete stochastic submodel into a multiscale modeling environment for studying P. aeruginosa swarming. At the continuum level, thin liquid film submodel is used to describe the hydrodynamics of mixt
A CRM is embedded at each grid column ( 100s km) of the host GCM to . A. and K.-M. Xu, 2011: Improved low-cloud simulation from a multiscale modeling framework with a third-order turbulence closure in its cloud-resolving model . Off-line sensitivity test to vertical resolution ATEX Cumulus ASTEX Stratocumulus . Low-level (sfc - 700 hPa .
incomplete mathematical formulations, and numerical implementations that are inconsistent with both the mathematical and physical properties of the system. In general, multiscale research ef-forts remain largely disjoint across disciplines and typically focus on only one of the two multiscale categories.
proposed algorithms. Matlab toolboxes are online accessible to reproduce the results and be implemented for general multiscale denoising approaches provided by the users. Index Terms—image denoising, multiscale analysis, cy-cle spinning, translation invariant, Gibbs phenomenon, Gaussian noise, Poisson noise, 2-dimensional image, 3-dimensional .
1) Multiscale Collisional-Radiative Kinetics ("chemistry") - Hierarchy of complexity: Single-fluid / MHD Multi-fluid / MHD / Maxwell Non-Maxwellian (discrete & MC) 2) Multiscale Transport ("dynamics") - LHS of FP/Boltzmann streaming operator - Different directions: PIC, Vlasov, Hybrid - Anisotropic effects (magnetized) 3) Multiscale Couplings (additional physics)
Quality level according to API 6A - PSL 1, 2 or 3. 1. In the trunnion mounted design configuration, the ball is supported by bearing, held in position by the valve closures. This configuration allows to discharge any side loads on the valve body, enabling a smoother operation of the ball, minimizing the operating torque and reducing seat seal wear. 2. Anti-Blow Out stem design. 3. Standard .
