Numerical Simulation of Foaming Processes

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.


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  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.

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 3 ,0 e ges,0 V,0 Figure 1.Spherical foam model (left) and growth of a single spherical bubble (right) [9].
g,0 g,0 V,0 ges,0 g,0 f 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: 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.
The constant liquid volume e f V 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].

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 cur- 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) ( ) Similarly, the simplified momentum equation in radial direction r for the single velocity component u is given by ( ) 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 ( ) 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 ( ) ( ) ( ) 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η 0 D, where S is the extra stress tensor, D is the deformation velocity tensor and η 0 is the zero-shear viscosity of the liquid.The normal stress difference in (9) reads then At time t = 0, the temperature ϑ 0 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 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 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]: ( ) 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 ( ) , 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 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.
( ) 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: ,begin2 ,end1 16π 3 15 7 5 The liquid volume e f V 4 is transformed into a thin layer normal to the bubble surface with thickness The second condition (ii) for the model change leads to the time derivatives

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) ( ) ( ) with the elongational rate The momentum equation in axial direction z simplifies to 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 ( ) ( ) in accordance to [9].For a Newtonian liquid with zero-shear viscosity η 0 , the normal stress difference in (24) reads 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 ( ) The temperature equation and first boundary condition are given by

Model Equations in the Dimensionless Form
The model Equations ( 9)-( 16) for the spherical foam are subjected to the Lagrange coordinate transformation ( ) and transformed into dimensionless equations [9].Similarly, the model Equations ( 24)-( 29) for the polyhedral foam are subjected to the Lagrange coordinate transformation 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: 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 ( ) ( ) 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 ( ) , where ( ) Concentration equations ( ) ,0 0, 0 : Mass balance equations Temperature equation ( ) The pressure inside a bubble is determined by Polyhedral foam: Equation of motion at the phase boundary ( ) Newtonian liquid ,0 : , 0 :

Numerical Method for Solving the Model Equations
The models developed in [9]  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 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 N S h Y = , N ∈  , the set of all spatial grid points is defined: with the constant discretization parameter h is always World Journal of Mechanics 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, ( ) For the time discretization, the grid selected, where the discretization parameter end1 K v T = , K ∈  , is the tem- poral step size.The system ( 46) is integrated by one-step-θ methods [16] 46) is referred to as ( ) F Y T , the integration within the limits k T T = and For [ ] 0,1 θ ∈ , 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: ( ) ( ) : Inserting the expressions for ( ) ( ) gives the discretized equations for the variables ( ) For each discrete time point T k , the grid function , : The equations for the grid function are obtained when the error term in (51) is omitted: 1 2 1 1 Hence, the discrete solution is determined by a linear system of equations that has to be solved in each time step.The initial and boundary conditions are taken into account by the vector , k i h b , the first and the last component of which are only non-zero.The matrix L 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 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 2) The matrix L 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 Hence, theorem 3 is applicable, and the following estimation is valid: Estimation ( 55) is relevant for the convergence considerations in the next section.The Equation (46), 1 N 1 n ≤ ≤ − 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 −2 ).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.

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 ( ) , , , R h denotes the restriction operator that gives the exact solution of (45) in the spatial grid points at time T k , whereas , k i h ξ is the approximated solution.From (51)-( 52), it follows that the finite difference method is consistent.The maximum norm of the local truncation error , Now, the global discretization error is studied.Note that, from (51)-( 52), it follows , , By recursive application, the following in equation is obtained: The matrix norm of Using ( 55) and ( 62) in (61), it follows ( ) ( ) ( ) max due to .
The finite difference method is consistent and stable and thus convergent.The stability results from the fact that L are M-matrices and ( ) 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).

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 where K exp = 0.6 s −1 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.Table 2. Process data for the system polystyrene/styrene/nitrogen at 240˚C [9].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.
( ) ( ) and on the basis of a single bubble and the constant liquid volume around it in the second model (see Equation ( 5)), i.e.
The local considerations suggest the definition of a local Reynolds number that characterizes the foaming process.The Reynolds number is defined by , where q is the volume flow rate, ν η ρ = 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 The volume flow rate q is determined by the time derivative of the foam volume e ges V in Equation (66).The maximal bubble radius can be calculated ana- lytically, 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: ( ) 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  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.      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       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  ii a > for all I i ∈ and 0 ij a ≤ for all i j ≠ .
2) A is not singular and A −1 ≥ 0. , E i j ∈ , then i is called directly connected with j.Moreover, i is called connected with j if there is a chain Theorem 2 [17].If the matrix A satisfies the sign condition in definition 1 and if A is diagonally dominant or irreducibly diagonally dominant, then, A is an M-matrix.
Theorem 3 [17].Let A be an M-matrix and let w be a vector with Aw ≥ 1. Then,

K.
Gladbach et al.DOI: 10.4236/wjm.2017.711024301 World Journal of Mechanics mass transport of the volatile components from the liquid into a bubble is based
dimensionless quantities used for the polyhedral foam model are the thickness of a liquid lamella

3
[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 concrete difference scheme.If the right side of Equation (

L
satisfies the sign condition in definition 1.

1
World Journal of Mechanics

Figure 3 .
Figure 3. Bubble radius R B (red) in the first model, lamella thickness H (green) in the second model and foam volume V (blue) as functions of the time T for the two-phase system polystyrene/styrene/nitrogen.

Figure 4 .
Figure 4. Concentration profiles for nitrogen (left) and styrene (right) as functions of the normalized liquid volume around a bubble for the system polystyrene/styrene/nitrogen.

Figure 5 .
Figure 5. Bubble radius R B (red) in the first model, lamella thickness H (green) in the second model and foam volume V (blue) as functions of the time T for different two-phase systems without mass and heat transport during foaming.

Figure 6 .
Figure 6.Bubble radius R B (red) in the first model, lamella thickness H (green) in the second model and foam volume V (blue) as functions of the time T for different two-phase systems with mass and heat transport during foaming.
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 ϑ 0 = 30˚C[3] (p.E139).The remaining data not determined are taken from the

Figure 7 .
Figure 7. Concentration profiles for nitrogen (left) and a flavoring substance (right) as functions of the normalized liquid volume around a bubble for Newtonian liquids with the viscosity 500 Pas (blue) and 0.05 Pas (black), respectively.
using the generalized contravariant Maxwell-Oldroyd model based on a discrete relaxation time spectrum λ R,m and η m with η 0 = Ση 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.World Journal of Mechanics

Figure 8 .
Figure 8. Bubble radius R B (red) and foam volume V (blue) as functions of the time T for liqiuid tempered chocolate with dissolved gases.
called graph of the matrix A. If ( )

Definition 3 .z
A matrix A is called irreducible if each i I ∈ is connected with each j I ∈ .Definition 4. A matrix A is called diagonally dominant if ∈  with radius r ∈  and let ( ) { }

Table 3 .
Two-phase system polystyrene/styrene/nitrogen in equilibrium after foaming.

Table 4 .
1, Table2, unless otherwise specified in this subsection.In Figure6, the simulation results are plotted for the systems with the viscosity 500 Pas, 5 Pas and 0.05 Pas.The same observations Two-phase systems in equilibrium after foaming (without mass and heat transport).

Table 5 .
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, At the end, the simulation results are compared with the experimental results [3]sented 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

Table 5 .
Two-phase systems in equilibrium after foaming (with mass and heat transport).