## 1. Introduction

Modeling of chemical reactors attempts to solve both conservation (mass, energy, and momentum) and chemical kinetics equations [1]. The complexity of the mathematics involved can be drastically reduced by considering that convection dominates the diffusion, by assuming a unidimensional scenario or by simplifying the momentum transport equations [2]. Nevertheless, these assumptions may oversimplify the mathematical model by neglecting mixing problems. Mixing plays a fundamental role in reaction engineering. For instance, the kinetics, the molecular weight, and the composition of polymers can be altered due to local concentration gradients as a consequence of bad mixing [3].

In this study, a batch chemical reactor is analyzed. This type of reactor is defined as a closed and spatially uniform system where the chemical species are transformed only as a function of time. The transformation of chemical species can be quantified by following any physicochemical property associated with either reagents or products. During free radical polymerizations, the viscosity of the medium increases dramatically while products are formed [4]. The kinetics of polymerization can be followed from the change of viscosity.

Rapid computational development has made the numerical analysis of phenomena associated with stirred tanks easier [5]. For example, through CFD, Patel [6] studied the mixing process on a continuous stirred tank reactor and how the thermal polymerization of styrene is affected.

Computational analysis in stirred reactors has to consider at least two models: one for turbulence and the other for stirring. The turbulence model describes the random and chaotic movement of a fluid [7], while the stirring model describes the displacement of the fluid as a consequence of the local movement of mechanical parts.

The study of batch reactors with a tracer is the basis for understanding flow behavior [8]. Tracer evolution curves allow to identify regions with turbulence, dead zones, recirculation cycles, closed circuits, or even to determine the mixing time of the reactor [9, 10].

In this work, a tracer test was used to validate a mathematical model. The mixing process was analyzed by using both the experimental and simulated behavior of the tracer. The experimental kinetics of polymerization was obtained by a multiparametric nonlinear regression of viscosity-time data.

## 2. Problem definition

The uses of polyacrylamide have been extended to different applications in the oil industry, such as water conformance, fracking, and enhanced oil recovery (EOR) processes.

In EOR applications, acrylamide (AAm) polymers are dispersed in water to increase their viscosity. However, at high temperatures the viscosity of AAm polymer decreases due to hydrolysis [11]. This can be mitigated by using co-monomers such as sodium 2-acrylamido-2-methylpropane sulfonate (AMPSNa) [12]. This sodium sulfonate monomer is well known because it confers stability to the polymer against divalent cations and high temperatures (above 90°C). In view of the benefits, it is necessary to develop a process for the synthesis of the AAm-AMPSNa copolymer that guarantees product quality and synthesis reproducibility in order to properly design the polymer. The next sections will be focused on studying the relation between the mixing time and the mixing process during the synthesis of copolymer in a batch reactor.

## 3. Rheokinetics and inverse problem

The kinetics of polymerization was analyzed through the evolution of rheological behavior of the reactive system. Rheology has entered into science fields, such as biology and polymer science [13]. Historically, polymer science and rheology converge in what is known as rheokinetics. This field was created more than 30 years ago to have a better understanding of the phenomenological nature of polymerizations.

There are two main problems related to rheokinetics: the direct and inverse problems. The inverse problem, which is the principal focus of this work, deals with the determination of kinetic parameters given by the experimental data of viscosity-time curves. Equation (4) reproduces the viscosity-time curve behavior (rheokinetic model) and it was obtained from Eqs. (1–3). This equation assumes a linear free radical polymerization and it does not consider mass and energy transport effects. Equations (1) through (4) are deeply analyzed by Malkin, see [14].

where*η*is the viscosity,

*t*is the reaction time,

*x*is the conversion, [

*M*]

_{o}is the initial monomer concentration, [

*I*]

_{o}is the initial initiator concentration,

*K*,

*a*,

*b*, and

*f*are system parameters and

*k*

_{p},

*k*

_{t},

*k*

_{d}are the rate constants of propagation, termination, and initiation, respectively.

The ratio **Figure 1**. To guarantee the quality of the adjustment, *k*_{d} must be a number between 0.01 and 1 [4]. To know which process dominates the polymerization, the magnitude of

## 4. Computational fluid dynamics (CFD)

CFD is concerned with the numerical solution of the following partial differential equations that express the conservation principles of mass, energy, and momentum transport (5–7) [2].

where*x*

_{A}is the concentration of “

*A*” species,

*T*is the temperature,

**v**is the velocity,

*ρ*is the density,

*C*

_{p}is the heat capacity,

*D*

_{AB}is the diffusion coefficient,

*k*is the thermal conductivity,

**g**is the gravity,

*P*is the pressure,

**τ**is the stress tensor,

*R*

_{A}

*y*

*q*

_{R}are terms associated with the chemical reaction,

*S*

_{A}

*y*

*q*

_{I}are source terms and

Equations (5–7) are supported in two assumptions. The first one is the conservation principle which states that mass, energy, and momentum are transformed without creating or destroying themselves and the second one is the continuum hypothesis which considers continuity of its physical properties [15].

When simplifying Eq. (7) by considering a fluid of constant density and viscosity and a linear relation between the shear rate and the shear stress, the Navier-Stokes equations can be obtained Eq. (8):

Navier-Stokes equations are the basis of CFD and its numerical solution is fundamental to understand and describe the phenomena of fluid flow.

A CFD simulation is limited by the data processing rates and storing capacity. Nevertheless, improvements of computers capacity have stimulated the growth and diversification of CFD applications [16]. Nowadays, there are several CFD software tools as COMSOL^{®} and Fluent^{®}.

## 5. Fluent^{®} simulation

A typical simulation comprises the formulation of the problem, physical assumptions to simplify the mathematical model, the numerical solution of the conservation equations, data processing, and the discussion of results. The mathematical models, initial conditions, and other adjustments can be implemented through Fluent^{®} software.

In **Figure 2**, there are at least two critical stages: convergence of the numerical solution and validation of the mathematical model.

## 6. Turbulence and RANS equations

Turbulence is described as a random and chaotic movement of a fluid. Mathematically, a turbulence model is a nonlinear system in which a minimum modification on its boundary conditions produces severe alterations in the global behavior of the system [17]. Ranade proposes three approaches to model turbulence in fluids: statistical, deterministic, and structural [2]. Reynolds-average Navier-Stokes (RANS) equations are part of the statistical approximation in which turbulence is described as a combination of average variables _{f} [18]:

Semi-empirical *k*- is a turbulence model commonly used in stirred tanks. The model assumes complete turbulence and neglects molecular viscosity effects. *k*- is part of the RANS equations and it is composed by a two-equation system with two parameters to solve: *k* (turbulence kinetic energy) and * (turbulence dissipation rate). Standard **k*- was the first model; RNG (renormalization group theory) and realizable models were developed from subsequent modifications [19]. In contrast to standard *k*-, RNG improves flow with eddies and it adds a term to the equation (*R*_{}).

Realizable *k*- has a superior performance for rotational flows and for boundary layers under adverse conditions like high-pressure gradients, separation, and recirculation [18].

*G*

_{k}is the generation of “

*k*” by velocity gradients and

*G*

_{b}by buoyancy,

*S*

_{k}and

*S*

_{ϵ}are source terms,

*C*

_{2}and

*C*

_{1ϵ}are constants, and

*σ*and

_{k}*σ*are

_{k}*σ*

_{ϵ}are the turbulent Prandtl numbers for

*k*and

*ϵ*, respectively. Turbulent viscosity (

*μ*) is calculated as indicated by Eq. (17):

_{t}In contrast to RNG, realizable model uses a variable *C _{μ}* that satisfies the “realizability” through Schwarz shear rate inequality and by making the normal stress tensor positive.

*U*^{*} is calculated by Eq. (19), where **ω _{k}**:

## 7. Stirring model

CFD simulation of moving parts, e.g. impellers and turbines, requires approximations that consider the displacement and rotation of mechanical parts on a computational grid. The most used models for stirred tanks are the moving reference frame (MRF) and the sliding mesh (SM). In contrast to MRF, SM requires more computational resources and its convergence time is higher.

MRF is defined by a rotational and a stationary region. The equations are solved on a reference frame that rotates with the impeller and the problem is solved on a stationary grid [20]. When the momentum equation is solved, an additional acceleration term is incorporated in the velocity vector formulation as relative Eq. (22) or absolute Eq. (23).

The term *ρ*[**2w** × **v**_{r} + **w** × **w** × **r**] is composed by the Coriolis acceleration (**2w** × **v**_{r}) and the centripetal acceleration (**w** × **w** × **r**). The stress tensor **τ**_{r} keeps its mathematical structure, but it uses relative velocities.

SM models the rotation of the grid by adding a source term as a function of time in Eq. (8) allowing a relative movement of the adjacent grids among themselves. The SM equation is formulated for a scalar (*ϕ*) as follows:

*V*is the control volume,

*D*is the diffusion coefficient,

*S*

_{ϕ}is a source term,

**u**is the flow velocity vector,

**u**

_{g}is the velocity of the moving grid, and

*∂V*is the control volume interface.

A third manner to model the movement of mechanical parts is through the boundaries of the walls; this approximation was named tangential rotation (ROT) and holds true only for viscous flows (non-slip condition).

## 8. Numerical method

Two numerical solvers can be selected in Fluent^{®}. The first is a pressure-based solver that was initially developed for high-speed incompressible flow. The second is a density-based that was developed for high-speed compressible flow. Regardless of the solver being used, the velocity field is calculated from the momentum equations [19]. The general solution algorithm can be divided in three stages:

1. Generation of the discrete volumes (computational grid).

2. Discretization of the conservation equations.

3. Linearization of the discretized equations and solution of the resulting system.

The pressure-based solver is established from the pressure-correction equation obtained from the momentum and continuity equations. Convergence is reached when the estimated velocity field satisfies the continuity condition:

where*J*

_{f}is the mass flux and

*A*

_{f}is the surface area of face “

*f*”. In the pressure-based methods, there are segregated algorithms based on corrector-predictor approximations (e.g. SIMPLE, SIMPLEC and PISO). SIMPLE algorithm or semi-implicit method for pressure-linked equations satisfies Eq. (25) by correcting the flux

*J*

_{f}through

*J′*and by the corrected pressure

_{f}*p*’. The algorithm postulates that

*J′*follows Eq. (26) [19].

_{f}*J**is calculated by using the pressure field

_{f}*p**.

## 9. Methodology

The following materials were obtained from Sigma Aldrich (Munich, Germany) and were used without any further purification: AAm (99.8%) and AMPS (99%). Ammonium persulfate or APS (98%) was obtained from Tecsiquim (State of Mexico, Mexico). AMPSNa preparation is reported elsewhere [21]. The molar ratio of the monomers (1:1) was constant with a total initial concentration of 10.6% wt. The initial concentration of APS was kept constant at 0.5% wt. Polymerization progress was followed in an Anton Paar^{®} MCR 301 Rheometer with a concentric cylinder geometry (CC27/CX) coated with polytetrafluoroethylene. A batch of reagents was prepared and then divided into seven samples that were polymerized at different shear rates as shown in **Table 1**.

In relation to the CFD simulation, the geometry and the grid were constructed in Gambit^{®}. The dimension of the geometry described the Parr^{®} batch reactor used in the experimental work. Geometrically, the computational model is composed by a cylinder in whose interior a stirrer with rectangular impeller blades is located.

After designing the grid, sensitivity analysis was carried out to compare the velocity field magnitude in two grids with different cell sizes. The first grid (M_{200k}) contains 201,927 cells, while the second (M_{400k}) holds 482,312 cells. Pure water was used for both simulations.

Afterwards, Fluent^{®} simulations were run to select a turbulence and stirring model. Experimental validation of the computational model was done by injecting 1 mL of 1 M sodium hydroxide solution (tracer). The response of the tracer was quantified with the multiparametric device OAKTON^{®} PC 2700.

To correlate the shear rate and the stirrer speed in a batch stirred reactor, Eqs. (29) and (30) are used. This allows the comparison of stirring between the two systems used in this work, the reactor and the rheometer.

where ^{−1}] and *N* [rpm] is the stirrer speed. This relation was studied theoretically and experimentally by Sanchez, see [22]. The Re number was used to verify the turbulence.

The number 117.39 was calculated for the reactor filled with liquid water (viscosity of 0.001 Pa·s and density of 998.2 kg/m^{3}) and by using the geometrical dimension of the system.

The methodology used for the simulations is presented in **Figure 3**. The diagram contains the experimental test used to validate the simulation model.

## 10. Kinetics of polymerization

For the synthetized copolymers of **Table 1** the experimental viscosity-time curves are presented in **Figure 4**, an adjustment of experimental data was done by using Eq. (4).

The numerical results of this regression are presented in **Table 2**.

It was found that *k*_{d} values are kept constant through all experiments, concluding that the initiation kinetics is independent of the shear rate. The copolymers were characterized to obtain a better understanding of the chemistry involved in the synthesis.

All copolymers were characterized by the following techniques: Fourier transform infrared spectroscopy (FTIR) in a Cary 600 Series spectrometer from Agilent^{®}; differential scanning calorimetry (DSC) in an HP DSC 1 STAR^{®} from Mettler Toledo^{®} and rheology with a Physica MCR-301 Rheometer from Anton Paar^{®}. The results of all these characterizations are presented in **Table 3**.

Based on the experimental results, the values of the M_{v} (molecular weight), T_{g} (glass transition temperature), and T_{f} (fusion temperature) increase according to the shear rate. M_{v} of polymers obtained under C6 and C7 conditions increased 281 and 317%, respectively, compared with C2. In all experiments, the AMPSNa molar composition of the copolymer chains was relatively constant (between 42 and 50%) according to the calculated F_{2} parameter. This result was supported by the FTIR results.

Copolymers synthetized at 150 and 200 s^{−1} increased their T_{g} by 12% and 29%, respectively, considering a reference value of T_{g}=245°C. Increased values of T_{g} and T_{f} are consequences of both M_{v} [23] and stiffness of the chains. The latter is a consequence of the incorporation of sulfonate groups (e.g. AMPSNa) into the polymer [24].

## 11. Geometry, grid, and boundary conditions

The configuration of the grids for each stirring model (ROT, MRF, and SM) is shown in **Table 4**. A non-structured grid was adapted to generate tetrahedral and hybrid volumes (Tet/Hybrid).

Model | ROT | MRF | SM | |
---|---|---|---|---|

Cells | 381,352 | 201,927 | ||

Faces | 783,158 | 439,109 | ||

Nodes | 75,560 | 40,052 | ||

Volume elements | Mixed | |||

Elements in face | Triangular and quadrilaterals | |||

Volume (cm^{3}) | 3757.353 | 3758.577 | ||

Fluid regions | 1 | 2 |

The geometry generated in Gambit^{®} is presented in **Figure 5**. Zone A was defined as the stirred region, while zone B as the stationary region.

The boundary conditions were defined as “walls” for the impeller and the reactor surfaces, which implies a no-slip condition. On the top face of the reactor, a zero-shear stress boundary condition was established, implying a free fluid movement. Region A was defined as “interior” during SM and MRF simulations. Only for SM an “interface” was defined in the boundary between the stirred and stationary volumes.

Using the original geometry, the effect of resized grid cells (M_{400k} and M_{200k}) on the mathematical simulation was analyzed through the velocity variation over a defined position within the reactor; such results are presented in **Figure 6**.

## 12. Turbulence and stirring model selection

Standard, RNG and realizable *k*- were compared in terms of the convergence time. Global residual tolerances for all variables (velocity, continuity, *k*-) were kept constant at 1 × 10^{−3}. Realizable *k*- was selected due to its reduced convergence time (1305 iterations), this being compatible with the literature recommendation for high-velocity rotational flows [25, 26].

Various simulations with liquid water were performed at 100 rpm (Re=11,739) to select a stirring model (ROT, MRF, or SM). All simulations were run using a 3D no-stationary solver. Cross-sections, as shown in **Figure 11**, were used to visualize velocity profiles in the stirred tank. **Figures 7**–**9** show a typical velocity behavior of stirred tanks: a mixing zone in the upper section, high velocity gradients close to the moving blades, and a dead recirculation zone below the impeller. The mathematical model was validated by comparing the mixing time obtained against experimental data.

The graphics presented in **Figures 7**–**9** were taken from the cross-section defined in **Figures 10** and **11**. **Figures 7**–**9** show the velocity profiles of MRF, SM, and ROT models. The color scale indicates the magnitude of each vector view from two planes: Y-Z and X-Y. All images enclosed in red squares are radial cross-sections near the blades.

## 13. Tracer analysis in Fluent^{®}

Tracer simulations were developed considering an unsteady state, using as initial condition, the data obtained at the stationary state solution (*t*=0). The stirrer speed was set at 100 rpm using liquid water and tracer as materials. Initial injection positions are shown in **Figure 12**. This configuration was used in all tracer simulations. Monitors were defined as spheres with a radius of 0.8 cm and were used to track the tracer concentration in a specific region (**Figure 13**).

Each monitor registers a mass fraction of tracer every 0.55 s. Results for all simulations are presented in dimensionless concentration.

Asymptotic behavior of the tracer concentration is a criterion to define the mixing time. Tracer curves for each model are shown in **Figures 14**–**16**.

Mixing time of ROT simulations diverges drastically between injections. Tracer 1 curves start to become asymptotic at 80 s, while there is no clear tendency when the injection is made in region 2.

Results from **Figure 16** show a significant tracer concentration variation above 30 s. SM predicts that mixing time of all monitors is beyond 80 s.

SM tracer curves of **Figure 16** differ qualitatively and quantitatively from curves obtained in **Figures 14** and **15** (ROT and MRF). In general, the monitored tracer behavior is in agreement with the velocity field of a stirred tank.

Concentration contours are presented in **Table 5** for tracer 1 and 2. Mixing profiles are shown on a radial cross-section over the X-Y plane.

## 14. Validation of the computational model

The experimental mixing curves for Monitor C (Tracer 1) and E (Tracer 2) were obtained from pH data plots against time. Experimental and simulation data are compared in **Figure 17**.

Simulated curves in **Figure 17** were obtained at a stirrer speed of 100 rpm with liquid water, using MRF and realizable *k-* as working models. The injections were made in the regions established in **Figure 12**.

The mixing time for tracer 1 has a value of 15.4 s and for tracer 2, a total of 35.8 s (**Figures 18** and **19**). The criterion used to define the mixing time was the standard deviation of the concentration data registered in all monitors [26].

The experimental and numerical mixing times are shown in **Table 6**. The mixing time in Monitor C differs from the experimental value by 5 seconds, while Monitor E differs from it by 0.1 seconds. The resemblance between Monitors B-C and D-E is due to their spatial position.

Monitor | Tracer 1 (C) | Tracer 2 (E) |
---|---|---|

Experimental | 21 s | 27 s |

Simulated | 15.4 s | 26.9 s |

Global simulated | 15.4 s | 35.8 s |

The injection zone 1 (Tracer 1) allows a lower mixing time compared with the injection zone 2. However, mixing times of Monitors C and E are less sensitive to changes in stirrer speed.

The SM and ROT stirring models were discarded because the calculated mixing time does not correspond to the experimental data, as seen in **Figures 20** and **21**. Experimental mixing time is well fitted by MRF calculations.

## 15. Effect of the stirring speed

Once MRF model was selected, a stirring speed swept from 100 to 300 rpm was done. Monitor C mixing time presents minor variations between 100 to 200 rpm; however, at 300 rpm the mixing time is drastically reduced (**Figure 22**). In contrast, in Monitor E the mixing time is reduced by increasing the stirring speed. An increase in the stirring speed leads to a reduced mixing time, allowing better contact among chemical species. Nevertheless, in the case of polymerization, a high stirring speed (above 300 rpm) can produce mechanical degradation of the formed chains.

## 16. Conclusions

Chemical product design requires the development of standardized procedures to ensure reproducibility and quality of the synthetized product. If this is possible, then the impact of any experimental variation on the product properties can be properly analyzed and eventually, the optimization of the designed product can be reached.

To set up a procedure for the synthesis of an AAm-AMPSNa copolymer in a batch reactor, we have used CFD-simulations and rheokinetics. These tools were used to research the relationship between the polymerization reaction kinetics and the mixing process.

The AAm-AMPSNa copolymer properties (M_{v}, T_{g} and k_{p}/k_{t}^{1/2}) increase according to the shear rate (better mixing) in the synthesis. Specifically, the molecular weight of the polymer synthetized at the highest stirring speed (C7) increases up to 317% with respect to the lowest stirring speed in the stirred tank (C2), showing a direct relation between the mixing stirring times and the chemical kinetics.

MRF and realizable *k*-* satisfactorily model the mixing process in the stirred tank. The tracer curves obtained numerically from CFD were experimentally validated using a 1 M NaOH tracer. The simulated mixing time differs by 0.4% with regard to the experimental value of Monitor E (Tracer 2).*

According to the tracer analysis and the rheokinetics of the polymerization, it is recommended that reagents be injected (e.g. initiators or REDOX pairs) in the region defined as “Tracer 1,” operating the reactor at 217 rpm (200 s^{−1}) and controlling the temperature at 60°C.

To give continuity to this work, we suggest to include the rheokinetics model in the transport phenomena equations, to consider rheology progression and its effect on the flow pattern, as a consequence of the growing polymer chains.