| Bio-Brussels - iGEM 2022

Modeling

Throughout our project we modelled different parts to give us more insights into how our project would be needed, and how it could answer to current demand. Early models we made, were used in our human practices for example and tried to calculate the amount of salt we use in Etterbeek (where we’re from) and how much calcium carbonate we’re missing out on. To show and illustrate why water softeners are of great use in the industry we built the following model:

Model 1: Look at those pipes!

It’s one thing to say ‘calcium carbonate causes a decrease in energy efficiency and reduces lifespan of appliances”, but in this model we wanted to actually see what that meant. Usually when you heat your water (let’s say you run your dishwasher) that’s when the calcium carbonate buildup starts to come out. That’s why we’ve built a model that calculates from parameters of hardness and temperature how much of heat transfer efficiency will be lost for a family of 4. We’ve also visualized it, because models are just more fun that way.

If you just want the end result, you can play around with the model here.

However, we’ll walk you through our thinking process and how to build this model too!

A model always starts with some assumptions; the ones we made are listed here:

Assume an infinitely long pipe to make the heat loss at the extremities negligible.
Assume the outer layer of the pipe is set to a certain constant temperature. In our case it’s 100°C.
An approximation was used for the heat transfer coefficient inside the pipe because the exact calculation is too complex. This coefficient gives an idea of how efficient the heat is transferred from the inner side of the plate to the flowing water.
Let’s assume that the limescale is uniformly precipitated along the inside of the shown part of the pipe.
Assume that the hardness is caused solely by the presence of Calcium ions. Magnesium ions are also present, but in a much lower amount so we can make this assumption.

But not every pipe is the same of course; these are the parameters we used for our pipe. We modelled our pipe as a system that heats up the water inside it by having the outer layer of the pipe at a higher temperature. We will set the temperature of the outer layer of the pipe at 100 °C which means that water cannot be heated up to temperatures higher than 100 °C. The pipe is made of steel and has an outer diameter of 300 mm and an inner diameter of 240 mm. We assume the pipe to be infinitesimally long but we only visualize one meter in our model. Water consumption was set at 115m³, the average amount of water used annually by a family of four in Flanders (Vlaamse milieumaatschappij, 2022).

Fourier’s law of conductivity allows us to calculate the heat transfer across the pipe wall and the limescale layer when it is formed. What is Fourier’s law of conductivity you may ask? Here it is!

$$ Heat~flux~for~conductivity~(W) =-k*A*\frac{T_{inner} - T_{outer}}{L}= \frac{T_{outer} -T_{inner}}{R_{cond}} $$ $$ R_{cond} =\frac{L}{k*A} $$

$$ Heat~flux~for~conductivity~(W)= $$ $$ =-k*A*\frac{T_{inner} - T_{outer}}{L}= \frac{T_{outer} -T_{inner}}{R_{cond}} $$ $$ R_{cond} =\frac{L}{k*A} $$

with:

R_cond (K/W) = the resistance for heat transfer across the wall
k (W/m*K) = the material’s conductivity, an intrinsic characteristic of the material
A (m²) = the area of the pipe wall
L = the wall thickness
T_inner (K or °C) = temperature at the inner layer of the pipe, meaning the part that is in contact with the flowing water
T_outer(K or °C) = temperature of the outer layer of the pipe

Heat is also transferred via convection from the inside of the steel pipe to the flowing water, this heat flux can be calculated with Newton’s law of cooling:

$$ Heat~flux~for~convection~(W) = h*A*(T_{in} -T_{water}) $$ $$ R_{conv} = \frac{1}{h*A} $$

$$ Heat~flux~for~convection~(W) = $$ $$ = h*A*(T_{in} -T_{water}) $$ $$ R_{conv} = \frac{1}{h*A} $$

with:

R_conv (K/W) = the resistance for heat transfer through convection
h (W/m²*K) = the convective heat transfer coefficient, it’s quite complex to calculate its exact value but we can simply use an approximation for our example
A (m²) = the area of the pipe wall
T_water= the temperature at which you want to set your water to

Heat must overcome three resistances to reach the actual flow:

the conductive resistance of the steel pipe (R_cond,1)
the conductive resistance of the limescale layer (R_cond,2)
the convective resistance (R_conv)

The total resistance can be calculated as such: $$ \frac{1}{R_{tot}}= \frac{1}{R_{cond,1}}+\frac{1}{R_{cond,2}}+\frac{1}{R_{conv}}$$

Specific expressions for a cylindrical system were used to calculate the conductive resistances: $$ R_{cond,cylinder}= \frac{ln(\frac{r_1}{r_2})}{2*\pi*L*k}$$

The reduction of heat transfer caused by the limescale precipitation can be calculated by comparing the actual heat transfer (Q_real) with the heat transfer in the case that no precipitation was formed (Q₀). The only difference between the two is that in Q₀ we remove R_cond,2 out of the equation for the total resistance. The reduction in heat transfer can then be calculated as such: $$ Reduction= 1- \frac{Q_{real}}{Q_{0}}$$

Let’s talk limescale crystallization now

Calcium carbonate crystallization is a complex process which is dependent on many variables like temperature, pressure, carbonate concentration, pH and other salts. It was therefore quite the challenge to find calcium solubility values in the literature which mimicked our actual situation in Brussel. By comparing values in the literature (Coto et al.,2012) and our own EDTA titrations to calculate the hardness of our water, it was possible to deduce the solubility of calcium carbonate in function of temperature for our water. Although it is not so accurate, it does give us an idea of the upper limit of calcium carbonate solubility.

The experimental data of Coto et al. can be very nicely fitted with the following exponential function and parameters:

a=0.0007558

b=−0.01499

CaCO₃ solubility (M) =a∗exp(b∗T _water)

The above graph shows the solubility line of calcium carbonate. Above this line calcium carbonate will precipitate. However a problem was encountered when comparing these values with the hardness levels in Brussel. Our titration experiments showed our water to have more calcium ions than could be predicted by this graph. The issue lies in the fact that calcium carbonate solubility is also affected by other factors next to the temperature. Considering all possible factors would lead to a way too complex model which is why the exponential function was adjusted to better represent the measured calcium concentrations.

The model checks if the amount of calcium ions exceeds the theoretical upper limit of calcium carbonate solubility. If that is the case then the difference is converted into solid mass. By applying the density of calcium carbonate, one can estimate the thickness of the limescale layer that is precipitated at a specific temperature.

Then, knowing all that you can get to coding!

Sources

Coto, B., Martos, C., Peña, J. L., Rodríguez, R., & Pastor, G. (2012). Effects in the solubility of CaCO3: Experimental study and model description. Fluid Phase Equilibria, 324, 1-7.

Vlaamse Milieumaatschappij. 2022. Gemiddeld kraanwaterverbruik gezinnen. Available at: https://www.vmm.be/data/gemiddeld-leidingwaterverbruik-gezinnen [Accessed 11 October 2022].

Model 2: membrane reactor

Because in our proposed implementation we suggest using a membrane reactor we wanted to see how that would work for an ideal system (aka. Our protein is a super dense plate that just won’t stop removing calcium carbonate until it’s saturated, independent of the water going in). This model calculates a profile of the calcium carbonate decrease in function of the initial concentration and the flow rate.

Starting with our assumptions:

First order kinetics: Only one component is needed for the reaction to take place.
Steady State: After a short time the profile won’t change as the product will flow away and new reagent will flow into the disk. This assumption is only acceptable when there isn’t a lot of crystallization yet.
The peptide is homogenously spread in solution between the membranes, meaning that each part contains an equal amount of peptides.
No diffusion restrictions in the liquid/membrane interface: CaCO₃ will have no problems flowing into the protein suspension.
Enough carbonate ions are present in the liquid. This is needed to crystallize the calcium into calcite

Now let’s get into the modelling:

Our protein works by crystalizing calcium carbonate, with a reaction that kind of looks like this:

$$ [Ca^{2+}] + [CO_3 ^{2-}] ⇋ [CaCO_3]_{solution}$$ $$ [CaCO_3]_{solution} → [CaCO_3]_{crystal}$$

Let’s assume the second step is the rate limiting step. The kinetics could be modelled as such:

$$ r_{reac}= \frac{d[CaCO_3]_{crystal}}{dt} = k_1 * [CaCO_3]_{solution}$$

with:

k₁ (s^-1) a kinetic constant that gives an idea of how fast the crystallization process occurs thanks to the peptide

We can make a mass balance from that!

Let’s start by only focusing on an infinitesimally cylindrical disk and describing what happens inside it:

$$ Flux~of~CaCO_3~into~the~disk~(\frac{mol}{s})= -S*D*\frac{\delta c}{\delta x} + S*u*c $$ $$ Flux~of~CaCO_3~out~of~the~disk~(\frac{mol}{s})= -S*D*(\frac{\delta c}{\delta x} + \frac{\delta^2 c}{\delta x^2}) dx + S*u*(c + \frac{\delta c}{\delta x} dx) $$

$$ Flux~of~CaCO_3~into~the~disk~(\frac{mol}{s}) = $$ $$ = -S*D*\frac{\delta c}{\delta x} + S*u*c $$ $$ Flux~of~CaCO_3~out~of~the~disk~(\frac{mol}{s}) = $$ $$ = -S*D*(\frac{\delta c}{\delta x} + \frac{\delta^2 c}{\delta x^2}) dx + S*u*(c + \frac{\delta c}{\delta x} dx) $$

with:

D (m²/s) = the diffusion coefficient for calcium carbonate
S (m²) = the cross surface of the filter
u (m/s) = the speed at which water is pumped through the membrane

By constructing a mass balance, one obtains the following:

$$ \frac{\delta c}{\delta t}*V = -S*D*\frac{\delta c}{\delta x} + S*u*c +S*D*(\frac{\delta c}{\delta x} + \frac{\delta^2 c}{\delta x^2}) dx - S*u*(c + \frac{\delta c}{\delta x} dx) - r_{reac}*V $$

$$ \frac{\delta c}{\delta t}*V = -S*D*\frac{\delta c}{\delta x} + S*u*c + $$ $$ + S*D*(\frac{\delta c}{\delta x} + \frac{\delta^2 c}{\delta x^2}) dx - S*u*(c + \frac{\delta c}{\delta x} dx) - r_{reac}*V $$

A few terms will cancel each other out and the mass balance can be simplified by dividing by V= S*dx:

$$ \frac{\delta c}{\delta t}= D*\frac{\delta^2 c}{\delta x^2} - u*\frac{\delta c}{\delta x} -r_{reac}$$ $$ with~r_{reac} = k_1 * [CaCO_3]_{solution}$$

Assuming steady-state (aka no changes in time) we get:

$$ D* \frac{\delta^2 c}{\delta^2 x}= u*\frac{\delta c}{\delta x} + k_1 * [CaCO_3]_{solution} $$

That’s a differential equation! If we say that at the start of our membrane the concentration is just the original C₀ and that we have stabilized when reaching the end of the membrane. That way, we can use programming to solve it. Our solution shown below is just an example with different numbers plugged in, these values can always be changed according to the conditions under which the protein has to work.

Let’s try to solve that!

The following values were used in this example:

Diffusion coefficient D = 1.00 * 10^-09 m²/s; This is the typical value for molecules dissolved in liquids. However this value will be irrelevant in our example as the diffusion contribution is negligible compared to the convection contribution.
Initial calcium concentration C₀: 3.2 mmol/l; This is the typical value for water in Brussels.
Water velocity v: 2 m /s;
Kinetic constant k₁: As we don't know the kinetic constant, multiple values have been used. Determination of the kinetic constant requires kinetic experiment, which we didn't do yet. We did prove that our proteins work in a frame of 30 seconds.

By using this code (Athena):

Global Diff,c0,v,k1 As Real

Diff=1D-09
c0=3.2D-03
v=2

@Initial Conditions
U(1)=c0

@Model equations
F(1)=Diff*Uxx(1)-v*Ux(1)-k1*U(1)

@Boundary Conditions
If(LEFT)Then
F(1)=U(1)-c0
Else
F(1)=Ux(1)
EndIf

You get this!

This shows us that depending on the k₁ value, another thickness of the membrane filled with proteins is needed to decrease the calcium carbonate concentration to a certain threshold. Important to note is that we expect that the k₁ value is much bigger than one and thus a decrease can be observed when working with small membranes. This model doesn't tell anything about how fast the proteins saturate. We have found through experimentation that the assumption “there is plenty of calciumcarbonate present in solution” is not one that is accurate. We hope to improve that in future years by altering our designs.