Modeling Bubble Dynamics in a Stirred Tank Reactor Using ANSYS Fluent and Population Balance Model

Hello and welcome to this short tutorial on how to set up a population balance model in ANSYS Fluent. So, for this problem, I'm using a stirred tank reactor.

We are not really implementing any reactions, but we have a two-phase model, basically a multi-phase model using the Eulerian model in ANSYS Fluent. So, to give a brief overview, this is our geometry, which was adapted from ANSYS Fluent's PyAnsys library, actually.

And I added a few inlets for water flow. So, water is the primary phase in this model. It's allowed to move from the top, as you can see here, from this inlet. If I can properly select it, right here.

So, this is the inlet for water, and then the outlet for water is mainly at the bottom of the domain. And we also have a stream of air, which is our population balance phase model. The domain is an outlet boundary condition as well.

So, this is what we have in this sketch, in this schematic, showing where air is added into the domain. But here, we are not adding air as a continuous, I mean, it is a continuous phase, but it's not being added in the same way as water.

The population balance model is a popular option for modeling things like droplet dynamics, as well as particle dynamics, either solid or liquid or gaseous particles.

So, the population balance model that we are implementing has basically two phenomena that we are incorporating, the coalescence of bubbles, smaller bubbles forming bigger bubbles, as you can see here, and we can also have a breakage, which is when you have a high enough shear stress that a bubble, instead of just deforming, will actually break into two bubbles.

This is implemented within Fluent using the PBM population balance model. So, basically, the population balance model is a Lagrangian-based model, where you're tracking a population of particles at each cell. So, you're tracking concentration-based number of particles per unit volume.

And you basically have models to account for how they grow, how they form, how they shrink, or how they can even disappear. The frequency of aggregation is used as a synonym for coalescence or agglomeration.

Breakage or rupture or division of particles can be determined by the physical conditions in the fluid around it, the water phase, as well as the intrinsic properties of the bubbles.

The population balance model follows a discrete method, which is basically a model whereby we are tracking the bins of particles at different sizes. I'm going to show how that looks like in Fluent. The first thing is, once we load the model, and I've already done this part with the mesh.

This is a relatively simple design with one impeller. You might see situations where there might be several impellers, much more complex geometry than what we have here. So, basically, we have an impeller in the center near the bottom of the domain.

Let me show the tank as well, so you can see the impeller with the shaft.

And around the impeller, you have the MRF body, which we will use to implement a moving reference frame, so that we can implement motion, basically, within our model, and shear stresses that will affect the flow profile within the tank.

Once the model has been loaded in the mesh, and this is a relatively simple geometry to mesh, we can add the materials. First, I added Air, which was already the default.

I added WaterLiquid from the Fluent database, which you can go by clicking Air, and click on Fluent database, and then select WaterLiquid, and then you can just copy. So, I didn't change any properties for WaterLiquid or for Air; it's a relatively simple simulation, just to show an example.

And then I come into the model, the multi-phase model here, and then I select Eulerian.

I've already done that, and I left the default as implicit, and the number of Eulerian phases to 2. So, we have two Eulerian phases, and the first phase is water liquid; we can rename that to water, or we can keep it as it is.

So, in phase two, it is air, and here it defaults to a diameter of one meter, basically, when you activate the population balance model. Let's go to the Phase Interaction tab, where you can see the interaction between the two phases. There are many things that you can add to your phase interaction.

There are forces, heat, mass, and reactions as well. There are a lot of things that we can explore here. But for now, I kept it very simple; I kept the default Schiller drag coefficient, and I added a surface tension coefficient for air and water of 0.072 N per meter.

I could have added a lift coefficient as well, and a virtual mass coefficient. Ideally, for the initial stages of modeling, of running the simulation, it would be ideal to keep it simple, and then once we get closer to convergence, we can start activating more complexity to our model.

Finally, the population balance model. This is a tab that's available. It's set to off by default. To activate, you just select one of the methods. Here, I selected the discrete method, which I mentioned before; it's basically you're tracking bins of particles, and each bin has a specific size.

And to know basically what size you want, or what size is physically reasonable, you can select the size calculator here. So, I selected the value of the minimum and the maximum values for the bubble size in this case. It accounts for the properties of the primary and the secondary phase.

Then, I click calculate, and it will tell you about the minimum and the empirical limit for the calculated diameter. So, I selected that value here. I also used a value in the maximum diameter, somewhat closer to the maximum stable bubble size based on the capillary length.

However, it's a bit higher than that. To change the maximum diameter, you can change the number of bins, or you can change, alternatively, the ratio exponent. If I change the ratio exponent here, for example, if I go to 1.4, it's going to automatically change my maximum diameter as well.

So, I can go back to 1.2, and then we can activate the phenomena. The phenomena is basically all the kinds of changes that can happen to your bins. Aggregation Kernel is the coalescence between bubbles.

Here, I selected the Luo Aggregation Model, which is a popular model that accounts for a variety of physical phenomena to track the aggregation of particles, aggregation of bubbles, so the shear stresses and other effects.

The breakage model from Luo is also a reasonable approach; it's been used in literature quite extensively. Those are the two models I've selected, and there are a variety of other models as well that you can select.

In this case, I don't have nucleation; nucleation is basically the formation of new bubbles, or if you're looking at a different application, such as aerosol formation, nucleation would play a role.

So, here, to activate nucleation, you would have to check this box, and either set it as a constant, or create a user-defined function for nucleation specifically, and also growth rate.

So, growth rate is when basically the particles are growing, as they, for example, in the case of aerosol, or basically, they're gaining mass through condensation or other phenomena. So, in this case, I didn't activate it, because it's not applicable for this problem.

So, once all this is set, we can click Apply, and then Close. And then I can close this as well. I should actually go back there; I just want to mention something that's quite important when setting the bins. You can verify the bin size using the print bins.

More bins may be more descriptive, but more expensive computationally. So, I can close that. And I set energy as off. I also activated the SST-K-Omega turbulence model. Basically, kept all the defaults, except for, I activated the dispersed turbulence multiphase model.

Basically, this model is good for when you have one phase that's relatively small or dispersed, in this case, bubbles. So, it is a suitable model, and I've observed it helping in terms of convergence as well. So, we can click OK here, and that is basically it, in terms of the models and materials.

We've defined all the materials already, with their expected properties. So, now let's go into the cell zone conditions. So, the cell zone that we want to specify is the MRF cell zone, the zone in red that you can see here.

If you double-click on the impeller, you can check the frame motion, axis of origin, axis of rotation, and axis of rotation for the impeller itself. And so, that's set at 0 and negative 1, and I set a rotational velocity in radians per second.

The units you can actually change, if you go on domain units, if you select angular velocity, you can actually change to RPM or degrees per second. The default is radians per second. For this simulation, I set it to 30 radians per second.

And the rotation axis direction will help you define whether you're going counterclockwise or clockwise as well. The reason why we choose Z is because we are on the Z axis, that's where the impeller is aligned with the shaft. I also have to make sure that the shaft moves accordingly.

I went into wall shaft here. In wall shaft, you can basically do something very similar. And the rotation speed is 30 radians per second. So, that's all good for now. So, once we've defined all these parameters, what we need to do now is make sure that we define the boundary conditions.

However, Phase 2, I did specify a velocity of 0.5 m/s. We can play around with this velocity to see the effects of bubble size and bubble distribution as well.

And if you go into multi-phase, you can actually define what kinds of bins, or what kind of bubbles, you're going to add, in terms of the fraction of each.

So, the volume fraction for the bubbles, for the PB phase, is 1, and I set everything to 10%, so basically, we have a uniform injection of bubbles into the tank.

Also, the water injection, or the water inlet, is set at the top of the domain, as I mentioned, and we can go into phase 1, that's the water phase.

I set a mass flow rate of 6.57, because it matches closely with the expected flow rate at the bottom here of water, since we have water flowing at the bottom under a pressure outlet boundary condition.

So, basically, that helps keep the system, the continuity, the mass balance for water in the system. The next step is to define some report definitions.

It's good for checking convergence, and if you want to make some adjustments to the model, as well as it goes, and increase its complexity, and add more physical phenomena.

So, basically, what I'm doing is, I'm measuring the overall diameter of bubbles, as they emerge through this pressure outlet here at the top. So, I want to know how big my bubbles are, my air bubbles are, at the wall liquid level, basically, at the very top.

To do that, you can choose a field variable, properties, and then diameter, once you select phase 2. This is an area-weighted average. I can click OK. I can also create another report definition that's quite useful. For these problems, it's good to check on the shaft torque.

So, I will select moment, I will keep the moment axis as it is, I will write torque, and I'll select the impeller wall as my area, where I am evaluating the torque, and then I will click OK. That is pretty much it, in terms of models.

I should also mention that we are running a steady-state simulation. Once we are done, we can actually do some post-processing, looking at things like the contours for velocity, and for diameter. I created already some surfaces to do that. I'll initialize the problem, and then I'll run.

Basically, we can see the torque, as well as the diameter changing a bit. I initialize the domain with a small amount of bubbles, that's why we have some change here in the diameter, right in the beginning, to help with convergence, especially. Here are the residuals.

It can take a little bit of time. The population balance residuals can take a while, and they fluctuate quite a bit. I've already basically ran all these models, so I can actually go forward, and show what we are going to get, once we run them. I can compare the results for other angular velocities.

Here, we see the contours for diameter, at a plane crossing the tank. I can look at things like the path lines, as well, to see where my bubbles are going. So, I'm going to track in phase 2. You can see, I should actually draw a mesh, just so that we can figure out where we are at.

So, basically, you can see that there's a lot of movement here, as the gas bubbles get released from their injection point, they get affected quite significantly by the impeller.

One thing to note, if you go with a smaller angular velocity, in this case, I ran with 10 radians per second, so a third of the value, once the bubbles are released from their inlet point, they just go up immediately. They don't really get affected as much by the mixing action of the impeller.

I also want to briefly talk about comparing the diameter distribution of the bubbles, for the different cases. So, basically, what we have here is the profile, the contours. We can compare 30 radians per second, with 20 radians per second.

Again, it's a bit darker green, indicating that it seems to have smaller bubbles. And if you go with 10 radians per second, you actually get a lot more concentrated plume of bubbles, basically moving upwards. So, usually, the bigger the bubble, the faster they go up.

We can also look at things like the velocity profile. Here, I set a limit of 0 to 3 m/s, so that we can keep a consistent scale. Here is the velocity, as well. And here is the velocity, which I should name accordingly.

So, if you go from 30 to 20 to 10, you can clearly see that, as you increase the angular speed, you get more mixing, more shear. And that has an impact on the collision rate of bubbles, as well as on their probability of aggregation or breakage. And that's what we have here.

This is basically the 10 radians per second case. A useful way to evaluate the bubble distribution is using the results tab, in the model-specific section. You have the population balance option here, and we're not using moments.

If we were using moments, we could check this, but we can go directly to number density, since we're using a discrete method for population balance. So, here, what we can do is, we can look at the discrete number density, at a surface, or for a volume, for example.

For instance, if I look at the midplane, which is a plane, basically a plane that's cutting the domain. And I can actually show the midplane is here. So, we could actually cut the domains in different sections, like at different heights, and we could actually look at the distribution, for example.

We can either print, or we can plot, or we can write it as well, and do the plot in Excel. I can do the plot right here, and you can see, actually, the distribution of particles, the overall distribution of particles, which is somewhat representative of what we would have across the domain.

So, basically, we have a lot of smaller diameters, and I should also change the x-axis into exponential, so that we can read what the diameters are. Let's see how many particles are found.

Smaller diameters are 1 x 10 to the power of negative 3. That's the most abundant number of particles that we have. Please contact us at https://ozeninc.com/contact for more information.