COMSOL Under Pressure

Sometimes you need to model the concentration of a gas in a liquid or a gas-permeable solid. You can work out simple models by hand using mass transport principles, but more complicated models require finite element software. Unfortunately COMSOL is ignorant of Henry’s law, but you can bring it up to speed through the clever use of boundary conditions.

A Little Background

The idea that a gas exerts pressure is familiar to anyone who has ever over-inflated a balloon. The gas molecules bounce around, hitting each other and the wall of the balloon. A gas can dissolve in liquids and some solids (e.g., polymers) and because the gas molecules are bouncing around they exert pressure on the liquid or solid molecules as well. Even though the gas is dissolved, it still has a pressure.

Mixtures of gasses work the same way. Air is made up of a mixture of gases and each gas contributes its own partial pressure.

$$P_\mathrm{air}=P_\mathrm{O_2}+P_\mathrm{CO_2}+P_\mathrm{N_2}+P_\mathrm{H_2O}+P_\mathrm{other}=760~\mathrm{mmHg}$$

where $P_\mathrm{O_2}$ is the partial pressure of oxygen, $P_\mathrm{CO_2}$ is the partial pressure of carbon dioxide, and so on. The partial pressure of oxygen in air is about 160 mmHg.

Imagine that a beaker half full of water has been sitting in a vacuum chamber overnight so that no gas is dissolved in the water (Figure 1a). You take the beaker out and let it sit in air at normal temperature (20 °C) and pressure (760 mmHg) until it equilibrates, as shown in Figure 1b. Once the system reaches equilibrium, the partial pressure of oxygen will be the same in the air and the water, but the concentration will be different. The water molecules take up space, so there cannot be as many oxygen molecules per unit volume in the water even though the partial pressures are equal. At normal temperature and pressure (20 °C and 760 mmHg) the concentration of oxygen in air is 8.61 mol/m$^3$ and in water at is 0.21 mol/m$^3$.

Henry’s law is used to relate the concentration of a gas to its partial pressure. The affinity that a gas molecule has for a liquid or solid phase is called its solubility coefficient, or $\alpha$. Henry’s law is

$$P_A = \frac{C_A}{\alpha}$$

where $P_A$ is the partial pressure of gas $A$, $C_A$ is the concentration of gas $A$, and $\alpha$ is solubility coefficient of dissolved gas $A$. Some common solubility coefficients are listed below.¹

Why Use Henry’s Law in COMSOL?

Unfortunately COMSOL is unaware of the concept of solubility. If you model the beaker in Figure 1 in COMSOL using steady-state analysis you will get the results shown in Figure 2.

These results look good except for the fact that COMSOL is lying. COMSOL is telling us that there are 8.61 mol/m$^3$ of oxygen in the air above the beaker and in the water in the beaker. If you read the section above on Henry’s law, then you know that this is codswallop. It is true that the partial pressures of oxygen are the same in the air and the water, but not the concentrations. Instead, we have to use Henry’s law to figure out what the concentration is in the water. In this case the pressure of oxygen in the air, $P_\mathrm{O_2}^\mathrm{A}$, is equal to the pressure of oxygen in the water, $P_\mathrm{O_2}^\mathrm{W}$, or

$$P_\mathrm{O_2}^\mathrm{A} = P_\mathrm{O_2}^\mathrm{W}$$

Using Henry’s law we have

$$\frac{C^\mathrm{A}_\mathrm{O_2}}{\alpha_A}=\frac{C^\mathrm{W}_\mathrm{O_2}}{\alpha_\mathrm{W}}$$

where $C^\mathrm{A}_\mathrm{O_2}$ is the concentration of oxygen in the air, $C^\mathrm{W}_\mathrm{O_2}$ is the concentration of oxygen in water and $\alpha_W$ is the solubility of oxygen in water. The $\alpha_A$ constant does not really exist because solubilities are measured with respect to the gas phase, so we can assign it a value of one. We can find the concentration of oxygen in water with

$$C^\mathrm{W}_\mathrm{O_2}=C^\mathrm{A}_\mathrm{O_2}\alpha_\mathrm{W}$$

In other words, all we have to do to find the concentration of oxygen in water is scale the solution that COMSOL finds by the consant $\alpha_\mathrm{W}$. If we multiply 8.61 mol/m$^3$ by 0.024 then we get an oxygen concentration in water of 0.21 mol/m$^3$. You could take this approach in all your COMSOL simulations, however a much better way is to make COMSOL do it for you.

How to Use Henry’s Law in COMSOL

You can get COMSOL to implement Henry’s Law via boundary conditions. The trick is to (1) assign a physics to each material in your model and (2) assign a flux boundary condition at each interface between materials. In the beaker model above, there are two materials (water and air) and one interface between them. Since we have two domains in the COMSOl model we need to assign the flux boundary condition in each domain. The formula for inward flux at the boundary follows this form:

$$N=M\left(C_{adjacent}-C\frac{\alpha_{adjacent}}{\alpha}\right)$$

Here $C_{adjacent}$ is the gas concentration in the neighboring domain, $C$ is the gas concentration in the reference (i.e., current) domain, $\alpha$ is the gas solubility in the current domain, and $\alpha_{adjacent}$ is the gas solubility in the adjacent domain. You can generate the two boundary conditions for the beaker model using this pattern. In the air domain the boundary condition at the air/water interface is

$$N=M(C^\mathrm{W}_\mathrm{O_2}-C^\mathrm{A}_\mathrm{O_2}\alpha_\mathrm{W})$$

and in the water domain the boundary condition at the air/water interface is

$$N=M\left(C^\mathrm{A}_\mathrm{O_2}-\frac{C^\mathrm{W}_\mathrm{O_2}}{\alpha_\mathrm{W}}\right)$$

The variable $M$ is the stiff spring velocity and should be a fairly large value. I used 100000 m/s. (Update: a large value may not always work—see note below²)

Now when you run COMSOL you get the result shown in Figure 3. Each domain has a uniform concentration of O$_2$, but the concentrations in each domain are different, as expected (8.61 mol/m$^3$ in the air domain and 0.21 mol/m$^3$ in the water domain).

This approach to modeling solubilities in COMSOL can be extended to models with many different domains. Just be sure that you follow the correct pattern when assigning your flux boundary conditions.

A Quick Note on Plotting Your Results

You will need to assign a physics to each domain (material) in your model if you want to account for solubility in COMSOL. The water beaker model shown above consisted of two domains, each assigned the transport of dilute species (tds) physics. By default, COMSOL will produce two output plots for this model: one for the air domain and one for the water domain, shown below.

There are a variety of colors in these surface plots, but they are due to rounding errors and in fact represent the same concentration of oxygen in each domain. Incorporate both of these results into a single plot by following these steps:

1. Under Component -> Definition add two variables. In one of the variables enter a Name that will represent your composite plot. I used $c_{all}$. Under the Expression column enter one of your concentration variables. I used $c_w$ for water concentration in Variable 1. I did the same for Variable 2, but used $c_a$ as my Expression value.
2. Go to Results and add a new 2D Plot Group. Under this new group add a new Surface plot.
3. In the new surface plot settings go to the Expression heading and enter the variable name for your composite plot. In my case I entered $c_{all}$. Run the plot and you will have an output like that shown in Figure 3.

Supplementary Information

• Link to COMSOL file for simulation without solubility.
• Link to COMSOL file for simulation with solubility.