COMSOL Blog » Ed Fonteshttp://www.comsol.com/blogs
Tue, 31 Mar 2015 17:56:50 +0000en-UShourly1http://wordpress.org/?v=3.9.1Why Car Batteries Perform Poorly in Cold Weatherhttp://www.comsol.com/blogs/why-car-batteries-perform-poorly-in-cold-weather/
http://www.comsol.com/blogs/why-car-batteries-perform-poorly-in-cold-weather/#commentsThu, 05 Feb 2015 09:02:57 +0000http://com.staging.comsol.com/blogs/?p=61061Starting the car on a cold winter morning can be unpleasant if you have not been proactive the night before. When you are unable to start an engine, it is often the battery’s fault. Why is a battery more sensitive than other processes in a car? The answer lies in the battery’s ability to convert chemical energy into electrical energy, with a minimum of heat generation, and the relatively small amounts of thermal energy available at low temperatures.

Getting Started

I remember one autumn a few years ago during which I bought a new car. The following winter was one of the coldest in several years. For two weeks, the thermometer in the garden showed temperatures below -10°C (14°F).

One February morning, while on a ski vacation in the Swedish mountains, I stepped out to the cottage driveway to start the car, hoping to ensure a nice and comfortable short drive for the family on the way to the ski lift. Turning the ignition, the car barely started. The vehicle made a sound revealing that the six cylinders were not running as smoothly as usual. It took almost a minute before the engine sounded as it should. Since the car was new, I found this alarming. Very slowly, the LCD between the speedometer and the tachometer came to life, showing -35°C (-31°F). No skiing this morning!

As an electrochemical engineer, my thoughts drifted from hitting the slopes to good old lead-acid battery technology, which, at the time, was able to deliver the peak current to drive the starter and start the engine in the first short turn of the key.

This problem is not limited only to batteries — the internal combustion engine also runs into difficulties at extremely low temperatures. The lubricating oil becomes thicker, the combustion reactions become sluggish, and the condensation may freeze in critical parts of the fuel system. My car, however, did start. Any electric vehicle not plugged in during such a cold night would probably not have started at all.

What is the reason for this difference? The answer is found in the way chemical energy is converted to mechanical energy:

The internal combustion engine converts chemical energy stored in the fuel to heat, which is then converted to mechanical energy.

The electric vehicle’s engine converts the chemical energy in the battery to electric energy, which is then converted to mechanical energy by the electric motor. It generates a very small amount of heat compared to the combustion engine.

The process via thermal energy to mechanical energy in the internal combustion engine gives plenty of heat from the first stroke to quickly heat up the engine, allowing the car to drive off almost immediately. Yet, the slow heat generation that occurs in extreme temperatures in the electric vehicle does not provide this same experience. To quote Les Grossman, “That is physics, it’s inevitable.”

Note that the efficiency in the conversion of chemical to mechanical energy is much higher in the electric vehicle, since the losses in the battery and in the electric motor are relatively small.

Efficiency issues and heat generation aside — and before we discuss the battery — let us compare the processes that may cause difficulties in cold weather in electric and conventional cars.

Comparing Vehicle Processes

Let’s begin with the electric motor and the internal combustion engine. We can imagine that the electric motor is less affected by low temperatures compared to an internal combustion engine. It has fewer moving parts and, since the moving parts are mainly separated by air gaps, it should require less lubrication and be less sensitive to low temperatures.

The transmission of an electric vehicle is also less complicated than the transmission of an internal combustion car, since the electric motor can operate over a wide range of loads with excellent torque. In addition, the electric car may have several motors (for example, one in the front and one in the back) and thus avoid a lot of the transmission needed for a four-wheel drive operation. This implies that the electric car does not require a complicated gearbox that has to be lubricated. Subsequently, the electric vehicle should be less sensitive to temperature for these reasons too.

Finally, the electric vehicle does not require a complicated fuel system with pumps, valves, gauges, injectors, etc. This should also make it less sensitive to low temperatures compared to a conventional car, with less components to be impeded by ice build-up.

As expected, it is the battery that performs poorly at low temperatures. In fact, the impact of low temperatures on the operation of batteries can be observed in a variety of applications, from military equipment and space applications to cell phones and house alarm keypads. This component is obviously less critical in an internal combustion engine, which only needs a short peak current to start the engine. Compare this to an electric vehicle that needs a continuous current supply. Let us therefore look more closely at the battery performance and how it is influenced by temperature.

Temperature-Dependent Properties of a Battery

A battery consists of two porous electrodes, one positive and one negative. The electronic conducting electrode material consists of packed particles of electrode material. The void between these particles produces the porosity of the electrodes (see the figure below).

The two electrodes are separated by an electrolyte. In addition, both porous electrodes incorporate pore electrolyte in the void between the solid electrode material particles. The figure below shows the discharge process in a battery, with the particle size strongly exaggerated.

The losses in the battery at a given state-of-charge are depicted in the next figure, which shows the current-potential curves for the positive (red) and negative electrodes (blue), with an operating point given by i1 and -i1 at the respective electrode. We can assume that the positive and negative electrode potentials are measured using a reference electrode in the middle of the electrolyte (see the figure above). This is to obtain the two individual electrode potentials and to include the ohmic losses on both sides of the reference electrode.

The cell potential decreases compared to the open cell voltage (see below) because of activation losses (due to electrochemical reaction kinetics), mass transport losses, and ohmic losses. Note that the cathodic current at the positive electrode is defined as negative, while the anodic current at the negative electrode is defined as positive. This is because the polarity of the electrolyte inside the battery is opposite to that of the external circuit.

Open Cell Voltage

The difference in electrode potential at zero current density is called the open cell voltage at a given state-of-charge, as shown by the previous figure.

The open cell voltage for a battery as a function of temperature at a given state-of-charge is calculated by the following expression:

where E is the cell voltage, {\Delta S} is the entropy change of the battery reaction, z is the number of electrons transferred, and F is the Faraday constant. This means that for a battery with a net discharge reaction with a positive entropy change ({\Delta S}), the cell voltage increases with temperature. For a battery with a negative entropy change, the cell voltage decreases as temperature increases.

Most lithium-ion batteries used in modern electric vehicles have a slightly negative or a very small entropy change, which means that the open cell voltage increases slightly as temperature decreases. This alone would actually improve performance at lower temperatures. However, the change in open cell voltage as a function of temperature is relatively small compared to other parameters, around 0-0.4 mV/K — that is less than 30 mV over a range of very cold temperatures (-35°C, -31°F) to room temperature. We can therefore rule out the net discharge reaction thermodynamics as a reason for poor performance at low temperatures.

Physical Properties of the Electrolyte and the Electrodes

The physical properties of the electrolyte have a large impact on battery performance. Temperature influences the conductivity and diffusivities in the electrolyte, thus also impacting the effective conductivity and diffusivities in the pore electrolyte.

The conductivity of the electrolyte can increase by one or more orders of magnitude from very cold temperatures (-35°C, -31°F) to room temperature. If we plot the logarithm of the electrolyte conductivity as a function of 1/T, we obtain a linear relation, as shown in the figure below. This figure illustrates the low conductivity at low temperatures and its exponential increase when going to higher temperatures.

Therefore, the ohmic losses (resistive losses) in the battery electrolyte increase with decreasing temperatures, which results in a lower cell voltage at a given current at lower temperatures. In addition, the poor electrolyte conductivity results in a less uniform current density distribution in the porous electrodes, which in turn lowers the capacity of the battery. The capacity is defined as the amount of ampere-hours that can be taken out of the battery before the voltage drops rapidly. At lower temperatures, the capacity is there, but the low conductivity and subsequent nonuniform current density distribution makes it unavailable until after the battery heats up.

Also, the diffusivities of the chemical species in the electrolyte, which are vital for feeding the electrochemical reactions, are lowered to the same extent as the electrolyte conductivity. Decreased diffusivities increase concentration overpotential, which decreases the cell voltage. A lowered diffusivity also decreases the capacity of the battery, since a larger fraction of the particles in the battery’s electrodes is not accessible due to mass transport limitations.

The physical explanation for the lowered mobility is that there is less thermal energy available in the electrolyte, which makes it more difficult for ions and molecules to overcome their mutual interactions or “friction”. The mobility in electrolytes as a function of temperature is described by an Arrhenius equation, where the activation energy (Ea in the figure above) represents the energy required for molecules to overcome their interaction with neighboring molecules and move around in the electrolyte.

The solid electrode material usually has a conductivity that is several orders of magnitude larger than that of the pore electrolyte. The change in conductivity with temperature in the solid material is usually negligible for battery performance. However, in some batteries, recharge may be problematic at cold temperatures, since this may lead to dendrite formation that destroys the battery.

Electrode Kinetics

The last large contribution to the poor performance of batteries at low temperatures is the sluggish kinetics of the anodic and cathodic reactions, which result in increased activation overpotential. The physical explanation for the slow electrode kinetics is that the activation energy becomes more difficult to overcome due to the lower amounts of thermal energy available in the system at low temperatures.

The figure below depicts the total effect on battery performance on account of increased activation losses, ohmic losses, and mass transport losses. We can see how the increased total overpotential at the two electrodes results in a lowered cell voltage at a given current and state-of-charge.

Thermal Management

Modern battery systems in electric vehicles are equipped with advanced thermal management systems. These systems are able to cool the battery when it operates at high loads and heat it when it’s plugged in during cold winter nights.

Heating the battery at cold temperatures also means that the efficiency and range of the electric engine declines, since a part of the electricity or regenerative power has to be converted to heat to keep the battery operating in the optimal range. Additionally, some of this power may also be used to heat up the cabin, which also reduces the efficiency and the range of the car.

The figure above shows the results from a model of a Li-ion battery pack for automotive applications, equipped with cooling and heating channels. Such models are extensively used in the design of the battery’s thermal management system.

Concluding Thoughts

The inability of electric cars to quickly and spontaneously heat up their batteries after extremely cold winter nights is due to the electric motor’s high efficiency and the fact that it does not require the production of thermal energy to transform to mechanical work. Therefore, an electric vehicle should always be plugged in during nights before ski trips such as mine, so that the battery temperature is kept at a reasonable temperature range.

If these guidelines are followed, your electric vehicle will easily start — even in the Swedish mountains. In fact, most outside parking lots in the North (such as Alaska, Canada, Sweden, and Norway) have electricity outlets, and most conventional cars are also equipped with engine heaters. You do not want to take any chances at these temperatures, not even with internal combustion engines.

If you forget to plug in your car during your ski vacation, you may be tempted to return to the comfortable cottage and perhaps think about Svante Arrhenius, the Swedish scientist who developed the first quantitative description of the temperature dependence of chemical reaction rates and transport properties.

]]>http://www.comsol.com/blogs/why-car-batteries-perform-poorly-in-cold-weather/feed/0Modeling Approaches in Heterogeneous Catalysishttp://www.comsol.com/blogs/modeling-approaches-in-heterogeneous-catalysis/
http://www.comsol.com/blogs/modeling-approaches-in-heterogeneous-catalysis/#commentsTue, 03 Feb 2015 13:48:02 +0000http://com.staging.comsol.com/blogs/?p=60571Modeling of heterogeneous catalysis traditionally attracts great interest from the chemical engineering community, due to the many industrial processes that utilize this type of catalysis. Here, we discuss the procedure of starting with detailed micro-geometries and then proceeding with approximations through homogenization. By following this procedure, from the microscopic particle level to the macroscopic reactor level, we can design the catalyst in detail and study the influence of this design on the total reactor performance.

Introduction to Heterogeneous Catalysis

Motivated by the enormous amount of production of bulk and fine chemicals and the large-scale removal of environmentally harmful chemicals in processes based on heterogeneous catalysis, the science and engineering community is very interested in this subject. There is also a large focus on understanding these processes through mathematical modeling and simulations.

Heterogeneous catalysis concerns catalytic reactions at the interface between two different phases. Here, we limit our discussion to heterogeneous catalysis at solid-fluid interfaces, which is the most common form of heterogeneous catalysis in industrial processes. In such processes, reactants and products are present in a fluid mixture in contact with a solid catalytic surface where the catalytic reactions take place.

Common large-scale industrial processes based on solid-fluid heterogeneous catalysis include the production of ammonia in the Haber-Bosch process, nitric acid in the Ostwald process, ethylene oxide in the Wacker process, hydrogen by steam reforming, and polyethylene and polypropylene in Ziegler-Natta polymerization.

Electrochemical reactions at the surface of electrocatalysts are also categorized under heterogeneous catalysis. Examples of large-scale industrial processes based on electrocatalysis are the chlor-alkali, chlorate, and water splitting processes.

The catalytic converter used in the automotive industry, where one of the steps is the reduction of nitrogen oxides to nitrogen and oxygen, is another example of a widely used process based on heterogeneous catalysis. Other important processes from an environmental point of view are desulfurization of petroleum and the catalytic oxidation of sulfur dioxide, both with the ultimate goal of reducing the emissions of sulfur oxide to the environment during combustion of petroleum products.

Sensors that are able to detect very small amounts of adsorbed molecules are often based on heterogeneous catalysis. Hydrogen and oxygen sensors are examples of such sensors working through electrocatalysis, as are biosensors based on enzymatic reactions.

The Main Steps

Catalysts in heterogeneous catalysis act in the same way as in homogeneous catalysis — by lowering the activation energy for the catalyzed reactions.

We will not discuss how catalysts actually achieve this, but we can briefly mention that most processes involve two types of solid catalysts: Acid and metal catalysts (as per G. F. Froment, K. B. Bischoff, J. De Wilde in Chemical Reactor Analysis and Design). Acid catalysts, such as aluminosilicates, may act as Lewis or Brønsted acids and in this way form a surface complex with lower activation energy. Metal catalysts are mainly concerned with hydrogenation and dehydrogenation reactions involving a surface complex of adsorbed hydrogen.

The main steps for chemical species in heterogeneous catalytic reactions can be described as follows (as per Chemical Reactor Analysis and Design and H. S. Fogler’s book Elements of Chemical Reaction Engineering):

Transport of the reactants from the bulk of a mixture to a catalyst particle

Transport of the reactants in the pores of the catalyst particles to an active site

Adsorption of the reactants to the active site

Reaction of reactants to form an adsorbed product

Desorption of the product from the active site

Transport of the products in the pores of the catalytic particle out of the particle

Transport of the products from the particle to the bulk of the mixture

The figure below shows these seven steps schematically in a cross section of a catalyst particle. The catalyst particle consists of the grains of the catalyst support. The active sites are colored in apricot and are positioned at the surface of the catalyst support (blue pattern).

Let’s look at the possible modeling approaches for these seven steps. We can also study how the seven steps are included in a model of a whole catalytic reactor, including the material balance in the bulk of a reacting mixture.

Adsorption-Desorption Models

There are two main types of adsorption: Physisorption and chemisorption.

Physisorption occurs through van der Waals forces, while chemisorption involves covalent chemical bonds. In heterogeneous catalysis, chemisorption is usually a necessary part of the catalytic reaction. Therefore, this is the only type of adsorption that we discuss here.

The most common description of chemisorption in heterogeneous catalysis is through the Langmuir isotherm. The theory for the Langmuir isotherm is based on the following assumptions:

Adsorption is energetically uniform for different adsorption sites and at different surface coverages

The Freundlich isotherm is based on the assumption that adsorption is nonuniform with some sites having higher adsorption coefficients and so being covered by molecules first. These sites have a more exothermic heat of adsorption. The heat of adsorption depends logarithmically on the surface coverage. The Temkin isotherm can be obtained by assuming a linear dependency of the heat of adsorption on the surface coverage. Other adsorption models such as the BET isotherm also include multiple layers of molecules, which have to include physisorption in their description.

Surface Reactions

According to Elements of Chemical Reaction Engineering, three possible reaction mechanisms are often used in the modeling of heterogeneous catalysis:

Single-site reaction mechanisms describe the adsorption of a molecule and its isomerization or decomposition on this single site. The figure below shows an example of molecule isomerization.

Dual-site reaction mechanisms describe the adsorption of a reactant molecule that then interacts with another free site on which the product is adsorbed. A dual reaction mechanism may also describe the adsorption of two reactants on different sites, where the adsorbed species then react to form the product, as in the figure below.

Dual- and single-site reactions are usually considered to follow Langmuir-Hinshelwood or Hougen-Watson kinetics.

A third surface reaction mechanism describes the adsorption of a reactant molecule on a single site where the adsorbed species reacts with a molecule in the fluid (not adsorbed) to form the adsorbed product. This is usually referred to as the Eley-Rideal mechanism.

Even for the simplest possible single-site reaction mechanism, the overall rate expression may be relatively complex when expressed analytically as a single expression. We can look at the example of isomerization according to the following total reaction:

The reaction rate expression even for this simple mechanism, using the law of mass action and the simplification that reactions 1 and 3 are close to steady state, is relatively complicated (as shown by G. F. Froment, K. B. Bischoff, and J. De Wilde in Chemical Reactor Analysis and Design):

In this equation, \theta denotes the surface concentration of unoccupied active sites, {c_A}, and {c_B} are the concentrations of isomers A and B respectively. {K}, {k_1}, {k_2}, and {k_3} are constants.

When using numerical methods, we do not have to derive a total rate expression for the mechanism. In this case, we can just type in the chemical equations for the mechanism. Instead, we have to give estimates of the rate constants for each of the steps and also an estimation of the surface concentrations at the start of the process.

The combination of deriving an analytical expression and entering the full reaction mechanism, without any simplifications, is also a good alternative. The reason is that the analytical expression contains products of constants and intermediate concentrations that can be obtained from measurements of the total reaction rate of the reaction. We can use these products to estimate the reaction rates for all steps in the detailed mechanism. In addition, this combined approach gives us the possibility to investigate under which conditions the simplifying assumptions are valid and when we have to use the full set of reactions in a model.

The Transport Steps

The description of the transport in the bulk and in a porous catalyst may be done at the particle level at a microscopic scale, where the catalyst particles and pore geometry are described in detail. In addition, when the catalyst is very active, then a porous structure is not required and the geometry can also be described in detail.

However, when the pores are very small compared to the size of the typical representative unit cell of the catalyst geometry, or the number of the pores in this unit cell is too large, a detailed geometrical description of the porous catalyst is not viable. At this stage, we can use a homogenized description of the porous catalyst structure.

Let’s assume that we treat the microstructure of the catalyst with a homogenized model. We can still describe the structure made up of the particles using a detailed model when we study the transport in the bulk mixture in the whole reactor.

When the number of particles becomes too large and the geometry of the particle bed becomes too complex, homogenization at the scale of the bulk mixture has to be used as well.

Let’s have a look at these four possible modeling approaches:

Detailed description of the microscopic catalyst structure

Homogenized description of every particle

Detailed description of the structure of many particles where the description of the individual particles is homogenized

Homogenization of the particle bed in combination with homogenized description of the particles

Microscopic Scale

The possibility of modeling very complex irregular structures has increased with the development of X-ray microtomography and the computing performance available for engineers and scientists working with modeling and simulations.

A small part of a porous catalyst in a battery.

In addition, the ability to manufacture microstructures with well-defined catalyst and pore geometries has also increased the interest for studying these regular structures using their detailed geometry. Let’s look at a hypothetical example of a catalyst particle constructed of a regular geometry of channels, according to the figure above. The geometry consists of around a hundred overlapping catalyst grains that together form a spherical catalyst particle of about 10 micrometers in diameter. The voids between the catalyst’s grains form a tortuous microchannel system.

The equations for the material balances describe the accumulation and flux balance in the domain. The adsorption, reaction, and desorption processes are described as boundary conditions at the surface of the catalyst grains. The grains themselves are not part of the modeling domain; only the fluid outside the particles and the fluid in the pores inside the particles are part of the modeling domain. The material balance for a species, _i, in the fluid in the absence of homogeneous reactions is:

The flux of species at the boundaries that represent the catalyst surface is then given by the reaction rate, R_i, for the consumption or production of species _i in the catalytic reaction:

(5)

\[{{\bf{N}}_i} \cdot {\bf{n}} = - {R_i}\]

where {\bf{n}} represents the outwards normal vector to the surface of the catalyst.

The figure below shows the results of a simulation of the flow around a particle and the concentration of the reactant. The surface plot shows the concentration of a species that reacts at the surface of the catalyst grains. The streamlines show the velocity field.

We can see that the concentration in the bulk around the particles varies between 0.98 mol/m3 at the inlet (at the bottom of the domain) to around 0.4 mol/m3 at the outlet (at the top of the domain). Furthermore, the highest concentration at the catalyst grain surface, which is found at the bottom of the sphere, is 0.59 mol/m3. This means that even when the fluid first makes contact with the catalyst particle at the lower part of the sphere, there is a mass transport boundary layer where the concentration varies by about 0.30 mol/m3 (= 0.98-0.59).

The concentration difference between the bulk and the surface of the sphere varies around the spherical catalyst particle. This mass transport boundary layer corresponds to Step 1 in the list outlined in the “Main Steps” section above.

Not all of the reactant molecules are consumed at the outer surface of the catalyst particle. Some of the molecules diffuse into the porous structure to later adsorb and react at the surface of the pore walls, i.e., the surface of the catalyst grains inside the porous particle. This transport process corresponds to Step 2 in the list.

Homogenization at the Microscopic Scale

If we want to model a large number of particles, then we cannot afford to model the geometry in the same level of detail as shown above. In that case, we may model the porous structure through homogenization. Essentially, this means that we treat the porous particle as a homogeneous slab containing both the fluid and the catalyst. The fluid is modeled using the porosity in the description with effective transport properties that depend on the porosity and tortuosity of the particles.

After homogenization, the adsorption, reaction, and desorption processes become sources and sinks in the material balances. In other words, they are no longer described as boundary conditions to the fluid domain. Instead, the reactions are represented as sources or sinks in the porous domain of the catalyst particles.

In this equation, {S_a} represents the specific catalytic surface area of the particles, which is the surface area per unit volume.

In addition, we now have two different modeling domains with two different sets of domain equations: One for the bulk fluid surrounding the particle, without the catalytic reaction term, and one for the porous particle. In the internal boundaries between these two domains, we have continuity in concentration and in flux of chemical species.

We also find continuity in the fluid flow equations, which for the bulk fluid are described by the Navier-Stokes equations, while the Brinkman equations may describe the fluid flow in the porous domain.

The figure below shows the flow field and the concentration at the surface and around a homogenized particle. Note that the surface now only shows the concentration at the outer surface of the porous pellet. In the detailed case above, the surface plot represented the surface of the grains all the way into the center of the particle.

The results are close to those including the full geometrical description. The maximum concentration at the surface of the pellet is 0.59 mol/m3 and the minimum is almost 0 mol/m3 in the middle of the particle, just as in the case with the detailed description. Also, the outlet of the unit cell is about 0.40 mol/m3, as in the case for the detailed description.

Why is there decent agreement? Well, it’s because the catalyst grains in the detailed geometry are relatively small compared to the total size of the particle. When the catalyst grains become larger and more irregular in shape, homogenization becomes a less accurate description.

It may also be interesting to look at the concentration inside of the particle in the homogenized case and in the detailed case. The figure below shows that the concentration distribution in the particle is also decently described in the homogenized model. It is extremely important that the homogenized model gives a good description of the transport and reaction inside the particles if we want to use it for design purposes. The reaction distribution gives us a picture of how much of the internal parts of the catalyst are used in the adsorption, reaction, and desorption processes. If the internal parts are not used, then we are actually wasting a lot of catalyst material and we may instead use a nonporous catalyst, like in the monolithic reactors. For very active catalysts, we usually do not need the enlarged surface area provided by the porous particles.

Detailed Description in the Macroscopic Scale of Many Homogenized Particles

Once we know that a homogenized particle model gives us a decent approximation of the detailed geometry, we may use this approximation for describing a large number of particles in a fixed bed reactor. An example may be a micro-reactor, where the detailed geometry may be constructed using microtechnology. Such reactors may be used for detailed kinetic studies and for lab-on-a-chip (LOC) devices.

Also, in this case, we have two types of domains: One for the bulk fluid and one for the porous particles. We also find continuity in flux and concentration across the internal boundaries between the bulk fluid and the porous particles.

Now, the concentration decreases to small values in the bulk mixture outside of the catalyst particles as well. The maximum concentration at the surface of the particles is still the same, 0.59 mol/m3. That is to be expected, as the first row of particles experiences an inlet that is almost identical to the single-particle model above.

The figure below shows two plots of the concentration from the inlet to the outlet along two vertical lines. One line runs in a plane slightly above the particles and another runs in a plane along the position of the center of the particles. We can see that the curve along the line closer to the particles (green) reveals the position of the particles with its wavy pattern. The concentration curve along the line a bit further from the particles is smooth because it is smoothed by diffusion. The average concentration along the reactor, from the inlet to the outlet, is somewhere in between the green and the blue curves.

Homogenized Description in the Macroscopic Scale

In the example above, we introduce a first homogenization for the inner parts of the catalyst particles. See the figure below:

When there are too many particles in a fixed bed reactor, the particles are also very small in comparison to the dimension of the fixed bed. In that case, we cannot afford to describe the detailed geometry of the bed. Instead, we may use a second homogenization for that and the bulk fluid in the fixed bed.

Since the porous particles are very small compared to the size of the bed, we can collapse the particles to sources and sinks in the porous bulk domain. However, in this homogenization, we have to keep the description of the transport across the boundary layer that surrounds every particle and the transport-reaction process inside the porous particles. We also have to keep the approximate description of the transport and the flow in the porous bulk domain. See the figure below for the principle depicted on the non-homogenized particles. We can achieve this by introducing an independent variable, r, in every point in the fixed bed, which corresponds to the radius of the porous particles, including the boundary layer just outside of this particle.

The transport in the boundary layer and the transport reaction process inside the particles is then approximated with a spherically symmetrical material balance along the radius. The advective flux in the detailed models in our example above mainly runs axially (in the direction down-up across the particle), while the second homogenization only accounts for radial variations in the particles. However, some of this advection may be accounted for as a slightly increased effective diffusivity, which gives dispersion in addition to molecular diffusion in the simplified description.

If we model the fixed bed in three dimensions, this actually means that we obtain four independent space variables: x, y, z, and r for the variable along the particle radius. The fifth independent variable is time. In addition, we have to use a so-called bimodal porosity description, which describes the micro-porosity inside the catalyst particles and a macro-porosity for the void between the particles where the bulk mixture flows (see the figure above).

The coupling between the 1D transport-reaction equations at the particle level and the bulk of the fixed bed is obtained by introducing the flux in and out of the porous particles as sources and sinks in the porous bulk domain. Conversely, the local concentration in the bulk domain is the boundary condition for the concentration at the boundary layer in the particle model in every point in the porous bed.

For steady-state models with relatively simple reaction kinetics described as a polynomials of one concentration, the transport reaction process along the radius of the particle (including the boundary layer) may be solved analytically. This solution is tabulated for different types of particles and kinetic reaction orders by using the Thiele modulus and the so-called effectiveness factor. In such cases, the analytical expressions for the flux in and out of the porous particles are used as sources and sinks in the material balance for the bulk of the fluid in the fixed bed.

Also, at the macroscopic scale, the accuracy of the homogenized approximation depends on the relative size of the particles compared to the size of the fixed bed. The accuracy of the approximation increases with a decreasing particle size.

We may also note something interesting. If we have small variations along the width and depth in a fixed bed reactor (or radius if it is a cylindrical reactor), we can describe the macroscopic model with a 1D model along the height of the reactor, from the inlet to the outlet. We obtain a model with two independent space variables, one in the z-direction and one for the porous catalyst particle radius, r. In the case of our microreactor above, such a simplification gives a decent agreement with the more detailed model for the average concentration along the height and for the total conversion of the reactants.

The figure below shows a case different from the one that we modeled above. Here, the reactor is in the meter scale and there are millions of catalyst particles in the fixed bed. The catalyst is less active than the cases we studied so far. The fixed bed reactor in this figure cannot be simplified with a 1D model, since the inlet nozzles (in red) influence the concentration distribution in the reactor. In this case, the model is 3D in the macroscale with an additional fourth dimension for the particle radius and a fifth independent variable for time.

Concluding Remarks

We have seen that we may describe fairly complex catalytic geometric structures at the microscopic level including the transport, adsorption, reaction, and desorption steps.

This description was approximated at a microscopic level using homogenization, which depicts the catalyst particles as porous slabs.

In the next step, we used the simplified particle model to model a large number of particles, but then describing the position and shape of every particle in detail, only the internals of the particles were simplified.

Then, when the number of particles became too large as well, we also homogenized the fixed bed domain at the macroscopic reactor level.

By following the procedure of homogenization and approximation, from the microscopic particle level to the macroscopic reactor level, we can design the catalyst in detail and study the influence of this design on the total reactor performance. In addition, this methodology can also be used ”backwards”. In other words, we can study how the macroscopic operation may influence the operating conditions at the microscopic level, thus allowing us to adapt the catalyst and bed design for the relevant operating conditions at different positions in the reactor.

]]>http://www.comsol.com/blogs/modeling-approaches-in-heterogeneous-catalysis/feed/0Geometry Modeling in Simulation Appshttp://www.comsol.com/blogs/geometry-modeling-in-simulation-apps/
http://www.comsol.com/blogs/geometry-modeling-in-simulation-apps/#commentsThu, 08 Jan 2015 09:02:57 +0000http://com.staging.comsol.com/blogs/?p=54131COMSOL applications created with the new Application Builder will make sophisticated simulations based on parameterized CAD models more accessible than ever before. A COMSOL application allows easy access to not only parameterized models but also completely different geometry configurations, such as a mixer with a variable number of impeller blades or a variable number of impellers. To make this easy for the application developer, we have made available cumulative selections and geometry subsequences. See how these tools work.

Cumulative Selections

When we create a simulation app based on a parameterized geometry model, we probably don’t want the users of our application to have to manually update all the materials, physics, and mesh settings for each new parameter value they enter. Instead, we need tools that automatically update these operations behind the scenes. Cumulative selections will help us with this — they are used to group geometric entities such as domains, boundaries, edges, and points.

For example, when a parameter value is changed, the geometry model may change in such a way that splits a boundary. What was once a single boundary may now be split in two or more parts. The application now “needs to know” what boundary conditions to apply to these newly created boundaries. If the boundary condition is tied to a cumulative selection, then this tells the application what new boundary condition to apply to the newly created boundaries.

Cumulative selections are also used to update material properties, other physics settings (e.g., source terms), mesh settings, and results. For example, different domains created in different geometry operations can be gathered in the geometry sequence according to the material of the domains. That way, all domains with the same material are available in a single cumulative selection. The assignment of a material to such a selection can be updated automatically when domains are added or removed in the geometry sequence during a parameterization.

The Geometry node in the model tree of the COMSOL Multiphysics simulation software contains an ordered sequence of geometry operations. These include the creation of a block or spheres, Boolean unions, intersections, or more high-level operations, such as fillet and loft. This sequence of operations is referred to as the geometry sequence.

Cumulative selections can also be used within the geometry sequence for pure geometry modeling tasks. For instance, a cumulative selection can gather objects to be subtracted in a Boolean difference operation. The difference operation is then automatically updated even in the case where the number of objects to be subtracted changes.

Geometry Subsequences

Geometry subsequences can be viewed as custom composite geometry objects. They can be used in a geometry sequence in the same way as primitive objects like blocks and spheres.

Geometry subsequences are defined globally for a model and can be called from the geometry sequence of any model component in order to create a geometry object. For example, a model may have different geometry objects as input to the geometry sequence depending on the type of design that the user wants to investigate. The desired object can be selected from a list of possible objects by making a call to the corresponding geometry subsequence.

A geometry subsequence can receive model parameters as input, called arguments. The output is one or several geometry objects and one or several selections that can be used as input selections by geometry operations and the material, physics, mesh, and results settings. Geometry subsequences can use cumulative selection internally within the subsequence, but they can also contribute with their output selections to cumulative selections in the main geometry sequence.

Cumulative selections and geometry subsequences are important concepts in the ability to create applications that use parameterized geometries. In fact, the geometry subsequences and the cumulative selections were both built with the Application Builder and the creation of apps in mind.

Using Cumulative Selections

The geometry sequence is based on primitive objects — such as cylinders and blocks — and composite objects — such as objects that are a result of a Boolean operation or a sweep of a cross-sectional profile. Operations that create composite objects may require selections as input.

For example, the union of two objects requires that the two objects in question are specified in a selection list. Cumulative selections make it possible to organize this process and to automatically update geometry operations that require input selections created in previous operations.

A Rushton impeller consists of a number of impeller blades attached to a disk with an impeller hub that can be attached to the impeller shaft. We can exemplify the use of cumulative selections in the geometry sequence that creates a Rushton impeller.

Geometry of a Rushton impeller.

The first step in the geometry sequence is to create the impeller disk with a parameter that controls the diameter of the disk. For this step, a work plane is used with a series of geometry operations defined with respect to this work plane. The second step is to create the first impeller blade, which is created in a work plane positioned perpendicular to the disk. The third step is to use a rotate operation in order to rotate the blade and create the desired number of blades.

The first step, with geometry operations in Work Plane 1, creates the impeller disk and an “impeller_disk” selection. The second step creates an impeller blade and a cumulative selection that we can name “impeller_blades”, as shown in the figure below.

The “impeller_blades” selection is used as input to the rotate operation. This automates the creation of the impeller blades and does not require manual input from the graphics window.

The output from the rotate operation, which creates the other impeller blades in the Rushton turbine, can also contribute to the “impeller_blades” cumulative selection. We can specify the number of blades by updating a parameter for the number of blades. The rotation angle, “angle_ib”, is then created by dividing 360° with the number of blades. The selection “impeller_blades” is then automatically updated. The figure below shows the settings window for the corresponding rotate operation.

Cumulative selections also allow for the automatic update of input selections in the materials, mesh, physics, and results settings, once the geometry sequence is finalized.

Let’s assume that we want to set the same material or the same boundary conditions on the impeller disk and the impeller blades. We can then use the output selection of the rotate operation and the disk in a pure selection operation. Doing so will make the impeller disk and the impeller blades available in the same cumulative selection. For this purpose, we can add a union selection that unites the selections “impeller_disk” and “impeller_blades”.

Note that the union selection operation is different from the geometric union operation — it makes it possible to refer to the blades and the disk with a single selection.

When the geometry sequence is finalized, the “rotating_internal_walls” selection will be automatically updated for the materials, physics, mesh, and results settings that use this selection as input.

For example, if we are modeling turbulent flow in a mixer with a rotating impeller, the input selection to the rotating interior wall boundary condition can be automatically updated by the “rotating_internal_wall” selection created in the geometry sequence. This can be done for different choices in the number of impeller blades and for different dimensions.

Using Geometry Subsequences

It is easy to imagine that we would like to have the possibility to change the type of impeller completely, depending on the type of fluid that the impeller stirs. If we are building a mixer application, we would like to provide different alternative impeller designs that the users can choose from and dimensions for their specific processes and conditions. In addition, we may also want the ability to create several impellers along the length of the impeller shaft. To make this easy, we can use geometry subsequences. They allow us to build a geometry object library that can be accessed from the geometry sequence.

In a way, the geometry subsequences can work like our own custom-created objects, which we can use in a similar fashion to the primitive objects like blocks and cylinders. However, in addition to the input dimensions and the output geometry objects, the geometry subsequences can also output selections.

The geometry subsequences are created in the Global node. The reason for having them in the Global node is that they can be called from several geometry sequences in different components. The figure below shows eight different geometry subsequences for seven different impellers and one impeller shaft.

The geometry subsequence that defines the Rushton impeller has a number of arguments that can receive input parameters that control the geometry object parameterization. Additionally, the geometry subsequence can also include a number of local parameters, which are only defined in the subsequence. For example, the number of impeller blades is received as an input parameter, denoted “n_ib”, through the arguments of the subsequence. The angle between the blades — which is obtained by dividing 360° with the number of blades — is a local parameter since it is enough to receive the number of blades to calculate the angle.

The settings in the subsequence call contain a reference to the selected subsequence — in this example, the Rushton Turbine subsequence (see the figure below). The argument expressions in the settings window now refer to global parameters in the Global Parameters node in the model tree. In this case, most of these arguments have the same name as the arguments in the geometry subsequence. However, note that the argument expressions can have other parameter names and numbers as input.

The selections created in the subsequence can be mapped to cumulative selections in the main geometry sequence. For example, the “impeller_blades” selection now contributes to the “rotating_internal_wall” cumulative selection defined in the main geometry sequence. The “impeller_hub” selection contributes to the “rotating_wall” selection in the main geometry sequence. The cumulative selections created in the main geometry sequence can then be used as input in other geometry operations, such as the union or difference Boolean operations. They can also be used in the materials, physics, mesh, and results settings.

Also, note that the subsequence call further allows for the positioning and orientation of the geometry objects created in the geometry subsequence call. Two different subsequence calls to the same subsequence can in this way get different positions and orientations in the main geometry sequence.

The figure below shows the unlikely case where seven different impeller types are positioned along a single impeller shaft. This, perhaps a bit unusual, example shows that the output selections from the geometry subsequences and the cumulative selection in the main geometry sequence can group geometric entities that are completely different in configuration. The cumulative selections can then automatically update the materials, physics, mesh, and results settings.

For example, the cumulative selection “rotating_internal_wall” automatically selects all the impeller blades and impeller disks for all different impellers in the rotating interior wall boundary condition for turbulent flow. The same functionality also groups the hub and the shaft surfaces in the “rotating_wall” cumulative selection, which is used in the rotating wall boundary condition.

Using Cumulative Selections and Geometry Subsequences in Simulation Apps

The geometry subsequences and the cumulative selections allow us to create libraries of complex and parameterized geometry objects. In addition, the ability to create cumulative selections allows us to build applications that are very user-friendly. The user of the application does not need to be bothered with irrelevant technicalities in the design of a process or a device.

As a final example, we can look at a mixer application that could be a result of the operations exemplified above.

In this app, the mixer designer can select impeller type and dimensions, vessel type and dimensions, and operating conditions for the mixer. The output of the simulations is the mixing efficiency in the process. The user of this application can focus on changing different model parameters to increase mixing efficiency using available components, with only a minimal knowledge of CFD simulation technicalities.

Additional Resources

If you enjoyed this blog post, you might also like to learn about COMSOL Server™ — your platform for distributing, managing, and running your apps.

]]>http://www.comsol.com/blogs/geometry-modeling-in-simulation-apps/feed/0Apps for Teaching Mathematical Modeling of Tubular Reactorshttp://www.comsol.com/blogs/apps-teaching-mathematical-modeling-tubular-reactors/
http://www.comsol.com/blogs/apps-teaching-mathematical-modeling-tubular-reactors/#commentsFri, 21 Nov 2014 09:02:57 +0000http://com.staging.comsol.com/blogs/?p=40057The Tubular Reactor application is a tool where students can model a nonideal tubular reactor, including radial and axial variations in temperature and composition, and investigate the impact of different operating conditions. It also exemplifies how teachers can build tailored interfaces for problems that challenge the students’ imagination. The model and exercise are originally described in Scott Fogler’s book Elements of Chemical Reaction Engineering. I wish I had access to this type of tool when I was a student!

Apps Simplify Teaching and Learning Mathematical Modeling Concepts

I still remember the calculus classes at engineering school where we first encountered partial differential equations. Despite the teacher’s efforts in trying to exemplify diffusion with the distance and the time it takes for a shark to detect your blood in the water if you cut yourself while diving, the rest of the course was mostly overshadowed by theorems. Theorems that could prove existence and uniqueness, for relatively simple problems, and by techniques such as variable separation and conformal mapping.

Apart from math theory and solving techniques, I realize now that what we really needed, in order to understand mathematical models, was to study the solution to the model equations and investigate this for different assumptions and conditions.

The Tubular Reactor with Jacket application in the COMSOL Multiphysics® software version 5.0 gives students the possibility to go from a mathematical model of a nonideal tubular reactor straight to the solution of the corresponding numerical model. The model is taken from an exercise in Scott Fogler’s book Elements of Chemical Reaction Engineering, which is one of the most popular books in undergraduate and graduate courses in chemical reaction engineering.

Value to the Student

The mathematical model consists of an energy balance and a material balance described in an axisymmetric coordinate system. As a student, you can change the activation energy of the reaction, the thermal conductivity, and the heat of reaction in the reactor (see Step 2 in the figure above). The resulting solution gives the axial and radial conversion and temperature profiles in the reactor. For some data, the results from the simulation are not obvious, which means that the interpretation of the model results also becomes a problem-solving exercise.

Value to the Teacher

The Tubular Reactor app can be accessed by a teacher in the Application Builder. As a teacher, you can investigate how to include model and application documentation in an application’s user interface. You can also learn how to include user interface commands that allow the students to generate a report from each simulation. In addition, the application accessed in the Application Builder also shows you how to create menu bars, ribbons, ribbon tabs, form collections, and forms in an application’s user interface and how to link these user interface components with settings and results in the underlying embedded model.

The Tubular Reactor Application

The different steps in the exercise for the tubular reactor problem are reflected in the ribbon on Windows® operating systems or in the main toolbar on Linux® operating systems and Mac OS in the application’s user interface.

The natural first step is to read the documentation (see Step 1 in Figure 1 above). The students can then change the activation energy and the heat of reaction, as well as the thermal conductivity in the reactor in Step 2.

The third step is to compute the solution to the model equations (Step 3). This makes it possible for the students to analyze the solution in four different plots (Step 4): Two surface plots that show the temperature and the conversion of the reactant in the reactor, and two cut line plots that show the temperature and conversion of the reactant in the reactor along three different lines placed at three different z-positions (see Figure 3 further down the page for an example). The four different plots are found under their respective tabs in a so-called form collection.

The last step is to generate a report (Step 5) that documents the model and the results from the simulation. In this case, the output is in Microsoft® Word® format, but you may also generate HTML reports.

For the teacher, the application builder tree and the member form preview in the Application Builder reveal the structure of the app (see Figure 2 below). The Main Window node (labeled 1 in the screenshot below) contains the child nodes that describe the file menu (2) and the ribbon (3). In Linux® operating systems, the ribbon is shown as a toolbar. It also contains a reference to the main form. The Form node (4) contains five forms in this case: One form that describes the main form and four forms that describe the different members in the graphics form collection. These four graphics form members correspond to the four plots mentioned above.

Figure 2. The Application Builder user interface that includes the application builder tree to the left and the preview of the included forms to the right. In between is the settings window for each selected form, declaration, method, library, or model nodes.

The text input widgets for the activation energy, the thermal conductivity, and the heat of reaction (5) are linked to the corresponding parameters in the embedded model. The range of values is also limited in order to provide a safe input range that does not produce garbage.

The Declarations node (6) includes the declaration of variables that are not defined in the embedded model. For instance, you can declare a string variable that displays a message in the user interface when the app is run based on a selection by the user. In this example, a string variable is created to show if the simulation results are updated or not (i.e., if the student changes the activation energy without re-solving the model equation, a string variable displayed in the graphics window is set to “*Not Updated”).

The application further contains a set of methods (7) that correspond to loading the model documentation, computing the results in a simulation, and generating the report. These methods are linked to the corresponding menus in the ribbon or in the main menu. The methods are graphically generated, but can then be edited manually using the method editor for further flexibility.

The Library node (8) contains files that are embedded in the application. In this example, we have a PDF-file that contains the application’s documentation linked to the corresponding ribbon menus.

The Tubular Reactor Model

The process described by the model is that for the exothermic reaction of propylene oxide with water to form propylene glycol. This reaction takes place in a tubular reactor equipped with a cooling jacket in order to avoid explosion (see the figure in the “Model Results” section below).

The reaction takes place in the liquid phase and in the presence of a solvent. The density of the reactor solution is therefore assumed to vary to a negligible extent despite variations in composition and temperature. Under these assumptions, it is possible to define a fully developed velocity profile along the radius of the reactor.

The model equations describe the conservation of material and energy. The dependent variables are the concentration, c, and the temperature, T, in the reactor. The material and energy equations are defined along two independent variables: the variable for the radial direction, r, and along the axial direction, z. These equations form a system of two coupled partial differential equations (PDE), along r and z.

The boundary conditions define the concentration and temperature at the inlet of the reactor. At the outlet, the outwards flux of material and energy is dominated by advection and is described accordingly. At the reactor wall, the heat flux is proportional to the temperature difference between the reactor and the cooling jacket.

Model Results

The results from the simulation are quite interesting. For example, the conversion profiles along the radial cut lines display a minimum and a maximum, as seen in Figure 3 below. In Fogler’s book, one of the tasks for the student is to explain these profiles.

Here, we can reveal that the profile is explained by the combination of the exothermic reaction, the advective term, and the cooling from the jacket.

In the middle of the reactor, the large flow velocity reduces the conversion, since the reactants reach far into the reactor before they react. This is labeled 1 in Figure 3.

Figure 3. Cut lines plot of the conversion in the reactor along the radial direction at different axial positions: Inlet, half axial location, and outlet.

Closer to the wall, the flow rate decreases and the conversion then increases, since the temperature is still relatively high far from the jacket wall, which also gives a high reaction rate (2).

However, as we get even closer to the wall, the conversion starts to decrease due to the cooling of the jacket, which decreases the reaction rate (3) in the figure above.

At the reactor wall, the cooling is very efficient, which should decrease the conversion even more. However, the conversion increases slightly, since there is no advection of reactants at the wall. In other words, the space time for the volume elements that travel at the wall is very high, since the flow is zero at the wall (4). The reactants are therefore consumed to a larger extent.

Applications in Teaching

The Tubular Reactor example shows how to create a dedicated user interface based on a model — an application — where students can build an intuitive connection between a physical description of a reactor and the implications of this description in the operation of the reactor. An important component in this exercise is that the results are not obvious; the interpretation of the results requires some thinking.

The Application Builder provides a user-friendly tool for the teacher to graphically create application interfaces. It allows teachers to concentrate on the exercise itself rather than investing time in explaining software tools or programming interfaces in the traditional way. They can focus on generating simulation results that trigger thinking.

The students get more challenging and entertaining exercises that focus on the problem, not on the technicalities of running simulation software.

Next Steps

Microsoft and Windows are either registered trademarks or trademarks of Microsoft Corporation in the United States and/or other countries.

Linux is a registered trademark of Linus Torvalds.

Mac OS is a trademark of Apple Inc., registered in the U.S. and other countries.

]]>http://www.comsol.com/blogs/apps-teaching-mathematical-modeling-tubular-reactors/feed/0Sometimes a Cigar Is More Than Just a Cigarhttp://www.comsol.com/blogs/sometimes-cigar-just-cigar/
http://www.comsol.com/blogs/sometimes-cigar-just-cigar/#commentsWed, 10 Sep 2014 17:22:54 +0000http://com.staging.comsol.com/blogs/?p=36897As a chemical engineer, I can’t just smoke a cigar and leave it at that. Here, I investigate the anatomy, structure, and chemical process zones of a cigar and show you a simple model of the temperature distribution of the smoke in a cigar as well as the concentration of oxygen.

A Cigar’s Anatomy, Structure, and Process Zones

Anybody who smokes a cigar every now and then recognizes the feeling of sharpness and alertness from the first few puffs, enhanced later by the calm and tranquility supplied by the following ones. Health issues aside, smoking a cigar can, in my opinion, sometimes be the best way to achieve serenity and contemplation.

In those moments of maximum alertness, I usually check that I get a straight burn line. Sometimes, this makes me think about combustion and the transport and reactions in the cigar. When calmness sets in, I usually admire the pattern of the smoke and, then, sometimes reflect about its fluid mechanics.

The schematic below shows the cigar’s anatomy with the head, body, foot, and cap; its structure with the filler, binder, and wrapper leaf; and the zones for the different processes that take place in the cigar during smoking, such as the combustion, pyrolysis, and condensation zones (references in the “Further Reading” section).

The anatomy, structure, and different process zones in a cigar. The core of the cigar is the filler, and the binder helps to keep the structure of the cigar, while the wrapper leaf envelops it. The cap is cut off before smoking. A thin straight burn line is desired.

After ignition, when the cigar has reached a pseudo-steady state, a combustion zone forms. A pyrolysis zone forms just behind the combustion zone.

Chemical Processes in a Burning Cigar

Oxygen that is transported to the combustion zone, through convection and diffusion during a puff and through diffusion between puffs, is consumed in the combustion process. Therefore, only oxygen-depleted gases reach the zone behind the combustion zone, where pyrolysis occurs. In this zone, thermal decomposition of the organic compounds in the tobacco leaves takes place in the absence of oxygen. Some of the products of pyrolysis vaporize, while others remain to be further oxidized later in what is going to become the combustion zone, as the combustion and ash zones move toward the head of the cigar. This means that distillation also occurs in the pyrolysis zone.

During a puff, vaporized products from the pyrolysis zone follow the smoke towards the head of the cigar. As the smoke travels to the head, it is cooled and filtered by the cool tobacco. Condensation of some of the compounds in the smoke is triggered by the decreased temperature, forming an aerosol. The small droplets that are formed contain the nicotine and a lot of the flavors of the mainstream smoke. Some dilution of the smoke occurs due to air leaking through the wrapper as well as some leakage of carbon monoxide out of the cigar. However, both of these processes are very slow, since the wrapper leaf is not very permeable to gas.

Also, note that the porosity is different for the filler tobacco, binder, and wrapper leaf. The porosity determines the shape of the combustion and pyrolysis zones.

What Affects the Cigar’s Shape

The shape of the combustion zone further determines the shape that the cigar takes on as you smoke it. During a puff, the combustion process close to the wrapper of the cigar is accelerated due to the high transport of oxygen around the edges. Between puffs, the edges are better cooled by the ambient air, the temperature distribution becomes more uniform, and the distribution of combustion to the middle of the cigar makes the combustion zone more uniform.

If the cigar is smoked too quickly, you’ll end up with coning, where the edges of the cigar are consumed faster than the middle. On the other hand, if the cigar is smoked too slowly, then tunneling occurs. Why? Because the middle of the cigar burns quicker between puffs due to the ambient air cooling the edges better. Imperfections in the ignition and structure of the cigar may also cause other effects of nonuniform combustion, such as canoeing and runners (read more on burn issues here.)

Personally, I think that tunneling is the worst problem, since it seems to lower the temperature, producing less and poorer quality smoke.

Cigar Simulation

Below are the results from a simulation of a relatively simple model of a cigar during a puff and between puffs. The model describes the combustion and formation of carbon, in addition to the pyrolysis of tobacco. We also included the transport of oxygen, smoke, and nitrogen in the model. Of course, in reality, the process is much more complex from a chemical engineering perspective, but the model provides a nice and simple view of the temperature distribution in the smoke and oxygen concentration.

Temperature distribution of the smoke in a cigar (top) and concentration of oxygen (bottom). The oxygen concentration dramatically decreases in the combustion zone during a puff. Leakage through the wrapper contributes only to a small extent. The combustion of carbon is a classic example of a shrinking core problem in reactions engineering.

Do Not Inhale Cigar Smoke

When the constituents of cigar smoke dissolve in the mouth, the saliva becomes alkaline. In contrast, cigarette smoke makes saliva acidic. At low pH, nicotine is protonated and is, therefore, neutral in charge. On the otherhand, at high pH, nicotine loses the proton and becomes negatively charged. The negatively charged ion of nicotine is able to pass through the membranes in the mouth much easier than the neutral form, which is repelled. For this reason, cigar smoke does not need to and should not be inhaled; nicotine is absorbed and transported into the blood just by the uptake in the mouth.

As mentioned earlier, once the first dose of nicotine has done its work by accelerating and stimulating thought, then the increased levels of nicotine replace the initial activation with serenity. This chain of events is referred to as Nesbitt’s paradox, after the scientist who first put his finger on it.

More Than Just a Cigar

In case you didn’t pick up on it, the title of this blog post is a play on the quote (or misquote, to be exact) “sometimes a cigar is just a cigar”. The original quote refers to Sigmund Freud’s answer to the insinuation that his constant cigar smoking had phallic symbolic meaning, as implied by Freud’s own theory (if you’re unfamiliar with the theory, read more here and here.)

The complex processes involved in cigar smoking and the beautiful smoke patterns show that, phallic symbolism aside, sometimes a cigar actually is more than just a cigar…

The Beckham and Maradona curl obtained with the inside of the soccer cleat (football boot), and the curl by Eder, Nelinho, and Roberto Carlos with the outside of the cleat, is due to the Magnus effect. The effect is named after the scientist who first observed it in a laboratory in the 1850s. The Magnus effect explains the side-force on a sphere that is both rotating and moving forward. Here, we use it to analyze the World Cup™ match ball.

Where Sports Meet Engineering

Like many kids around the world, my dream was to become a professional soccer player, or football as it’s called outside the U.S., Canada, and Australia. However, as a nerdy kid, I also had two other main interests: cars and science.

I remember being obsessed with the aerodynamics of cars in the early 1980s. In particular, the drag coefficient in competition and revolutionary aerodynamic designs between Audi and Ford. I was also keen to understand and achieve the specific curl of the ball when violently hitting it with the outside of my cleat. Engineering science was the glue that unified these interests later in life. Now, in preparation for the World Cup™, I will share some of my CFD analyses of the match ball with you.

Spin and Lack of Spin

The spinning of the ball has a stabilizing effect on the flow around it and thereby on its trajectory. But, let’s begin with a case where there is little or no spin.

If there is no spin, a Karman vortex street is formed behind the ball. The ball will then be subjected to force fluctuations with the detachment of vortices behind it. The wakes formed behind the ball not only increase the drag, they also contribute to the swerves that can be observed by anybody who has kicked a beach ball or been in the flight path of a knuckleball thrown by a baseball player. This semi-chaotic pattern can be partially explained in a transient simulation using the CFD Module.

The figure and animation below show the Karman vortices behind a ball rotating counter-clockwise with the same velocity at the equator as its forward movement, i.e., at the relatively low spin value. The animation is recorded on a corresponding 2D problem of a cylinder, which qualitatively shows the same effect.

Spin and the Magnus Effect

As the rotation velocity increases, the stagnation points on the ball move together and eventually outside its surface. At this point, the velocity due to the rotation of the ball is perfectly balanced by the ball’s forward movement [1]. If it were not for the fact that the ball loses some of its momentum due to friction, there would be a steady solution for this problem, which is in contrast to the case of the lower spin value discussed above. At this stage, the flight of the ball is stabilized and is much more predictable — at least for the player shooting the ball or the goalie.

The figures below show the velocity and pressure fields around the rotating forward-moving ball and a rotating cylinder. The velocity at the equator is much higher on the side of the ball that rotates and slides the air past its surface. On the other side of the ball, the ball’s rotation and the air that has to pass the ball work against each other.

Due to the difference in speed and shear between the two sides of the ball, a pressure difference is also built up between the two sides. This causes a force that sucks the ball towards the side where the air velocity is higher, which is the Magnus force acting on the ball. This is also reflected in an increased lift coefficient with the spin.

Turbulence and World Cup™ Match Ball Design

Although the simulations I shared above help explain the Magnus effect, there is more to the flight of a soccer ball than the simulations under ideal laminar conditions. Being the center of attention in the most popular sport in the world, the ball and its design have been the subjects of plenty of research. Since the release of the Adidas® Jabulani for the FIFA World Cup 2010™ in South Africa, this research has become even more intense due to the ball’s revolutionary design.

The high drag coefficients for laminar flow are caused by boundary layer separation forming a low-pressure wake that slows down the ball in this flow regime. For higher flight velocities, the boundary layer becomes turbulent ahead of the separation point and remains attached further downstream on the rear side of the ball before the flow reverses. This leads to a narrower wake and a correspondingly lower drag coefficient. The phenomenon is generally referred to as the drag crisis, which is illustrated below.

The traditional soccer ball (shown above) has 32 panels: 20 regular hexagonal and 12 regular pentagonal panels. The Jabulani, on the other hand, had eight panels, as you can see in this finite element impression of the match ball:

The reduced number of the seams, shown in black, was compensated with grooves that roughened the surface. However, the aerodynamic properties of the Jabulani were quite different from the traditional ball.

The fact that the panels and seams were fewer and shallower on the Jabulani, compared to those on a traditional ball, increased the laminar flow region with high drag coefficients, while it decreased its drag coefficients at higher speeds. The relatively large width of the laminar regime compared to that of a traditional ball gave it the behavior of a beach ball for a wider velocity range, which was the subject of complaints from many goalies. Also, the pattern that the ball presents to the wind has an impact on the sudden changes in direction for these “knuckleball” kind of shots [2].

The new Adidas Brazuca®, the official ball for the FIFA World Cup 2014™ in Brazil, has only six panels. However, the total length of the seams is comparable to that of a traditional ball. Additionally, the depth of the seams is larger than for the Jabulani.

The drag coefficient as a function of the Reynolds number for the Brazuca® is, therefore, more similar to that of a traditional ball, as you can see in the graph below. Thus, the ball is expected to have a stable flight over a wider range of velocity due to the turbulence caused by the seams.

Enjoying the Effects of the Lack of Spin, Spin, Turbulence, and the Magnus Effect

In association football, players like Ronaldo can hit a ball hard and consistently without spin, giving it a straight path far from the goal. That is due to turbulence and the swerving chaotic trajectory closer to the goal, when the flow starts to become laminar.

In contrast, the stable curl of the spinning ball, due to the Magnus effect, allowed players like Beckham and Maradona to repeatedly place a 30-meter cross within a radius of half a meter from the target.

The very hard hit together with the spin that Eder, Nelinho, and Roberto Carlos use, in combination with the transition between laminar and turbulent flow, can give the ball an ordered but peculiar trajectory, almost like a guided missile.

Right after impact with the outside of the cleat, when the ball has been accelerated to its top speed, turbulence around the ball and its small drag coefficient gives it a relatively straight trajectory. When the ball slows down, the relative spin becomes stronger and the Magnus effect becomes more pronounced. In other words, the ball first goes straight and then curves suddenly closer to the goal.

This combination of turbulence and the Magnus effect is observed in Roberto Carlos’ famous free kick in the game between Brazil and France in 1997. France’s goalkeeper, Barthez, did not even move until it was too late and a ball-boy standing several meters from the goal ducked his head. Both keeper and ball-boy thought the ball was going far outside the post!

More incredible goals featuring the Magnus effect can be found in this clip, where players such as Messi, Ronaldo, Ibrahimovic, Ronaldinho, Beckham, Eder, Cruyff, and many others make use of this effect to cheat the goalkeepers. One minute into the clip, you can also see Nelinho’s goal against Italy that, in 1978, got me hooked on hitting the ball with the outside of my cleat. It was so fascinating that I would practice that one shot several hours per week…

Cars, Science, and the World Cup™ Match Ball

In the early 1980s, car ads always featured the car’s drag coefficient. I have always wondered why this relevant measure disappeared from the published specifications for cars. Now, I can instead think of the curves for the drag and lift coefficients for the ball and picture their combination with the Magnus effect. Imagine that while you watch the incredible shots and goals in this year’s World Cup™.