ABSTRACT

The literature model studied in this article describes bubble formation and growth in a highly viscous polymer liquid with support of a gaseous matter dispersed under pressure before foaming. The foam growth is induced by the application of vacuum and mass transport of volatile components dissolved in the polymer liquid. Based on this literature model, aeration processes are calculated for intermediate viscosity and low viscosity biological systems, as they are of interest for biomatter foams, in particular for food foams in industrial processes. At the end of this article, the numerical results are presented and discussed.

Keywords:

Foaming Process, Conservation Equations, Finite Difference Method, Convergence Analysis, Reynolds Number, Newtonian Liquids

1. Introduction

Foams are of essential importance in nature and in a wide range of technical fields that deal with matter of biological origin such as biotechnology as well as food production and beverage industry. There are several methods to incorporate bubbles within the biomatter structures: mixing, steam generation, gas injection, vacuum expansion and fermentation to name the most important procedures [1] . Recent experimental studies [2] - [7] have established relationships between process parameters and dispersion characteristics. However, there is still a lack of mathematically well-founded physical models that describe bubble formation and behavior in biomatter adequately.

Piesche, Nonnenmacher and Schütz [8] [9] [10] developed models to calculate foaming processes of polymer liquids in a static degassing system, which is supported by a gas dispersed into the polymer liquid under pressure before foaming. The foam growth is induced by a pressure drop and diffusion of volatile components dissolved in the polymer liquid. The related fluid mechanical transport processes are divided into two consecutive models that concern spherical foam and polyhedral foam, respectively. The foam growth is calculated by solving problem-specific conservation equations for mass, momentum and thermal energy. In the present article, the models proposed in [8] [9] [10] are adopted to calculate aeration processes of biomatter. Usually, aeration occurs by applying vacuum or including the gas phase into the liquid phase with a high pressure. Since material parameters for modeling diffusion processes inside biomatter foams do not exist in the literature, they are estimated to illustrate the influence of possible variations.

The model equations [9] are first specified, and the numerical procedures used for solving the model equations in Matlab [11] [12] [13] [14] [15] are outlined. A convergence analysis [16] [17] confirms that the numerical method converges with order 2 in space and at least with order 1 in time. In contrast to the literature [8] [9] [10] , where foaming processes of highly viscous polymer melts are calculated and the Reynolds numbers are Re $\ll$ 1 so that inertia terms in the momentum equation can be neglected [8] (p. 63), foaming processes of food systems [1] - [7] with Reynolds numbers Re ≈ 1 are now calculated. The Reynolds numbers Re ≥ 1 require that inertia terms in the momentum equation have to be taken into account since inertial forces and viscous forces are likewise relevant. At the end of this article, the numerical results are plotted, discussed and compared both with the results presented by Piesche, Nonnenmacher and Schütz [9] [10] for polymer foams and the results given by Haedelt, Beckett and Niranjan [2] [3] for aerated chocolate as a representative of matter of biological origin.

2. Model Equations for Foaming Processes

For the sake of completeness, it is stated that the following considerations postulate physically existent, stable and mathematically unique solutions of the non-linear transport equations for mass, momentum and thermal energy [8] [9] [10] . As mentioned above, mimicking foaming processes occurs by employing two models that describe spherical foam and polyhedral foam, respectively. The model for the spherical foam is based on the assumption that, at the beginning of a foaming process, small monodisperse spherical bubbles are already present in the highly viscous liquid. The bubbles are arranged in a regular cubic spatial structure, see Figure 1 (left). Each elementary cell includes four bubbles, which are composed proportionately. The foam volume of an elementary cell at time t = 0 is given by

${\text{V}}_{\text{ges},0}^{\text{e}}=\frac{16}{3}\text{π}\frac{r_{B,0}^3}{{\text{c}}_{\text{V},0}}$ , (1)

where c_V,0 is the gas volume fraction and r_B,0 is the bubble radius before foaming. The gas volume fraction at time t = 0 is given by

${\text{c}}_{\text{V},0}=\frac{{\text{V}}_{\text{g},0}}{{\text{V}}_{\text{ges},0}}=\frac{{\text{V}}_{\text{g},0}}{{\text{V}}_{\text{g},0}+{\text{V}}_{\text{f}}}$ . (2)

The total foam volume V_ges,0 at time t = 0 is composed of the gas volume V_g,0 of the bubbles and a liquid volume V_f. The liquid volume that surrounds the bubbles of an elementary cell is constant, regardless of the time:

${\text{V}}_{\text{f}}^{\text{e}}=\frac{16}{3}\text{π}{r}_{B,0}^3\left(\frac{1}{{\text{c}}_{\text{V},0}}-1\right)$ . (3)

The first model is used to describe the initial phase of foaming until the bubbles begin to influence each other and polyhedral bubbles are formed. In the second model, the monodisperse spherical bubbles are substituted by monodisperse pentagon dodecahedron bubbles. The surface O_D and volume V_D of a dodecahedron bubble are calculated from the edge length a_D of the equilateral pentagon, i.e.

${\text{O}}_{\text{D}}=3\sqrt{25+10\sqrt{5}}\cdot a_D^2$ , (4)

${\text{V}}_{\text{D}}=\frac{1}{4}\left(15+7\sqrt{5}\right)\cdot a_D^3$ . (5)

The constant liquid volume ${\text{V}}_{\text{f}}^{\text{e}}/4$ is arranged in a thin layer normal to the surface of the respective bubble. The next three subsections are dedicated to the model equations, which are adopted from the literature [8] [9] [10] .

2.1. Model Equations for the Spherical Foam

The conservation equations for the liquid phase and the transport equations for mass and thermal energy are formulated in spherical coordinates (r, θ, φ). A pure elongational flow is assumed to arise. The velocity of the liquid in radial direction is only relevant and is denoted by u. The velocities in polar and azimuthal direction are zero. Figure 1 (right) illustrates the growth of a single spherical bubble during foaming. At time t = 0, a bubble consists completely of one gaseous matter, for times t > 0, the volatile components dissolved in the liquid diffuse into the bubble. In every time step, the outer radius r_S(t) of a constant liquid volume that surrounds the bubble is determined as a function of the current bubble radius r_B(t):

$r_S={\left(r_{B,0}^3\left(\frac{\text{π}}{3\sqrt{2}}\frac{1}{{\text{c}}_{\text{V},0}}-1\right)+r_B^3\right)}^{1/3}$ . (6)

Equation (6) determines the integration limits for the conservation equations. As shown in [9] , the conservation equation for mass in the liquid phase simplifies to an expression for the radial velocity component (see also Figure 1, right)

$u\left(r,t\right)=\frac{\text{d}r}{\text{d}t}=\frac{r_B^2}{r^2}\cdot \frac{\text{d}{r}_B}{\text{d}t}$ . (7)

Similarly, the simplified momentum equation in radial direction r for the single velocity component u is given by

$\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}\right)=-\frac{\partial p}{\partial r}+\frac{\partial {\tau }_{rr}}{\partial r}+\frac{2}{r}\left({\tau }_{rr}-{\tau }_{\theta \theta }\right)$ , (8)

where p is the pressure, τ_rr and τ_θθ are the normal stresses in radial and polar direction, σ is the interface tension and ρ is the density of the liquid [9] . Inserting the continuity Equation (7) into the momentum Equation (8) and integrating the momentum equation within the limits r = r_B(t) and r = r_S(t) gives the non-linear, integro-differential equation of motion

$\begin{array}{l}\rho \left[r_B^2\left(\frac{1}{r_B}-\frac{1}{r_S}\right)\frac{{\text{d}}^2{r}_B}{\text{d}{t}^2}+\frac{1}{2}\left(\frac{r_B^4}{r_S^4}-4\frac{r_B}{r_S}+3\right){\left(\frac{\text{d}{r}_B}{\text{d}t}\right)}^2\right]\\ =p_g-\frac{2\sigma }{r_B}-p_{\infty }+2{\displaystyle \underset{r_B}{\overset{r_S}{\int }}\frac{1}{r}}\left({\tau }_{rr}-{\tau }_{\theta \theta }\right)\text{d}r\end{array}$ (9)

to be solved. Equation (9) assumes that due to the large density differences between the liquid and gaseous phase, no tangential momentum exchange occurs between both phases. The ambient pressure and pressure inside a bubble are interlinked by

$p_{\infty }\left(t\right)=p\left(r_S\right)-{\tau }_{rr}\left(r_S\right)$ , (10)

$p_g\left(t\right)=p\left(r_B\right)-{\tau }_{rr}\left(r_B\right)+\frac{2\sigma }{r_B}$ . (11)

The integral on the right side of Equation (9) characterizes the tensile stresses around the bubble. With regard to the rheological behavior of the liquid, Newtonian, i.e. purely viscous liquids are only considered here. The rheological equation of state for an incompressible Newtonian liquid is S = 2η₀D, where S is the extra stress tensor, D is the deformation velocity tensor and η₀ is the zero-shear viscosity of the liquid. The normal stress difference in (9) reads then

${\tau }_{rr}-{\tau }_{\theta \theta }=2{\eta }_0\left(\frac{\partial u}{\partial r}-\frac{u}{r}\right)=-6{\eta }_0\frac{r_B^2}{r^3}\cdot \frac{\text{d}{r}_B}{\text{d}t}$ . (12)

At time t = 0, the temperature ϑ₀ and concentrations ξ_i,0 of the volatile components dissolved in the liquid are postulated to be constant. For t > 0, temperature and concentration profiles arise in the liquid at the phase boundary. The mass transport of the volatile components from the liquid into a bubble is based on the assumptions that the ideal gas law is valid, the components do not influence each other and the concentrations ξ_i of the components are low so that the Henry law can be applied [9] . Hence, the concentration equation for each component i reads

$\frac{\partial {\xi }_i}{\partial t}+u\frac{\partial {\xi }_i}{\partial r}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\text{D}}_i\frac{\partial {\xi }_i}{\partial r}\right)$ , (13)

where D_i is the diffusion coefficient [9] . Due to the mass transport, a mass balance equation for each component i in the gas phase has to be taken into account [9] :

$\frac{\text{d}}{\text{d}t}\left(r_B^3{\rho }_i\right)=3\rho {\text{D}}_i{r}_B^2\cdot {\frac{\partial {\xi }_i}{\partial r}|}_{r=r_B}$ . (14)

By assuming a homogeneous temperature field inside the bubble, characterized by ϑ(r = r_B), the partial pressure p_i in the gas phase is calculated by $p_i={\rho }_i\cdot {\text{R}}_i\cdot \vartheta \left(r=r_B\right)$ , where R_i and ρ_i, respectively, are the specific gas constant and the density of the component i. The total pressure p_g inside the bubble is the sum of all partial pressures. According to Henry’s law, the concentration of a gas dissolved in a liquid is proportional to the pressure of the respective gas in the gaseous phase. The equilibrium gas concentration and thus the first boundary condition to solve the concentration equation have to be formulated adequately, that is ξ_GL,i = p_i/H_W,i, where H_W,i is the Henry coefficient of the component i. Due to the mass transport, a transport of thermal energy between the phases also takes place. Further simplified assumptions apply: (i) adiabatic foaming process, (ii) when considering a single bubble, pure one-dimensional heat conduction in radial direction, (iii) negligible production of internal energy by viscous dissipation or diffusion. Under these premises, the transport of thermal energy is expressed by

$\text{c}\rho \left(\frac{\partial \vartheta }{\partial t}+u\frac{\partial \vartheta }{\partial r}\right)=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\lambda }_{\text{W}}\frac{\partial \vartheta }{\partial r}\right)$ , (15)

where λ_W is the heat conductivity and c is the specific thermal capacity of the liquid [9] . Solving the temperature Equation (15) requires knowing the temperature gradient at the phase boundary as the first boundary condition, i.e.

${\lambda }_{\text{W}}\cdot {\frac{\partial \vartheta }{\partial r}|}_{r=r_B}=\rho {\displaystyle \underset{i}{\sum }{\text{D}}_i}\Delta {\text{h}}_{\text{S},i}\cdot {\frac{\partial {\xi }_i}{\partial r}|}_{r=r_B}+\frac{1}{r_B^2}\frac{\text{d}}{\text{d}t}\left(\sigma \cdot r_B^2\right)$ , (16)

where Δh_S,i is the absorption enthalpy of the component i [9] . Equation (16) is completed by work, which is performed to enlarge the bubble surface.

The validity range of the spherical foam model ends if the gas volume fraction exceeds a limit value of about 0.7. In the second model, the bubbles are assumed to have the idealized shape of a pentagon dodecahedron. The transformations at time t*, where the model change occurs, are subject to the following conditions [9] :

(i) The gas volume and liquid volume remain unchanged. (ii) The time derivative of the volume of a bubble is maintained. (iii) The partial pressures in the gaseous phase remain unchanged. (iv) The concentration and temperature profiles are transformed with respect to the constant liquid volume. According to the first condition (i), the edge length a_D,begin2 of a pentagon dodecahedron bubble at the beginning of the second model is calculated from the radius r_B,end1 of a spherical bubble at the end of the first model:

$a_{D,\text{begin2}}={\left(\frac{16\text{π}}{3\left(15+7\sqrt{5}\right)}\right)}^{1/3}{r}_{B,\text{end1}}$ . (17)

The liquid volume ${\text{V}}_{\text{f}}^{\text{e}}/4$ is transformed into a thin layer normal to the bubble surface with thickness

$h_{\text{begin2}}=\frac{{\text{V}}_{\text{f}}^{\text{e}}/4}{3\sqrt{25+10\sqrt{5}}\cdot a_{D,\text{begin2}}^2}$ . (18)

The second condition (ii) for the model change leads to the time derivatives

$\frac{\text{d}{a}_D}{\text{d}t}=-\frac{1}{2}\cdot \frac{a_D}{h}\cdot \frac{\text{d}h}{\text{d}t}$ , (19)

$\frac{\text{d}{h}_{\text{begin2}}}{\text{d}t}=-2\cdot \frac{h_{\text{begin2}}}{r_{B,\text{end1}}}\cdot \frac{\text{d}{r}_{B,\text{end1}}}{\text{d}t}$ . (20)

2.2. Model Equations for the Polyhedral Foam

In the second model, the bubbles are assumed to have the idealized shape of a pentagon dodecahedron. The liquid volume is transformed into a thin layer normal to the surface of the bubbles. According to the model assumption, two neighboring bubbles share a liquid lamella, i.e. the calculations are limited to one half of the liquid lamella. The conservation equations for a liquid lamella are formulated in cylindrical coordinates (r, φ, z) with the velocities u in radial direction and w in axial direction. The geometrical formulation in cylindrical coordinates is based on the assumption that the stress state within a side surface of a pentagon dodecahedron bubble does not differ significantly from the stress state within a circular side surface of same area [9] . Figure 2 depicts schematically a liquid layer with height h(t), which is surrounded by the Plateau edges, the gas phase of a dodecahedron bubble and the liquid layer of the neighboring dodecahedron bubble.

The following additional model assumptions apply: In the selected cylindrical coordinate system, a biaxial elongational flow is formed, where the liquid lamella is shortened in axial direction and stretched in the other two directions. The normal stresses τ_rr and τ_φφ have the same size. Interface tension forces do not occur due the negligible curvature of the liquid lamella. Under these premises, the conservation equation for mass for the prevailing two-dimensional velocity field in the liquid phase (u, 0, w) reads (see also Figure 2)

$u\left(r,t\right)=-\frac{1}{2}\cdot \frac{r}{h}\cdot \frac{\text{d}h}{\text{d}t}=-\frac{r}{2}\stackrel{\dot{}}{\epsilon }\left(t\right)$ (21)

with the elongational rate

$\frac{\partial w\left(z,t\right)}{\partial z}=\frac{1}{h}\frac{\text{d}h}{\text{d}t}=\stackrel{\dot{}}{\epsilon }\left(t\right)$ . (22)

The momentum equation in axial direction z simplifies to

$\rho \left(\frac{\partial w}{\partial t}+w\frac{\partial w}{\partial z}\right)=-\frac{\partial p}{\partial z}+\frac{\partial {\tau }_{zz}}{\partial z}$ . (23)

Inserting the relation (22) into the momentum Equation (23) and integrating the momentum equation within the limits z = 0 and z = h(t) gives the equation of motion to be solved. The balance of forces at the phase boundary to the bubble z = h(t) and at the Plateau edges of the liquid lamella z = 0 are taken into account so that

$\rho \frac{h}{2}\frac{{\text{d}}^{2}h}{\text{d}{t}^2}=-p_g+p_{\infty }-\left[{\tau }_{zz}\left(0\right)-{\tau }_{rr}\left(0\right)\right]$ (24)

in accordance to [9] . For a Newtonian liquid with zero-shear viscosity η₀, the normal stress difference in (24) reads

${\tau }_{zz}-{\tau }_{rr}=3{\eta }_0\stackrel{\dot{}}{\epsilon }$ . (25)

The assumptions already made in the first model for the mass and energy transport also apply here in the second model [9] . After the transformation $x=h-z$ , the concentration and mass balance equation for each component i are given by

$\frac{\partial {\xi }_i}{\partial t}+\frac{\text{d}x}{\text{d}t}\cdot \frac{\partial {\xi }_i}{\partial x}=\frac{\partial }{\partial x}\left({\text{D}}_i\frac{\partial {\xi }_i}{\partial x}\right)$ , (26)

$\frac{\text{d}}{\text{d}t}\left({\rho }_i{\text{V}}_{\text{D}}\right)=\rho {\text{D}}_i{\text{O}}_{\text{D}}\cdot {\frac{\partial {\xi }_i}{\partial x}|}_{x=0}$ . (27)

The temperature equation and first boundary condition are given by

$\frac{\partial \vartheta }{\partial t}+\frac{\text{d}x}{\text{d}t}\cdot \frac{\partial \vartheta }{\partial x}=\frac{{\lambda }_{\text{W}}}{\text{c}\rho }\cdot \frac{{\partial }^2\vartheta }{\partial x^2}$ , (28)

${\lambda }_{\text{W}}\cdot {\frac{\partial \vartheta }{\partial x}|}_{x=0}=\rho {\displaystyle \underset{i}{\sum }{\text{D}}_i}\Delta {\text{h}}_{\text{S},i}\cdot {\frac{\partial {\xi }_i}{\partial x}|}_{x=0}$ . (29)

2.3. Model Equations in the Dimensionless Form

The model Equations (9)-(16) for the spherical foam are subjected to the Lagrange coordinate transformation $y=\left(r^3-r_B^3\right)/3$ and transformed into dimensionless equations [9] . Similarly, the model Equations (24)-(29) for the polyhedral foam are subjected to the Lagrange coordinate transformation $q=x{a}_D^2$ and transformed into dimensionless equations [9] . As shown in [9] , the coordinate transformations ensure that the convective terms in the concentration and temperature Equations (13), (15), (26) and (28) are eliminated. From a dimensional analysis, a set of dimensionless π-numbers is obtained:

${\pi }_1=\frac{\sigma }{r_{B,0}{p}_{\infty ,0}},{\pi }_3=\frac{{\eta }_0}{\rho r_{B,0}}{\left(\frac{\rho }{p_{\infty ,0}}\right)}^{0.5}$ , (30)

${\pi }_{4,i}=\frac{{\text{D}}_{i,0}}{r_{B,0}}{\left(\frac{\rho }{p_{\infty ,0}}\right)}^{0.5},{\pi }_5=\frac{{\lambda }_{\text{W}}}{\text{c}\rho r_{B,0}}{\left(\frac{\rho }{p_{\infty ,0}}\right)}^{0.5},{\pi }_6={\vartheta }_0\text{c}\frac{\rho }{p_{\infty ,0}}$ , (31)

${\pi }_{7,i}={\text{R}}_i{\vartheta }_0\frac{\rho }{p_{\infty ,0}}$ . (32)

The index 0 refers in each case to the time point t = 0, whereas the index i refers to a volatile component i dissolved in the liquid. The dimensionless quantities are indicated by the corresponding capital letters.

The relevant quantities are the time $T=\left(t/r_{B,0}\right)\cdot {\left({\rho }_{\infty ,0}/\rho \right)}^{0.5}$ , the radius of a spherical bubble $R_B=\left(r_B/r_{B,0}\right)$ , the outer radius of the liquid volume that surrounds a bubble $R_S=\left(r_S/r_{S,0}\right)$ , the transformed spatial coordinate $Y=\left(y/r_{B,0}^3\right)$ , the normal stresses $S_{YY}=\left({\tau }_{yy}/{\rho }_{\infty ,0}\right)$ and $S_{\theta \theta }=\left({\tau }_{\theta \theta }/{\rho }_{\infty ,0}\right)$ , the temperature $\Theta =\left(\vartheta /{\vartheta }_0\right)$ , the pressure inside a bubble $P_g=\left(p_g/p_{\infty ,0}\right)$ and the ambient pressure $P_{\infty }=\left(p_{\infty }/p_{\infty ,0}\right)$ . The additional dimensionless quantities used for the polyhedral foam model are the thickness of a liquid lamella $H=\left(h/r_{B,0}\right)$ , the edge length of a dodecahedron bubble $A_D=\left(a_D/r_{B,0}\right)$ , the transformed spatial coordinate $Q=\left(q/r_{B,0}^3\right)$ and the corresponding normal stresses $S_{ZZ}=\left({\tau }_{zz}/{\rho }_{\infty ,0}\right)$ and $S_{RR}=\left({\tau }_{rr}/{\rho }_{\infty ,0}\right)$ .

The model equations for the spherical and polyhedral foam in the dimensionless form with appropriate initial and boundary conditions [9] [10] are now specified.

Spherical foam:

Equation of motion at the phase boundary

$\begin{array}{l}{R}_B^2\left(\frac{1}{R_B}-\frac{1}{R_S}\right)\frac{{\text{d}}^2{R}_B}{\text{d}{T}^2}+\frac{1}{2}\left(\frac{R_B^4}{R_S^4}-4\frac{R_B}{R_S}+3\right){\left(\frac{\text{d}{R}_B}{\text{d}T}\right)}^2\\ =P_g-\frac{2{\pi }_1}{R_B}-P_{\infty }+2{\displaystyle \underset{0}{\overset{Y_S}{\int }}\frac{S_{YY}-S_{\theta \theta }}{3Y+R_B^3}}\text{d}Y\end{array}$ (33)

$R_B\left(T=0\right)=1$ and $\frac{\text{d}{R}_B}{\text{d}T}\left(T=0\right)=0$ , where $Y_S=\frac{1}{3}\left(R_S^3-R_B^3\right)$ .

Newtonian liquid

$2{\displaystyle \underset{0}{\overset{Y_S}{\int }}\frac{S_{YY}-S_{\theta \theta }}{3Y+R_B^3}}\text{d}Y=4{\pi }_3\frac{1}{R_B}\left(\frac{R_B^3}{R_S^3}-1\right)\frac{\text{d}{R}_B}{\text{d}T}$ (34)

Concentration equations

$\frac{\partial {\xi }_i}{\partial T}={\pi }_{4,i}\frac{\partial }{\partial Y}\left[{\left(3Y+R_B^3\right)}^{4/3}\frac{{\text{D}}_i}{{\text{D}}_{i,0}}\frac{\partial {\xi }_i}{\partial Y}\right]$ (35)

$T=0,Y=0:{\xi }_{i=1}={\xi }_{\text{GL},i=1},{\xi }_{i>1}=0$ and $T=0,0<Y\le Y_S:{\xi }_i={\xi }_{i,0}$ ,

$T>0,Y=0:{\xi }_i={\xi }_{\text{GL},i}$ and $T>0,Y=Y_S:\frac{\partial {\xi }_i}{\partial Y}=0$ .

Mass balance equations

$\frac{\text{d}}{\text{d}T}\left(R_B^3\frac{{\rho }_i}{\rho }\right)=3\frac{{\text{D}}_i}{{\text{D}}_{i,0}}{\pi }_{4,i}{R}_B^4\cdot {\frac{\partial {\xi }_i}{\partial Y}|}_{Y=0}$ (36)

$\frac{{\rho }_{i=1,0}}{\rho }=\frac{1+2{\pi }_1}{{\pi }_{7,1}}$ and ${\rho }_{i>1,0}=0$ .

Temperature equation

$\frac{\partial \Theta }{\partial T}={\pi }_5\frac{\partial }{\partial Y}\left[{\left(3Y+R_B^3\right)}^{4/3}\frac{\partial \Theta }{\partial Y}\right]$ (37)

${\pi }_5{\pi }_6\cdot {\frac{\partial \Theta }{\partial Y}|}_{Y=0}={\displaystyle \underset{i}{\sum }\frac{{\text{D}}_i}{{\text{D}}_{i,0}}}{\pi }_{4,i}\frac{\rho \Delta {\text{h}}_{\text{S},i}}{p_{\infty ,0}}\cdot {\frac{\partial {\xi }_i}{\partial Y}|}_{Y=0}+\frac{2{\pi }_1}{R_B^3}\frac{\text{d}{R}_B}{\text{d}T}$ (38)

$T=0,0\le Y\le Y_S:\Theta =1$ ,

$T>0,Y=0:\frac{\partial \Theta }{\partial Y}$ (Equation (38)) and $T>0,Y=Y_S:\frac{\partial \Theta }{\partial Y}=0$ .

The pressure inside a bubble is determined by $P_g=\left({\displaystyle \underset{i}{\sum }{\pi }_{7,i}\frac{{\rho }_i}{\rho }}\right)\cdot {\Theta |}_{Y=0}$ .

Polyhedral foam:

Equation of motion at the phase boundary

$\frac{H}{2}\frac{{\text{d}}^{2}H}{\text{d}{T}^2}=-P_g+P_{\infty }-\left(S_{ZZ}-S_{RR}\right)$ (39)

$H\left(T=T^{*}\right)=H_{\text{begin2}}$ and $\frac{\text{d}H}{\text{d}T}\left(T=T^{*}\right)=\frac{\text{d}{H}_{\text{begin2}}}{\text{d}T}$ .

Newtonian liquid

$S_{ZZ}-S_{RR}=3{\pi }_3\frac{1}{H}\frac{\text{d}H}{\text{d}T}$ (40)

Concentration equations

$\frac{\partial {\xi }_i}{\partial T}={\pi }_{4,i}{A}_D^4\frac{\partial }{\partial Q}\left[\frac{{\text{D}}_i}{{\text{D}}_{i,0}}\frac{\partial {\xi }_i}{\partial Q}\right]$ (41)

$T=T^{*},0\le Q\le H{A}_D^2:{\xi }_i={\xi }_{i,\text{end1}}$ ,

$T>T^{*},Q=0:{\xi }_i={\xi }_{\text{GL},i}$ and $T>T^{*},Q=H{A}_D^2:\frac{\partial {\xi }_i}{\partial Q}=0$ .

Mass balance equations

$\frac{\text{d}}{\text{d}T}\left(A_D^3\frac{{\rho }_i}{\rho }\right)=12\frac{\sqrt{25+10\sqrt{5}}}{15+7\sqrt{5}}\frac{{\text{D}}_i}{{\text{D}}_{i,0}}{\pi }_{4,i}{A}_D^4\cdot {\frac{\partial {\xi }_i}{\partial Q}|}_{Q=0}$ (42)

Temperature equation

$\frac{\partial \Theta }{\partial T}={\pi }_5{A}_D^4\frac{{\partial }^2\Theta }{\partial Q^2}$ (43)

${\pi }_5{\pi }_6\cdot {\frac{\partial \Theta }{\partial Q}|}_{Q=0}={\displaystyle \underset{i}{\sum }\frac{{\text{D}}_i}{{\text{D}}_{i,0}}}{\pi }_{4,i}\frac{\rho \Delta {\text{h}}_{\text{S},i}}{p_{\infty ,0}}\cdot {\frac{\partial {\xi }_i}{\partial Q}|}_{Q=0}$ (44)

$T=T^{*},0\le Q\le H{A}_D^2:\Theta ={\Theta }_{\text{end1}}$ ,

$T>T^{*},Q=0:\frac{\partial \Theta }{\partial Q}$ (Equation (44)) and $T>T^{*},Q=H{A}_D^2:\frac{\partial \Theta }{\partial Q}=0$ .

3. Numerical Method for Solving the Model Equations

The models developed in [9] [10] have been implemented with the software Matlab. The integro-differential Equations (33)-(38) and differential Equations (39)-(44) present two non-linear systems of equations that are treated consecutively. To solve the equations numerically, a finite difference discretization method and a diagonally implicit Runge-Kutta method are applied. The discretized equations are solved separately from each other. The concentration and mass balance equations, the temperature equation and the equation of motion are solved in this order for each model. Since the equations are non-linear and coupled, they have to be solved iteratively in each time step. For the outer iteration process, a Picard iteration method (fixed-point method with successive substitution) is used. In this section, discretization and solution of the integro-differential Equations (33)-(38) in a time step for the first model are explained. The same considerations apply to solving the differential Equations (39)-(44) for the second model numerically.

The concentration and temperature Equations (35) and (37) are parabolic partial differential equations. They are discretized using the finite difference method. As an example, the concentration Equation (35) for a gaseous component i is studied. The diffusion coefficient is assumed to be constant, that is D_i = const. Applying the product rule gives the equation

$\frac{\partial {\xi }_i}{\partial T}={\pi }_{4,i}\left[4{\left(3Y+R_B^3\right)}^{1/3}\frac{\partial {\xi }_i}{\partial Y}+{\left(3Y+R_B^3\right)}^{4/3}\frac{{\partial }^2{\xi }_i}{\partial Y^2}\right]$ . (45)

Equation (45) is discretized in $\Omega \times \Gamma$ with $\Omega =\left[0,Y_S\right]$ and $\Gamma =\left[0,T_{\text{end1}}\right]$ , where discretization is done first in space and then in time. Note that, in the following, the indices n and k refer to discrete points in space and time. Note also that, in section 3 and section 4.1, the constant discretization parameter h denotes the spatial step size, not the thickness of a liquid lamella in the polyhedral foam. With the constant discretization parameter $h=Y_S/\text{N}$ , $\text{N}\in \mathbb{N}$ , the set of all spatial grid points is defined: ${\Omega }_h:=\left\{Y_n\in ℝ|Y_n=n\cdot h\text{and}0\le n\le \text{N}\right\}$ . The discretization of $\Omega =\left[0,Y_S\right]$ with the constant discretization parameter h is always equidistant.

The spatial derivatives in (45) are approximated by central difference quotients. Inserting the central difference quotients in (45) gives a coupled system of ordinary differential equations for the unknown functions of the gaseous components dissolved in the liquid, ${\xi }_i\left(Y_n,T\right)$ , $1\le n\le \text{N}-1$ , that is

$\begin{array}{c}\frac{\partial {\xi }_i\left(Y_n,T\right)}{\partial T}={\pi }_{4,i}{\left(3Y_n+R_B^3\right)}^{1/3}\cdot \frac{{\xi }_i\left(Y_{n+1},T\right)-{\xi }_i\left(Y_{n-1},T\right)}{2h}\cdot 4\\ +{\pi }_{4,i}{\left(3Y_n+R_B^3\right)}^{4/3}\cdot \frac{{\xi }_i\left(Y_{n+1},T\right)-2{\xi }_i\left(Y_n,T\right)+{\xi }_i\left(Y_{n-1},T\right)}{h^2}\\ +{\pi }_{4,i}\cdot O\left(h^2\right)\end{array}$ (46)

For the time discretization, the grid ${\Gamma }_k:=\left\{T_k\in ℝ|T_k=k\cdot v\text{and}0\le k\le \text{K}\right\}$ is selected, where the discretization parameter $v=T_{\text{end1}}/\text{K}$ , $\text{K}\in \mathbb{N}$ , is the temporal step size. The system (46) is integrated by one-step-θ methods [16] [17] . The parameter $\theta \in \left[0,1\right]$ determines the concrete difference scheme. If the right side of Equation (46) is referred to as $F_i\left(Y_n,T\right)$ , the integration within the limits $T=T_k$ and $T=T_{k+1}$ gives

${\displaystyle \underset{T_k}{\overset{T_{k+1}}{\int }}\frac{\partial {\xi }_i\left(Y_n,T\right)}{\partial T}}\text{d}T={\displaystyle \underset{T_k}{\overset{T_{k+1}}{\int }}{F}_i\left(Y_n,T\right)\text{d}T}$

$\begin{array}{c}\iff {\xi }_i\left(Y_n,T_{k+1}\right)-{\xi }_i\left(Y_n,T_k\right)=v\cdot \left[\left(1-\theta \right)\cdot F_i\left(Y_n,T_k\right)+\theta \cdot F_i\left(Y_n,T_{k+1}\right)\right]\\ +{\pi }_{4,i}\cdot O\left(v^2\right)\end{array}$ (47)

For $\theta \in \left[0,1\right]$ , the Equations (47) have always accuracy order 2, since R_B(T_{k + 1}) is not known yet and R_B(T_k) or an approximated value obtained after an iteration step has to be inserted. The following abbreviations are now defined:

$F_{i,n}^k:={\pi }_{4,i}{\left(3Y_n+R_B^3\left(T_k\right)\right)}^{1/3}\cdot 2h$ , (48)

$E_{i,n}^k:={\pi }_{4,i}{\left(3Y_n+R_B^3\left(T_k\right)\right)}^{4/3}$ , (49)

$\gamma :=\frac{v}{h^2}$ . (50)

Inserting the expressions for $F_i\left(Y_n,T_k\right)$ and $F_i\left(Y_n,T_k{}_{+1}\right)$ in (47) gives the discretized equations for the variables ${\xi }_i\left(Y_n,T_k\right),\text{1}\le n\le \text{N}-\text{1}$ , that is

$\begin{array}{l}\left(1+2\theta \gamma E_{i,n}^k\right)\cdot {\xi }_i\left(Y_n,T_{k+1}\right)-\theta \gamma \left(-F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_i\left(Y_{n-1},T_{k+1}\right)\\ -\theta \gamma \left(F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_i\left(Y_{n+1},T_{k+1}\right)\\ =\left(1-2\left(1-\theta \right)\gamma E_{i,n}^k\right)\cdot {\xi }_i\left(Y_n,T_k\right)+\left(1-\theta \right)\gamma \left(-F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_i\left(Y_{n-1},T_k\right)\\ +\left(1-\theta \right)\gamma \left(F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_i\left(Y_{n+1},T_k\right)+{\pi }_{4,i}\cdot O\left(v{h}^2+v^2\right)\end{array}$ (51)

For each discrete time point T_k, the grid function ${\xi }_{i,h}^k:{\Omega }_h\to ℝ$ is defined. The equations for the grid function are obtained when the error term in (51) is omitted:

$\begin{array}{l}\left(1+2\theta \gamma E_{i,n}^k\right)\cdot {\xi }_{i,h}^{k+1}\left(Y_n\right)-\theta \gamma \left(-F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_{i,h}^{k+1}\left(Y_{n-1}\right)\\ -\theta \gamma \left(F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_{i,h}^{k+1}\left(Y_{n+1}\right)\\ =\left(1-2\left(1-\theta \right)\gamma E_{i,n}^k\right)\cdot {\xi }_{i,h}^k\left(Y_n\right)+\left(1-\theta \right)\gamma \left(-F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_{i,h}^k\left(Y_{n-1}\right)\\ +\left(1-\theta \right)\gamma \left(F_{i,n}^k+E_{i,n}^k\right)\cdot {\xi }_{i,h}^k\left(Y_{n+1}\right)\end{array}$ (52)

Hence, the discrete solution is determined by a linear system of equations

$L_{i,h}^k{\xi }_{i,h}^{k+1}=M_{i,h}^k{\xi }_{i,h}^k+b_{i,h}^k$ , (53)

$L_{i,h}^k,M_{i,h}^k\in {ℝ}^{\left(\text{N}-\text{1}\right)\times \left(\text{N}-\text{1}\right)},b_{i,h}^k,{\xi }_{i,h}^k\in {ℝ}^{\text{N}-\text{1}}$ ,

that has to be solved in each time step. The initial and boundary conditions are taken into account by the vector $b_{i,h}^k$ , the first and the last component of which are only non-zero. The matrix $L_{i,h}^k$ is diagonal if θ = 0 and tridiagonal otherwise. The following one-step-θ methods are often used in practice: explicit Euler method of order 1 (θ = 0), implicit Euler method of order 1 (θ = 1) and, finally, Crank-Nicolson method of order 2 (θ = 0.5) [16] [17] . If I_h is the identity matrix associated with the discretization, the matrix $L_{i,h}^k$ (similarly $M_{i,h}^k$ ) has the form:

$\begin{array}{l}{L}_{i,h}^k=I_h\\ +\theta \gamma \cdot \left[\begin{array}{ccccc}2E_{i,1}^k& -F_{i,1}^k-E_{i,1}^k& 0& \cdots & 0\\ F_{i,2}^k-E_{i,2}^k& 2E_{i,2}^k& -F_{i,2}^k-E_{i,2}^k& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ 0& \ddots & F_{i,\text{N}-2}^k-E_{i,\text{N}-2}^k& 2E_{i,\text{N}-2}^k& -F_{i,\text{N}-2}^k-E_{i,\text{N}-2}^k\\ 0& \cdots & 0& F_{i,\text{N}-1}^k-E_{i,\text{N}-1}^k& 2E_{i,\text{N}-1}^k\end{array}\right]\end{array}$ (54)

The numerical analysis is based on the properties of M-matrices. In the appendix, relevant definitions and theorems relating to M-matrices [14] [15] [16] [17] are summarized. Since $\theta \in \left[0,\text{1}\right]$ , $\gamma >0$ and $0<F_{i,n}^k<E_{i,n}^k$ is satisfied for all spatial and temporal grid points if $0<h<0.\text{5}$ in (48)-(50), following properties can be shown:

1) The matrix $L_{i,h}^k$ satisfies the sign condition in definition 1.

2) The matrix $L_{i,h}^k$ is diagonally dominant and irreducible.

According to theorem 1, $L_{i,h}^k$ is regular, according to theorem 2, $L_{i,h}^k$ is even an M-matrix. Hence, the linear system of Equations (53) is uniquely solvable. For a vector $w=1$ , it can be shown that $L_{i,h}^{k}w\ge 1$ , i.e.

${\left(L_{i,h}^{k}w\right)}_n=1+2\theta \gamma E_{i,n}^k-\theta \gamma \left(-F_{i,n}^k+E_{i,n}^k\right)-\theta \gamma \left(F_{i,n}^k+E_{i,n}^k\right)=1$ for $\text{2}\le n\le \text{N}-\text{2}$ ,

${\left(L_{i,h}^{k}w\right)}_1=1+\theta \gamma E_{i,1}^k-\theta \gamma F_{i,1}^k\ge 1$ for $n=\text{1}$ ,

${\left(L_{i,h}^{k}w\right)}_{N-1}=1+\theta \gamma E_{i,N-1}^k+\theta \gamma F_{i,N-1}^k\ge 1$ for $n=\text{N}-\text{1}$ .

Hence, theorem 3 is applicable, and the following estimation is valid:

${\Vert {\left(L_{i,h}^k\right)}^{-1}\Vert }_{\infty }\le 1$ . (55)

Estimation (55) is relevant for the convergence considerations in the next section. The Equation (46), $1\le n\le \text{N}-1$ which is discretized in space, represents a stiff system of ordinary differential equations. The greatest eigenvalue compared to the smallest eigenvalue behaves like O(h⁻²). Therefore, A-stable numerical time integration procedures are required [16] [17] .

The temperature Equation (37) is discretized and solved similarly to the concentration Equation (35). The mass balance Equation (36) is integrated by the one-step-θ method, which is also used for the concentration Equation (35). The equation of motion (33) is transformed into a system of first order integro-differential equations. Since stiffness occurs, this system is solved by the TR-BDF2 method implemented in Matlab. The TR-BDF2 method is a diagonally implicit Runge-Kutta method of order 2 with associated error estimate of order 3 and is L-stable and also strongly S-stable [11] [12] [13] . For an incompressible Newtonian liquid, the integral in (33) can be calculated analytically so that the system depends only on time (see Equation (34)). The equation of motion (33) requires particularly a high computing time if a low-viscosity Newtonian liquid is given so that a small temporal step size has to be chosen.

4. Presentation and Discussion of Results

4.1. Considerations on Convergence

This section aims at proven convergence of the finite difference discretization method and its order [11] - [17] . The convergence considerations concern the transport processes of mass and thermal energy from the liquid in the gas phase. The convergence theory is based on the properties of the class of M-matrices and follows the numerical studies [16] [17] . It is stated that the outer iteration method gives convergent approximated values. The error function at time T_k is defined by

$e_{i,h}^k=R_h{\xi }_i\left(\cdot ,T_k\right)-{\xi }_{i,h}^k$ . (56)

R_h denotes the restriction operator that gives the exact solution of (45) in the spatial grid points at time T_k, whereas ${\xi }_{i,h}^k$ is the approximated solution. From (51)-(52), it follows that the finite difference method is consistent. The maximum norm of the local truncation error ${\eta }_{i,h}^k$ is

${\Vert {\eta }_{i,h}^k\Vert }_{\infty }={\pi }_{4,i}\cdot O\left(v{h}^2+v^2\right)$ . (57)

Now, the global discretization error is studied. Note that, from (51)-(52), it follows

$L_{i,h}^k{R}_h{\xi }_i\left(\cdot ,T_{k+1}\right)-M_{i,h}^k{R}_h{\xi }_i\left(\cdot ,T_k\right)-b_{i,h}^k={\eta }_{i,h}^k$ . (58)

Applying $L_{i,h}^k$ to $e_{i,h}^{k+1}$ gives

$\begin{array}{c}{L}_{i,h}^k{e}_{i,h}^{k+1}=L_{i,h}^k{R}_h{\xi }_i\left(\cdot ,T_{k+1}\right)-L_{i,h}^k{\xi }_{i,h}^{k+1}\\ =M_{i,h}^k{R}_h{\xi }_i\left(\cdot ,T_k\right)+b_{i,h}^k+{\eta }_{i,h}^k-M_{i,h}^k{\xi }_{i,h}^k-b_{i,h}^k\\ =M_{i,h}^k{e}_{i,h}^k+{\eta }_{i,h}^k\end{array}$ (59)

Hence, a recursion formula for the error function is obtained. Multiplying Equation (59) from left by the inverse matrix ${\left(L_{i,h}^k\right)}^{-1}$ and applying the maximum norm gives

${\Vert e_{i,h}^{k+1}\Vert }_{\infty }\le {\Vert {\left(L_{i,h}^k\right)}^{-1}\Vert }_{\infty }{\Vert M_{i,h}^k\Vert }_{\infty }{\Vert e_{i,h}^k\Vert }_{\infty }+{\Vert {\left(L_{i,h}^k\right)}^{-1}\Vert }_{\infty }{\Vert {\eta }_{i,h}^k\Vert }_{\infty }$ . (60)

By recursive application, the following in equation is obtained:

$\begin{array}{c}{\Vert e_{i,h}^k\Vert }_{\infty }\le \left\{{\displaystyle \underset{l=0}{\overset{k-1}{\prod }}{\Vert {\left(L_{i,h}^l\right)}^{-1}\Vert }_{\infty }{\Vert M_{i,h}^l\Vert }_{\infty }}\right\}{\Vert e_{i,h}^0\Vert }_{\infty }\\ +{\displaystyle \underset{l=0}{\overset{k-1}{\sum }}\left\{{\displaystyle \underset{j=l}{\overset{k-1}{\prod }}{\Vert {\left(L_{i,h}^j\right)}^{-1}\Vert }_{\infty }{\displaystyle \underset{j=l+1}{\overset{k-1}{\prod }}{\Vert M_{i,h}^j\Vert }_{\infty }}}\right\}}{\Vert {\eta }_{i,h}^l\Vert }_{\infty }\end{array}$ (61)

The matrix norm of $M_{i,h}^k$ can be determined:

$\begin{array}{c}{\Vert M_{i,h}^k\Vert }_{\infty }=ma{x}_n\{\left|1-2\left(1-\theta \right)\gamma E_{i,n}^k\right|+\left|\left(1-\theta \right)\gamma \left(-F_{i,n}^k+E_{i,n}^k\right)\right|\\ +\left|\left(1-\theta \right)\gamma \left(F_{i,n}^k+E_{i,n}^k\right)\right|,\left|1-2\left(1-\theta \right)\gamma E_{i,1}^k\right|+\left|\left(1-\theta \right)\gamma \left(F_{i,1}^k+E_{i,1}^k\right)\right|,\\ \left|1-2\left(1-\theta \right)\gamma E_{i,\text{N}-1}^k\right|+\left|\left(1-\theta \right)\gamma \left(-F_{i,\text{N}-1}^k+E_{i,\text{N}-1}^k\right)\right|\}\\ =ma{x}_n\left\{\left|1-2\left(1-\theta \right)\gamma E_{i,n}^k\right|+2\left(1-\theta \right)\gamma E_{i,n}^k\right\}=1,\end{array}$ (62)

hence ${\Vert M_{i,h}^k\Vert }_{\infty }=1$ if

$2\left(1-\theta \right)\frac{v}{h^2}{E}_{i,n}^k\le 1\iff v\le \frac{h^2}{2\left(1-\theta \right)E_{i,n}^k}=\{\begin{array}{ll}\infty ,\hfill & \theta =1\hfill \\ h^2/E_{i,n}^k,\hfill & \theta =0.5\hfill \\ h^2/\left(2E_{i,n}^k\right),\hfill & \theta =0\hfill \end{array}$ (63)

Using (55) and (62) in (61), it follows

$\begin{array}{c}{\Vert e_{i,h}^k\Vert }_{\infty }\le \underset{\le 1}{\underset{︸}{{\displaystyle \underset{l=0}{\overset{k-1}{\prod }}{\Vert {\left(L_{i,h}^l\right)}^{-1}\Vert }_{\infty }{\Vert M_{i,h}^l\Vert }_{\infty }}}}\cdot {\Vert e_{i,h}^0\Vert }_{\infty }\\ +{\displaystyle \underset{l=0}{\overset{k-1}{\sum }}\underset{\le 1}{\underset{︸}{{\displaystyle \underset{j=l}{\overset{k-1}{\prod }}{\Vert {\left(L_{i,h}^j\right)}^{-1}\Vert }_{\infty }{\displaystyle \underset{j=l+1}{\overset{k-1}{\prod }}{\Vert M_{i,h}^j\Vert }_{\infty }}}}}}\cdot {\Vert {\eta }_{i,h}^l\Vert }_{\infty }\\ \le {\Vert e_{i,h}^0\Vert }_{\infty }+k\cdot ma{x}_l{\Vert {\eta }_{i,h}^l\Vert }_{\infty }\\ ={\Vert e_{i,h}^0\Vert }_{\infty }+{\pi }_{4,i}\cdot O\left(h^2+v\right)\text{dueto}k=\frac{T_k}{v}.\end{array}$ (64)

The finite difference method is consistent and stable and thus convergent. The

stability results from the fact that $L_{i,h}^k$ are M-matrices and ${\Vert {\left(L_{i,h}^k\right)}^{-1}\Vert }_{\infty }\le 1$ for

each discrete time point. The convergence analysis also shows that the finite difference method has convergence order 2 in space and at least convergence order 1 in time. Analogous studies can be carried out for the concentration and temperature equations in the second model, where similar results are obtained. In the numerical simulations, the implicit Euler method is used because this method is A-stable and it is not subject to the restriction (63).

4.2. Results of the Numerical Simulations

As a reference example, the foaming process of a two-phase system consisting of polystyrene, styrene and nitrogen is considered. The nitrogen was dispersed into polystyrene under pressure and is in form of small spherical bubbles, but also in dissolved form in polystyrene, whereas styrene is already dissolved in polystyrene. The necessary material and process data are taken from the literature [9] [10] , see Table 1 and Table 2. The time process of the pressure drop is given by

$P_{\infty }=\frac{p_{\infty }\left(t\right)}{p_{\infty ,\text{max}}}=e^{-\left({\text{K}}_{\text{exp}}\right)t}$ , (65)

where K_exp = 0.6 s⁻¹ is the pressure drop coefficient. The coefficient K_exp is chosen so that the minimal pressure p_∞,min = 0.025 bar is reached from the maximal pressure p_∞,max = 10 bar after about 10 seconds. The heat conductivity and specific thermal capacity of polystyrene are specified by λ_W = 0.17 W/(m・K) and c = 1500 J/(kg・K). Note that, at time T = 0, the bubbles consist only of nitrogen. For T > 0, the gaseous components dissolved in polystyrene diffuse into the bubbles.

In summary, foam growth is induced by dispersion of nitrogen, vacuum application and diffusion of the volatile components dissolved in polystyrene.

The dimensionless foam volume V is determined on the basis of an elementary cell in the first model (see Equations (1)-(3)), i.e.

$\text{V}=\frac{{\text{V}}_{\text{ges}}^{\text{e}}}{{\text{V}}_{\text{ges},0}^{\text{e}}}=\frac{\frac{16}{3}\text{π}{r}_B^3+\frac{16}{3}\text{π}{r}_{B,0}^3\left({\text{c}}_{\text{V},0}^{-1}-1\right)}{\frac{16}{3}\text{π}\frac{r_{B,0}^3}{{\text{c}}_{\text{V},0}}}=1+{\text{c}}_{\text{V},0}\left(R_B^3-1\right)$ , (66)

and on the basis of a single bubble and the constant liquid volume around it in the second model (see Equation (5)), i.e.

$\text{V}=1+{\text{c}}_{\text{V},0}\left[\frac{3}{16\text{π}}\left(15+7\sqrt{5}\right)\cdot A_D^3-1\right]$ . (67)

The local considerations suggest the definition of a local Reynolds number that characterizes the foaming process. The Reynolds number is defined by $\text{Re}=\text{q}/\left(\nu \cdot l\right)$ , where q is the volume flow rate, $\nu =\eta /\rho$ is the kinematic viscosity and l is the length of the relevant flow region [18] . Selecting the volume flow rate, when the model changes from spherical to polyhedral foam, gives the local Reynolds number

$\text{Re}=\frac{16\text{π}{r}_{B,\text{max}}^2\cdot \frac{\text{d}{r}_{B,\text{max}}}{\text{d}t}}{\nu {\left({\text{V}}_{\text{S}}\right)}^{1/3}}=\frac{16\text{π}{R}_{B,\text{max}}^2\cdot \frac{\text{d}{R}_{B,\text{max}}}{\text{d}T}{\left(\frac{p_{\infty ,0}}{\rho }\right)}^{0.5}{r}_{B,0}}{\nu {\left(\frac{16}{3}\text{π}\left(\frac{\text{π}}{3\sqrt{2}}\frac{1}{{\text{c}}_{\text{V},0}}-1\right)\right)}^{1/3}}$ . (68)

The volume flow rate q is determined by the time derivative of the foam volume ${\text{V}}_{\text{ges}}^{\text{e}}$ in Equation (66). The maximal bubble radius can be calculated analytically,

$R_{B,\text{max}}={\left(\frac{{\text{c}}_{\text{V},0}^{-1}-1}{{\text{c}}_{\text{V},t*}^{-1}-1}\right)}^{1/3}$ , (69)

and is R_B,max = 6.1358 in the examples studied here, whereas the time derivative dR_B,max/dT has to be determined by simulation. The constant liquid volume V_S surrounds the bubbles of an elementary cell in the first model:

${\text{V}}_{\text{S}}=\frac{16}{3}\text{π}\left(r_S^3-r_B^3\right)=\frac{16}{3}\text{π}{r}_{B,0}^3\left(\frac{\text{π}}{3\sqrt{2}}\frac{1}{{\text{c}}_{\text{V},0}}-1\right)$ . (70)

If only spherical foam is formed, the volume flow rate q with the maximal time derivative dR_B/dT is chosen. In this case, both quantities, the bubble radius R_B and its derivative dR_B/dT, have to be determined by simulation.

Figure 3 illustrates the simulation results, where polystyrene is assigned a Newtonian rheology. At time T = 0, the gas volume fraction is c_V,0 = 0.01, the bubble radius is r_B,0 = 100 μm and the initial concentration of nitrogen is ξ_1,0 = 0.072%. The initial concentration of nitrogen is chosen so that nitrogen is just in balance. The ambient pressure P_∞ decreases exponentially with time according

to Equation (65) and is constant after having reached the minimum. In Figure 3 , the bubble radius R_B (red), the lamella thickness H (green) and the foam volume V (blue) are plotted as functions of the time T over a period of 60 seconds. As expected, the bubble radius in the first model increases, the lamella thickness in the second model decreases and the foam volume increases continuously until a maximum value is reached. The model change from spherical foam to polyhedral foam is indicated by a point of discontinuity (black dotted line parallel to the y-axis) in the red/green graph. After 60 seconds, the foam volume has reached the maximum value.

In Figure 4, the concentration profiles of nitrogen and the aromatic compound are plotted at different times for the example described above. The concentration profiles are shown as functions of the liquid volume around a bubble. The liquid volume is normalized by the constant liquid volume that surrounds a bubble. At the beginning, nitrogen and the aromatic compound are uniformly dissolved in polystyrene. The partial pressure of nitrogen decreases with the ambient pressure continuously. Hence, the nitrogen concentration decreases quickly at the phase boundary between liquid and gas so that large concentration differences occur in the liquid. Due to the larger diffusion coefficient, nitrogen diffuses anyway faster into the bubbles than the aromatic compound. After 10 seconds, nitrogen is almost completely desorbed. In contrast to nitrogen, the partial pressure of the aromatic compound decreases more slowly. After 60 seconds, finally, the hydrodynamic and thermodynamic equilibrium is reached, i.e. pressure differences between the bubbles and their environment and concentration differences in the liquid are balanced (constant in space and time), see also Table 3.

In the following, aeration processes at the temperature 20˚C are calculated for intermediate viscosity and low viscosity systems that are representative for food systems: It is assumed that the liquid has a Newtonian rheology and the density and surface tension of water, i.e. ρ = 1000 kg/m³ and σ = 0.072 N/m, for the purpose of comparison. Nitrogen is dispersed into the liquid under pressure before foaming. First, mass and heat transport are neglected so that the foam growth is induced only by a pressure drop after dispersion of nitrogen. The necessary material and process data are taken from the Table 1, Table 2, unless otherwise specified in this subsection. In contrast to the literature [8] [9] [10] , Newtonian liquids of intermediate or low viscosity are considered. Numerical simulations are carried out for two-phase systems with the viscosity 500 Pas, 50 Pas, 5 Pas, 0.5 Pas, 0.05 Pas and 0.005 Pas, respectively, over a period of 30 seconds. Reynolds numbers Re ≈ 1 require that inertia terms on the left side of the Equations (33) and (39) have to be taken into account. In Figure 5, the

simulation results are plotted for the systems with the viscosity 500 Pas, 5 Pas and 0.05 Pas. It becomes apparent that the two-phase systems behave similarly during foaming and that they attain the same equilibrium state after foaming. The hydrodynamic equilibrium is reached after 10 or 20 seconds. Table 4 specifies the values of relevant material and process quantities after foaming and shows that the less viscous the liquid is, the faster the change from spherical to polyhedral foam occurs.

In addition, aeration processes at the temperature 20˚C are calculated for intermediate viscosity and low viscosity systems, where mass and heat transport from the liquid in the gas phase are included. It is assumed that the liquid has again a Newtonian rheology and the density and surface tension of water, ρ = 1000 kg/m³ and σ = 0.072 N/m. Nitrogen was dispersed under pressure before foaming and is in form of small spherical bubbles, but also in dissolved form in the liquid, whereas a flavoring substance with the properties of styrene is already dissolved in the liquid. Hence, the foam growth is induced by dispersion of nitrogen, vacuum application and diffusion of dissolved components. Since Newtonian liquids of intermediate or low viscosity are considered, inertia terms in the Equations (33) and (39) are taken along. Numerical simulations are carried out for two-phase systems with the viscosity 500 Pas, 50 Pas, 5 Pas, 0.5 Pas, 0.05 Pas and 0.005 Pas, respectively, over a period of 30 seconds, according to the material and process conditions given in the Table 1, Table 2, unless otherwise specified in this subsection. In Figure 6, the simulation results are plotted for the systems with the viscosity 500 Pas, 5 Pas and 0.05 Pas. The same observations

described above can also be made here: The two-phase systems, which differ only in the viscosity of the liquid, behave similarly during foaming and attain the same equilibrium state after foaming. Due to the mass and heat transport, a higher foam volume is reached than in the two-phase systems without diffusion. The concentration profiles of nitrogen and the flavoring substance in Figure 7 confirm that the flavoring substance desorbs faster in a low-viscosity liquid. To clarify the numerical simulation results, Table 5 contains the values of relevant material and process quantities after foaming. The hydrodynamic and thermodynamic equilibrium between the liquid and gas phase is reached after 20 seconds, i.e. pressure differences between the bubbles and their environment and concentration differences in the liquid are balanced.

Table 3 shows that another equilibrium state for the polymer system is reached, which can be explained by the higher temperature 240˚C, compared to the lower temperature 20˚C of the systems characterized in Table 5 after foaming. A local Reynolds number Re ≈ 1 is obtained in each case. Generally, the definition of a global Reynolds number would give values Re ≳ 1 in the Table 4, Table 5.

At the end, the simulation results are compared with the experimental results presented by Haedelt, Beckett and Niranjan [2] [3] . The authors investigated bubble formation in liquid tempered chocolate, which is an example of an intermediate viscosity food system. After the tempering procedure, nitrogen (or another gas) was injected into the chocolate mass under controlled pressure and finely dispersed using a stator and rotor arrangement. When the mixture was returned to atmospheric pressure, the gases dissolved in the chocolate desorbed, and bubbles were formed [2] [3] . The operating conditions were chosen as follows: system pressure p_∞,max = 4.5 bar, p_∞,min = 1.0 bar, chocolate temperature ϑ₀ = 30˚C [3] (p. E139). The remaining data not determined are taken from the Table 1, Table 2, but now pressure drop coefficient K_exp = 0.15 s⁻¹, bubble radius

r_B,0 = 25 μm and chocolate viscosity η₀ = 10 Pas. Furthermore, ρ = 1000 kg/m³ and σ = 0.072 N/m.

The simulations results plotted in Figure 8 over a period of 30 seconds show that a fine spherical foam is obtained. The maximal bubble radius is r_B,max = 0.0882 mm, and the gas volume fraction is c_V,t = 0.3073 in agreement with the experimental results [3] (p. E141) for aerated chocolate produced by using nitrogen. The low gas volume fraction c_V,t is especially a result of the moderate pressures p_∞,max and p_∞,min chosen, see also Table 6.

The literature model [8] [9] [10] is also applicable if the liquid has a non-Newtonian rheology. The equations of motion (9) and (24) for the spherical and polyhedral foam have to be modified according to the constitutive equation used. In [9] , the rheological behavior of a viscoelastic polymer liquid is described using the generalized contravariant Maxwell-Oldroyd model based on a discrete relaxation time spectrum λ_R,m and η_m with η₀ = Ση_m [19] [20] [21] [22] . The simulation results are presented and discussed in [9] [23] [24] and show that the elastic material properties of a viscoelastic liquid do not influence significantly the foaming process under the given conditions. The same equilibrium state is reached, but with delay in time.

5. Conclusions

In the present article, aeration processes of intermediate viscosity and low viscosity systems are calculated, based on a literature model. The literature model describes growth of a polymer foam under vacuum with support of a gas dispersed into the polymer liquid before foaming. Numerical simulations are carried out for two-phase systems that are of intermediate or low viscosity and that differ only in the viscosity at a temperature of 20˚C. The necessary material and process data are taken from the literature. Reynolds numbers Re ≈ 1 require that inertia effects in these systems have to be taken into account. The simulations show that the two-phase systems behave similarly during foaming and that they attain the same equilibrium state after foaming. The less viscous the liquid is, the faster the hydrodynamic and thermodynamic equilibrium is reached, i.e. pressure differences between the bubbles and their environment and concentration differences in the liquid are balanced. The numerical method used for solving the model equations was analyzed with respect to stability, consistency and convergence in order to ensure accurate numerical simulation results, which are of interest e.g. for food foams in industrial processes.

Acknowledgements

This research project was supported by the German Ministry of Economics and Technology and the FEI (Forschungskreis der Ernährungsindustrie e.V., Bonn). Project AiF 17125 N. This project is part of the cluster “Proteinschäume in der Lebensmittelproduktion: Mechanismenaufklärung, Modellierung und Simulation” funded by the FEI and the AiF.

Cite this paper

Gladbach, K., Delgado, A. and Rauh, C. (2017) Numerical Simulation of Foaming Processes. World Journal of Mechanics, 7, 297-322. https://doi.org/10.4236/wjm.2017.711024

References

Appendix

The convergence analysis for the finite difference method in the sections 3 and 4.1 is based on the properties of a special class of matrices, the class of the M-matrices. They are briefly introduced here with reference to the literature [14] [15] [16] [17] [25] . Let $A\in {ℝ}^{\text{n}\times \text{n}}$ be a matrix with the elements $a_{ij}\in ℝ$ and $i,j\in \text{I}=\left\{\text{1},\cdots ,n\right\}$ . The set I is called index set. It is written $A\ge B$ if $a_{ij}\ge b_{ij}$ for all $i,j\in \text{I}$ .