# Oil-in-water emulsion stabilized by hydrolysed black soldier fly larvae proteins: Reproduction of experimental data via phase-field modelling

Lucas Sales Queiroz, Angelique Berthome, Hector Gomez, Betul Yesiltas, Antonio Fernandes De Carvalho, Federico Casanova, Brais Martínez-López

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Lucas Sales Queiroz, Angelique Berthome, Hector Gomez, Betul Yesiltas, Antonio Fernandes De Carvalho, et al.. Oil-in-water emulsion stabilized by hydrolysed black soldier fly larvae proteins: Reproduction of experimental data via phase-field modelling. Journal of Food Engineering, 2025, 387, pp.112331. 10.1016/j.jfoodeng.2024.112331. hal-04749222

HAL Id: hal-04749222

https://hal.inrae.fr/hal-04749222v1

Submitted on 22 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Oil-in-water emulsion stabilized by hydrolysed black soldier fly larvae proteins: Reproduction of experimental data via phase-field modelling

Lucas Sales Queiroz<sup>a,b</sup>, Angelique Berthome<sup>a</sup>, Hector Gomez<sup>c</sup>, Betul Yesiltas<sup>d</sup>, Antonio Fernandes de Carvalho<sup>b</sup>, Federico Casanova<sup>a,**</sup>, Brais Martinez-Lopez<sup>e,*</sup>

## ARTICLE INFO

**Keywords:**

Phase field model

Emulsion

Stabilizer

Phase separation

Cahn-Hilliard equation

Black soldier fly larvae

## ABSTRACT

A phase field model based on the Cahn-Hilliard equation was validated experimentally by comparison of the numerical data with experimental data of emulsions of oil in water stabilized with Black Soldier Fly Larvae Proteins (BSFLP) and protein hydrolysates. Among the model parameters, the surface tension was determined by the pendant drop method, and the mobility by two different experiments, one based on the Turbiscan Stability Index, and another one based on the history of the average $d_{3,2}$, both of them throughout a 48 h period. The initial condition was built from an experimental droplet size distribution measured prior to phase separation. Results show that the model is able to quantitatively reproduce the phase separation kinetics, and provide an intuitive graphical representation of the droplet growth. Application of this procedure to other systems will allow the generalization of phase field modelling as a predictive tool for food applications.

## 1. Introduction

Modelling of the phase separation of an emulsion requires a proper mathematical description of the droplet growth that occurs during the Ostwald ripening stage. This represents a significant challenge with respect to classical processes governed by mass and heat transfer (e.g. drying, baking, freezing) since, if a mechanistic approach is to be used, the problem turns de facto into multiphase fluid flow. In order to either circumvent, or reduce this complexity, several approaches have been developed (Vetter et al., 2013): i) population balance, which can be simplified into a lumped (0-d or position-independent) equation ii) kinetic equations similar to those used to model chemical reactions (1st, sedond order), and iii) mass transfer between the droplets and the dispersed phase (Yarranton and Masliyah, 1997). To the best of the author's knowledge, none of these approaches have undergone extensive experimental validation, especially for food applications. Another possibility are the phase-field models, such as the Cahn-Hilliard equation (Cahn and Hilliard, 1958). This equation is based on the principle that the droplet growth is the result of spontaneous movements

triggered by a chemical potential gradient. One of the main advantages of this approach is its mechanistic nature, which means that the parameters required to describe the phase separation (the interface mobility, the surface tension and the interphase thickness) have a clear physical meaning. Additionally, its standard 2D formulation allows a very intuitive graphical representation of the phase separation phenomenon. Unfortunately, the Cahn-Hilliard equation includes a fourth order spatial derivative, which is not standard in either heat/mass transfer of fluid mechanics. Hence, despite the increasing popularity of mathematical modelling in the last decades, its implementation can still be non-trivial, and its numerical resolution may require a significant amount of computing power. To the best of the author's knowledge, and despite the significant amount of theoretical work about this framework, and its potential applications (Chella and Viñals, 1996; Fan et al., 1998; Kim et al., 2016; Lee et al., 2014), this approach has not undergone experimental validation, at least for food related applications. Therefore, with the exception of the surface tension, it is not possible to find values of the other model parameters that would allow to use this approach e.g. to estimate the phase separation time.

<footnote><sup>a</sup> Research Group for Food Production Engineering, National Food Institute, Technical University of Denmark, Søltofts Plads, 2800, Kongens Lyngby, Denmark</footnote>

<footnote><sup>b</sup> Departamento de Tecnologia de Alimentos, Universidade Federal de Viçosa (UFV), 36570-900, Viçosa, Minas Gerais, Brazil</footnote>

<footnote><sup>c</sup> School of Mechanical Engineering, Purdue University, 585 Purdue Mall, West Lafayette, IN, 47907, USA</footnote>

<footnote><sup>d</sup> Research group for Bioactives – Analysis and Application, National Food Institute, Technical University of Denmark, Kemitorvet, 2800, Kongens Lyngby, Denmark</footnote>

<footnote><sup>e</sup> IATE Université de Montpellier, INRAE, Institut Agro, 34060 Montpellier, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>** Corresponding author.</footnote>

<footnote>E-mail addresses: feca@food.dtu.dk (F. Casanova), brais.martinez-lopez@umontpellier.fr (B. Martinez-Lopez).</footnote>
On the context of emulsions for food applications, proteins are known to behave as surfactants, forming a thick layer at the emulsion's interface, which decreases the interfacial tension (Berton et al., 2011; Queiroz et al., 2021a). Traditionally, food industry has used conventional protein sources, such as milk and eggs, or even synthetic compounds as surfactants. The growing interest for lower carbon footprint ingredients has originated a trend that is progressively replacing, among others, conventional surfactants by environment-friendly protein sources (Gravel and Doyen, 2020; McClements et al., 2017). To this respect, Black Soldier Fly Larvae (BSFL) protein concentrate has shown promising properties that make it a good surfactant: water and oil holding capacity (Bessa et al., 2020), foaming capacity and stability (Queiroz et al., 2021b) and emulsifying properties (Mshayisa et al., 2021). The latter may be further improved with emerging technologies as sonication and ohmic heating (Queiroz et al., 2021a). Whereas whole insects have low acceptability, the use of food processing methods such as protein hydrolysis, could be an interesting strategy to insert insect proteins into food (Gravel and Doyen, 2020; Nongonierma and Fitzgerald, 2017).

The goal of this work is to show, at which extent, a simple version of the phase field model is able to reproduce experimental phase separation data of oil water emulsions, stabilized by BSFLP protein. Successful completion of the aforementioned goal allows to show the following.

* Choice and sampling of variables that can be measured both experimentally and numerically, and that are representative of the state of the phase separation.

* Determination of the model parameters, as well as the initial condition by means of independent experiments.

* Comparison between the numerical results and experimental data from two different experiments.

* To present ways in which the mathematical formulation can be modified to better represent the physical phenomena that.

In this simple version, the equation that describes the variation of the phase-field over time is solved in a bi-dimensional squared domain, without being coupled to the Navier-Stokes equations. This means that the velocity field that would be generated by the buoyancy forces is neglected, and hence the only driving force is the gradient of chemical potential in the Cahn-Hilliard equation. In the latter, the free energy has been computed by means of a simple double well potential function. By showing how the entire parameter set can be determined experimentally, this work intends to encourage its application to other systems, so that in the future, phase field modelling can be used as a predictive tool for food applications.

## 2. Numerical description of the phase separation phenomenon

The phase separation kinetics of an oil/water mixture can be described mathematically by means of the Cahn-Hilliard equation (Cahn and Hilliard, 1958):

$$ \frac{\partial \varphi}{\partial t} = M \nabla^2 \mu \qquad (1) $$

where $\varphi$ is the phase field (-) that represents the oil/water mixture. Accordingly, $\varphi$ takes the value of 1 in the regions occupied by oil, -1 for the regions occupied by water. The regions where $-1 < \varphi < 1$, correspond to the interface between oil and water. The other symbols in equation (1) are the laplacian operator $\nabla^2 (\text{m}^{-2})$, the mobility $M (\text{m}^5 \text{J}^{-1} \text{s}^{-1})$ and the chemical potential $\mu (\text{J m}^{-3})$, which is given by the expression:

$$ \mu = \sigma \left( \frac{1}{\epsilon} F'(\varphi) - \epsilon \nabla^2 \varphi \right) \qquad (2) $$

where $\sigma$ is the surface tension (J m$^{-2}$), $\epsilon$ is the interface thickness (m)

and $F'$ is the derivative of a double well potential function, that in this case is given by the expression $F = 0.25(\varphi^2 - 1)^2$. In this work, the mobility, the surface tension and the interface thickness are considered independent of the phase field, which allows to write the equation in the following form (Bueno et al., 2018):

$$ \frac{\partial \varphi}{\partial t} = M \sigma \nabla^2 \left( \frac{1}{\epsilon} F(\varphi) - \epsilon \nabla^2 \varphi \right) \qquad (3) $$

Equation (3) is solved in a closed domain, hence the following boundary conditions apply in the whole boundary:

$$ M \nabla \mu \cdot \vec{n} = 0 \qquad (4) $$

$$ \nabla \varphi \cdot \vec{n} = 0 \qquad (5) $$

Equation (4) enforces mass conservation, and equation (5) imposes a 90° angle between the boundary and the oil/water interface. In both equations, $\nabla$ represents the gradient operator ($\text{m}^{-1}$) and $\vec{n}$, a unit vector perpendicular to the boundary.

For simplification purposes, equation (3) will be solved in its dimensionless form, given by:

$$ \frac{\partial \varphi}{\partial \tau} = \nabla^2 (F(\varphi) - \epsilon^2 \nabla^2 \varphi) \qquad (6) $$

Under this form, the model parameters are grouped into two dimensionless numbers, the dimensionless time $\tau$, and the dimensionless interphase thickness $\epsilon$, which are respectively given by:

$$ \tau = \frac{M \sigma t}{\epsilon L^2} \qquad (7) $$

$$ \epsilon = \frac{\epsilon}{L} \qquad (8) $$

The interested reader can refer to the supplementary material, where the procedure to recast equation (3) into its dimensionless form given by equation (6) is described.

### 2.1. Numerical solution

Equation (6) was solved on a two-dimensional square, by means of a semi-implicit time integration algorithm based on convex-splitting (Elliott and Stuart, 1993). The spatial discretization is based on a second-order finite-difference algorithm. The linear systems of equations are solved using a fast Poisson solver (Strang, 2007). This algorithm was implemented in Matlab (2021a). A grid of 200 x 200 points was used for all the simulations. Two-time scales were studied: i) a short time kinetics, where the phase separation was simulated within the $\tau$ interval between 0 and 0.01, with a step of $2.5 \cdot 10^{-5}$.

### 2.2. Initial conditions and simulation domain

The phase separation was simulated in a normalized squared domain equivalent to $2.83 \cdot 10^{-5}$ m. The initial size distribution of the oil droplets was built from the experimental measurements of the DH12.5 (see §3.5). This is equivalent to a diameter size of $1.6 \pm 0.0001 \cdot 10^{-7}$ m (based on the experimental values).

### 2.3. Graphical representation of the numerical results and numerical data extraction

The numerical results of every time step were represented graphically in contour plots, such as the ones shown in Fig. 1a. Then, the numerical droplet diameter distributions $d_{[3,2],num}$, and average droplet surface $A_{num}$ and the average were extracted from the contour plots created for each time step (Fig. 1b). This was done by measuring the size of every oil droplet present in the image, by means of an image analysis
<figure><img src="image_4.png" /><figcaption>Fig. 1. (a) contour plot showing the oil droplets (white) and the water (black) phases after resolution of equation (6) for τ = 0.058. (b) Measurement of the droplet characteristic dimension with matlab function imcircles. (c) resulting distribution of the characteristic dimension.</figcaption></figure>

algorithm based on the Matlab® function imfindcircles. Finally, three numerical variables were deemed of interest to characterize the degree of advancement of the phase separation: the numerical dimensionless turbidity index, TSI<sup>*</sup><sub>num</sub>, diameter distribution d<sup>*</sup><sub>[3,2],num</sub> and surface A<sup>*</sup><sub>num</sub>. Their definitions are:

$$TSI_{num}^{*} = \frac{d_{[3,2],num} - d_{[3,2],num,0}}{d_{[3,2],num,\infty} - d_{[3,2],num,0}} \qquad (9)$$

$$d_{[3,2],num}^{*} = \frac{d_{[3,2],num}}{d_{[3,2],num,0}} \qquad (10)$$

$$A_{num}^{*} = \frac{A_{num}}{A_{num,0}} \qquad (11)$$

## 2.4. Model parameters

The value of the dimensionless interphase thickness was set to 0.006, which, according to the size of the simulation domain, corresponds to an actual interphase thickness of 1.7·10<sup>-7</sup> m or 170 nm.

The surface tension σ was determined experimentally from the independent experiment described in ¶3.3. The mobility M was determined by minimization of the norm of the squared residuals (equation (12)) between the experimental and numerical i) dimensionless turbidity index TI<sup>*</sup> and ii) d<sup>*</sup><sub>[3,2]</sub>:

$$SSQR = \sum_{i=1}^{nExp} \left( d_{[3,2],exp}^{*} - d_{[3,2],num}^{*} \right)^2 \qquad (12)$$

$$SSQR = \sum_{i=1}^{nExp} \left( TSI_{[3,2],exp}^{*} - TSI_{[3,2],num}^{*} \right)^2$$

The goodness of fit was assessed by equations 1-5 calculating the Root Mean Squared Error (RMSE), between, according to:

$$RMSE = \sqrt{\frac{SSQR}{n_{Exp}}} \qquad (13)$$

Notice that, since the SSQR is calculated from the dimensionless variables, the resulting RMSE is also dimensionless and within the interval 0-1. Hence, a better fit would be represented by an RMSE value close to 0.

# 3. Experimental materials and methods

## 3.1. Preparation of BSFL protein and hydrolysates

Black Soldier Fly larvae partially defatted flour was provided by nextProtein, (Paris, France). BSFL larvae were fed on biowaste before being separated from their feeding substrate by sieving when harvesting stage was reached. Then, a partially defatted BSFL flour was recovered after larvae processing: blanch, drying and pressing (Queiroz et al., 2021a). Alcalase® was provided by Novozymes (Krogshoejvej 36, Bagsvaerd, Denmark). For all dispersions, distilled water was used. BSFL defatted flour was obtained by mixing twice the BSFL flour for 2 h at 40 °C with ethanol (ratio 1:20, w/v), recovered with a filter paper and evaporated under fume hood overnight. Then, BSFL protein (BSFLP) was extracted according to (Queiroz et al., 2021a) method with slight modifications. Briefly, BSFL defatted flour was stirred in 0.25 M NaOH (ratio 1/20, w/v) during 1 h at 40 °C, 300 rpm. The resulting mixture was centrifuged at 2500×g, 4 °C for 20 min. This alkaline extraction was repeated on the recovered pellet. The supernatant from both alkali extractions were combined and proteins were isoelectric precipitated by adding 2 M HCl. The precipitate was centrifuged at 1300×g, 4 °C for 15 min, washed twice with distilled water and freeze-dried. A 1 g mL<sup>-1</sup> BSFLP (50 mL) dispersion was homogenized by stirring at 40 °C for 1 h at 450 rpm and overnight at room temperature. Temperature, and pH were adjusted by heating at 50 °C during 1 h and adding 0.5 M NaOH
until pH8 respectively. The dispersion was hydrolysed at 50 °C and pH8 with Alcalase (E/S ratio 3:100) at two DH: 4% and 12.5% by adding 0.337 mL and 1.165 mL of 0.5 M NaOH respectively, according to pH stat method. The DH was defined in equation (14):

$$DH = B \times x \times \frac{1}{\alpha} \times x \times \frac{1}{MP} \times x \times \frac{1}{htot} \times 100\% \quad (14)$$

Where B (mL) was the amount of base added during hydrolysis, Nb (N) the normality of the base, MP (g) was the mass of protein, htot (mEq/g) the total peptide bonds equals to 10.53 mEq/g for BSFLP according to amino acid composition found by (Queiroz et al., 2021a). 1/ α was the average degree of dissociation of the α-amino groups, equals to 1.13 at pH8 and 50 °C (the extrinsic conditions used in this study) for protein according to (J, 1986) and used in other studies adapted to their own conditions (Spellman et al., 2003). Hydrolysis was stopped by heating the samples at 85 °C for 20 min. The hydrolysates were cooled in ice prior to be centrifuged at 10,000 g, 10 °C during 15 min. Finally, supernatants were freeze-dried, DH4 and DH12.5 refer to the resulted powders at respective degree of hydrolysis and BSFLP hydrolysates (BSFLPH) to both hydrolysates. Prior to all experiment, pH of solutions made for BSFLP was adjusted to pH8 using 0.5M NaOH and the solution was vacuum filtrated using qualitative filter paper (Advantec). All BSFLP and BSFLP hydrolysed (BSFLPH) dispersions were homogenized by stirring (450 rpm) at 40 °C for 1 h and then overnight at room temperature.

## 3.2. Emulsion preparation

BSFLP and BSFLPH, 0.2 g mL⁻¹ dispersions, were emulsified with 0.05 m³ m⁻³ of sunflower oil. First emulsification was applied using Ultra-Turrax (IKA) dispersing device at 16,000 rpm for 3 min during which oil was slowly poured within the first minute resulting in a pre-emulsion then a high-pressure (12 kpsi) Microfluidizer® homogenizer (Microfluidics) was applied where samples were passed 3 times according to (Queiroz et al., 2021b).

## 3.3. Surface tension and dilatational properties

Surface tension (σ) required in equation (3) was measured by the pendant drop method using an optical contact angle and contour analysis OCA 20 (Data physics instruments GmbH, Filderstadt, Germany) instrument. 0.1 g mL⁻¹) BSFLP and BSFLPH dispersions were dispensed at the lower end of a dosing syringe (diameter 1.83 mm) in a cubic glass cell filled with sunflower oil (density: 0.918 g/cm³). σ was calculated according to Young-Laplace equation and plotted as a function of time:

$$\Delta P = \sigma \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \quad (15)$$

Where ΔP is the Laplace pressure (Pa), σ is the surface tension (J m⁻²) and R₁ and R₂ are the principal curvature radii of the pendant drop (m).

## 3.4. Multiple lights scattering measurements

Emulsion physical stability was measured using Turbiscan Lab Expert (Formulaction) and determined with Turbiscan Stability Index (TSI) according to (Queiroz et al., 2021b). A volume of 20 mL of BSFLP and BSFLPH emulsions were poured in tubes, closed with a lid, and gently poured in the device. TSI was measured for 2 days, and the physical stability of the emulsion was checked by the scanning program as follows: 3 scan every 2 min and 30 s, 61 scans every 2 min, 25 scans every 30 min, 25 scans every 1 h and 10 scans every 1 h. The dimensionless experimental turbidity index over the 48-h period was then calculated as below:

$$TSI_{exp}^* = \frac{TSI_{exp} - TSI_{exp,0}}{TSI_{exp,48h} - TSI_{exp,0}} \quad (16)$$

## 3.5. Droplet size

Droplet size was determined at day 0, day 1 and 2 after emulsification using Mastersizer 2000 (Malvern Panalytical) laser particle size analyzer. Several droplets of the emulsion were poured in a dispersion unit filled with distilled water and stirred at 2000 rpm. The experimental surface area moment diameter d[3,2],exp (Eq. (17)) was calculated based on droplet diameter d (nm):

$$d_{[3,2],exp} = \frac{\sum d^3}{\sum d^2} \quad (17)$$

## 3.6. Statistical analysis

For all statistical analyses ANOVA was applied, followed by pair comparisons using Tukey's test with a significance level p < 0.05 were applied using Statistical Package for the Social Sciences software (SPSS 22.0, SPSS Inc).

# 4. Results and discussions

## 4.1. Parameter determination and model validation

Simulation of the phase separation phenomenon by resolution of the system of equations (1)–(8) requires the parameter set formed by σ, M and ε and an initial condition. One of the goals of this work is to show that determination of these quantities experimentally allow to obtain adequate numerical reproductions of the phase separation. The whole procedure has been represented by means of a diagram in Fig. 2: the surface tension σ has been determined from an independent experiment (¶3.3 and ¶4.1.2). The interface thickness ε has been selected following the practical and physical criteria given in ¶2.4 and discussed in ¶4.1.5. Finally, the initial condition is a numerical representation of the droplet size distribution at time 0, and has also been determined experimentally (¶2.2). The mobility M has been determined by comparison of the experimental and numerical emulsion stability with two different techniques TSI* (¶4.1.3) and d[3,2]* (¶4.1.4).

### 4.1.1. Initial droplet size distribution

The initial droplet size distribution has a significant influence on the phase separation speed. More specifically, the presence of droplets larger than the average, even in a relatively small proportion, can make the separation process much faster. Because of this, for a reliable model validation, it is necessary to construct a numerical droplet size distribution that can be fed to the model as initial condition that is representative of the emulsion system. As it has been described in ¶2.2, the initial condition used for the simulations was built from the experimental droplet size distribution of the emulsion stabilized with DH12.5. Since the standard deviation is small, the distribution was homogeneous around the average. This average value, together with all the experimental d[3,2] determined from the experiments described in ¶3.5 are displayed in Fig. 3.

### 4.1.2. Experimental determination of surface tension

As shown by equation (3), the mathematical description of phase separation considers the surface tension as constant, which decreases significantly the numerical complexity of the problem. To this respect, the surface tension histories for the sunflower oil-water interface stabilized with BSFLP and BSFLPH, as determined from the experiments described in ¶3.3 are shown in Fig. 4.

As previously noticed by (Wang et al., 2021), BSFLP drops below 0.011 J m⁻² within the first 10 min and eventually reaches a
<figure><img src="image_5.png" /><figcaption>Fig. 2. diagram showing the model validation procedure with numerical (black) and experimental (purple) inputs and outputs.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. BSFLP and BSFLPH emulsion droplet sizes at day 0,1 and 2 measured with a mastersizer laser particle analyser as described in ¶3.5. Different letter indicates a significant difference (p < 0.05) between samples for the same day.</figcaption></figure>

quasi-equilibrium state around 0.087 J m⁻². This is probably caused by its surface hydrophobicity and less charged protein, which allows a denser surface coverage (Gould and Wolf, 2018). Use of BSFLH resulted in a σ drop within the first 10 min, similar to BSFLP, but the

<figure><img src="image_7.png" /><figcaption>Fig. 4. Interfacial tension of BSFLP and BSFLPH at sunflower O/W interface as a function of time measured by the pendant drop method described in ¶3.3.</figcaption></figure>
quasi-equilibrium state took longer, and its value was higher than with BSFLPH. Between DH4 and DH12.5, the former shows a sharper drop during the first 10 min, but they both stabilize around comparable values after 90 min This fact could be due to the higher surface hydrophobicity and possibility favourable amino acid disposition to better interact with both phases. For the sake of comparison, σ between water and sunflower oil, in the absence of emulsifiers was ~0.024 J m⁻².

In all the three tested cases, the variation of the surface tension was observed over a period of ca 100 min, with respect to a separation time of at least 48 h, which means that the assumption of a constant surface tension is reasonable. Table 1 shows the constant values that were selected, which correspond to the observed quasi-equilibrium states for each of the proteins reached after 48 h shown in Fig. 4

### 4.1.3. Determination of interface mobility from TSI

The Turbiscan Stability Index is displayed as function of time in Fig. 5 in two forms: the experimental data directly measured from the turbiscan (Fig. 5a), as well as the normalized TSI*, as calculated from equation (16) (Fig. 5b,c and d). As it can be seen (Fig. 5a), the BSFPLH emulsions are more stable than the BSFLP, since their TSI curves remain, during the whole 48 h period of observation, below the TSI BSFLP curve. More specifically, DH4 shows the greatest stabilizing capability, with a final TSI value of 5.1, followed by DH12.5 and BSFLP with 6.2 and 7.3 respectively. Moreover, emulsions stabilized by BSFLPH keep a TSI value under 3 for a longer time. More specifically, DH4 and DH12.5 only reach a value of 3 after ~19 h and ~23 h respectively, whereas BSFLP reaches it after ~10 h. The shape of the curves is in agreement with the observed droplet size increase. Estimation of the mobility *M* from the normalized TSI* curves (equations (12) and (16)) confirms these observations partially, as the emulsions stabilized by BSFPLH show lower mobility values than by BSFLP, yet not very significantly. As it is shown in Table 1, all the values are on the vicinity of 5·10⁻²² m⁵ J⁻¹ s⁻¹. The three curves are comparable with respect to their ability to reproduce the experimental data, with RMSE of 0.058, 0.057 and 0.062 for the emulsions stabilized with BSFLP, DH4 and DH12.5 respectively. A less quantitative analysis reveals that, among the three, the model reproduces the shape of the DH4 curve best, followed by the BSFLP and DH12.5. The latter shows an inflection point around 16 h that is not really reproduced by the model. This could be explained by either numerical or experimental reasons. Numerical reasons include i) a lack of precision on the numerical droplet sampling (¶2.3) and ii) absence of velocity fields, which cannot be established because the momentum and the continuity equations have not been included in the model formulation. With respect to reasons related to experimentation, the data used to represent stability is an average of the whole sample. The use of averaged or local measurements, may result in different parameter values, or even to deviation of the model predictions from experimental data (Anagnostara et al., 2022; Martinez-Lopez and Mauricio-Iglesias, 2022).

In order to obtain more realistic mobility values that reflect the different stabilization indexes discussed at the beginning of this section, the TSI should have had to be normalized with respect to a value that can

<caption>**Table 1**</caption>

Dimensionless interphase thickness used for the simulations, interphase thickness, surface tension determined experimentally, mobility determined by comparison between the numerical results and the experimental TSI* (a) and *d*[3,2] (b), as well as RMSE between the numerical and the experimental data.

<table><thead><tr><th>Sample</th><th>ε (−)</th><th>ε (m)</th><th>σ (J m<sup>-2</sup>)</th><th>M (m<sup>5</sup> J<sup>-1</sup> s<sup>-1</sup>)</th><th>RMSE (−)</th></tr></thead><tbody><tr><td>BSFLP</td><td></td><td></td><td>0.0087</td><td><sup>a</sup>7.0·10<sup>-22</sup></td><td><sup>b</sup>5.5·10<sup>-23</sup></td><td><sup>a</sup>0.058</td></tr><tr><td>DH4</td><td>0.006</td><td>1.7·10<sup>-7</sup></td><td>0.0095</td><td><sup>a</sup>7.0·10<sup>-22</sup></td><td><sup>b</sup>7.6·10<sup>-23</sup></td><td><sup>a</sup>0.057</td></tr><tr><td>DH12.5</td><td></td><td></td><td>0.0096</td><td><sup>a</sup>4.5·10<sup>-22</sup></td><td><sup>b</sup>3.5·10<sup>-22</sup></td><td><sup>a</sup>0.062</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Turbiscan stability index (TSI) as a function of time (a) for BSFLP (blue), DH04 (yellow) and DH12.5 (orange). Experimental (●) and predicted (continuous black line) normalized Turbiscan Stability Index (TSI*) as a function of time for BSFLP (b), DH4 (c) and DH12.5 (d) emulsions. Experimental determination by light scattering measurements described in ¶3.4.</figcaption></figure>
be deemed representative of a complete separation. This is rather inconvenient, and one of the major points of improvement for future work, since one of the possible applications of this model is the prediction of phase separation time.

### 4.1.4. Determination of interface mobility from droplet size history

This stabilization capability suggested by the TSI values after 48 h is reflected on the values of $d[3,2]$ reported in Fig. 3. Initially, smaller droplets were formed with BSFLPH, with a similar average size around 0.16 nm compared to BSFLP (0.27 nm). The $d[3,2]$ values reported in Fig. 3 were used to calculate the experimental values of $d[3,2]^*$ (equation (10)). These values are shown in Fig. 6 and give a more quantitative idea of the droplet size increase. The droplet sizes of emulsions stabilized with BSFLP and DH4 do not increase significantly ($d[3,2]^*$ values of 1.3 and 1.4 respectively after 48 h), whereas DH12.5 almost triplicates its size ($d[3,2]^* = 2.89$ after 48 h). Differently from the mobility values determined from the TSI*, the values determined from the comparison between experimental and numerical $d[3,2]^*$ (Table 1) do reflect the faster kinetics of DH12.5, with a mobility value of $3.5 \cdot 10^{-23} \text{ m}^5 \text{ J}^{-1} \text{ s}^{-1}$, almost 10 times greater than BSFLP and DH4, which show mobility values of $5.5 \cdot 10^{-23}$ and $7.6 \cdot 10^{-23} \text{ m}^5 \text{ J}^{-1} \text{ s}^{-1}$ respectively. An evaluation of the predictions solely based on the RMSE indicates that the predictions are much better for the BSFLP and the DH04 than for DH12.5, since their RMSE values are much lower. Nevertheless, and taking into consideration that the variation of $d^*[3,2]$ determined experimentally for the DH12.5 within the 0–48 h interval is much more significant than for BSFLP and DH04 (Fig. 3), in the author's opinion, the mobility coefficient determined for DH12.5 is more reliable despite the larger RMSE.

### 4.1.5. Interface thickness

As previously discussed, the value of $\epsilon$ was selected on the basis of physical and practical reasons. As this is the only parameter required to simulate the system in its dimensionless form, a lower value of $\epsilon$ results in a slower phase separation. It can be argued that i) an interphase thickness equivalent to roughly 1/3 of the average droplet diameter is too high, and ii) that ideally, the system should be simulated with $\epsilon$ calculated from a value of the interphase thickness $\epsilon$ obtained from experimental measurements. With respect to i), simulation of the phase separation with slightly larger values of $\epsilon$, more specifically within the interval 0.008–0.012 (not shown), results in the model depicting a stable mix between both phases after just a few hours, which is not physically realistic for oil and water. This behaviour at larger $\epsilon$ is partly due to the simplified double-well potential function $F = 0.25(\phi^2 - 1)^2$ whose derivative is required in equation (6). Because of this, and also considering

that the interphase thickness should be smaller than the radius, a value of $\epsilon = 0.006$ was selected. Should this result in a higher $\epsilon$ and hence a faster system, a more realistic double-well potential function would be necessary. With respect to ii) as of today, to the best of the author's knowledge, there are no literature values, either of experimental or numerical nature to compare with.

### 4.2. Numerical reproduction of the phase separation kinetics

Once the parameter set has been determined, the numerical results can be represented as it is shown in Fig. 7. This is the classical representation of the solution of the equation system 1–8 in a bi-dimensional squared domain and is usually used in the literature dealing with the Cahn-Hilliard equation. The dimensionless time $\tau$ has been converted into time units by means of the parameter sets given in Table 1. This way, the times expressed in hours in Fig. 7 would be equivalent to the $\tau$ values of 0.0014 (a), 0.0027 (b), 0.0041 (c) and 0.0055 (d). The numerical results in Fig. 7a–d depict the increases on the $d[3,2]$ shown in Fig. 6.

Although it has been shown that the model can reproduce the experimental data in an acceptable manner (Figs. 5 and 6), results are restricted to the 0–48 h time period. Still, with the parameter set determined in this work, it is indeed possible to simulate the process until complete separation. Having said this, it is the authors opinion that such simulations would largely overestimate the time required for complete separation. Mostly, this is due to the assumptions used in the model formulation and are part of this work's limitations. They are discussed in the following section.

### 4.3. Applicability and limitations of this work

The following subsections explain ways in which the model could be made more realistic and generalizable. More specifically 4.3.1 discusses what would be required for the application of the present method/model to other systems. 4.3.2 deals with the lack of forces other than the potential of chemical gradient in the model formulation; 4.3.3 discusses the influence of a more realistic geometry, notably the extension to 3D, and 4.3.4 the use of an interface mobility that is a function of the local concentration of stabilizer. Finally 4.3.5 elaborates on how the free energy model used in the chemical potential could be made more thermodynamically consistent.

#### 4.3.1. Application to other systems

Basically all the parameters required by classical heat- or mass transfer model formulations, (e.g thermal conductivity, mass

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** History of $d^*[3,2]$ (a) and $A^*$ (b). Experimental data (scatter), predictions, and RMSE bounds from Table 1 (–). BSFLP ■ and blue, DH04 ◆ and yellow, DH12.5 ● and orange.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Numerical results of equation (3) with the parameter set corresponding to DH12.5 shown in Table 1a. 12 h (a), 24 (b), 36 h (c), 48 h (d) of separation. The x and z axis have been normalized with respect to the size of the simulation domain. Animated results available via this DOI link.</figcaption></figure>

diffusivity), are system or composition specific. This means that, application of these models to specific systems should only be done when reliable parameter values, or models that can estimate them as a function of e.g. their composition can be found in literature. The model parameter set made up by the surface tension $\sigma$, the interface mobility $M$ and the interface thickness $\varepsilon$ is specific to the sunflower oil – water system stabilized by the surfactants used in this work since their chemical structure is numerically represented in i) the parameters values, and ii) the initial droplet size distribution. Since these parameters, specially the interface mobility cannot be found in literature for other systems to compare with, simulation with other systems with this parameter set will lead to unverifiable results. On the other hand, the methodology used in this work to determine the parameter set from experiments, is indeed, translatable to other systems, by following the procedure summarized in Fig. 2. In a nutshell, to increase the likelihood of a successful parameter estimation, the initial droplet size distribution must be known, and the mobility must be determined from emulsion stability data.

### 4.3.2. Forces driving the phase separation

In this work, the phase separation has been modelled by means of the Cahn-Hilliard equation on a squared domain. As it can be seen in equations (1)–(3), a direct consequence of this is that the only force driving the phase separation is the gradient of chemical potential $\mu$, which means that the advective transport of the phase field variable $\varphi$ has been neglected. In a classical heat- or mass-transfer formulation, this would be akin to a physical phenomenon being exclusively governed by thermal conduction or mass diffusion respectively. The reason behind this, is that the Navier-Stokes equations are not part of the model formulation. Inclusion of the Navier-Stokes equations would have allowed to calculate the velocity field generated by the spontaneous movement of the fluids due to buoyancy. The buoyancy term $\rho \vec{g}$, acts as a force in the conservation of linear momentum equation, and accounts for the velocity field being generated due to the density difference between oil and water. The resulting velocity field would most likely lead

to laminar or creeping flow regimes, but it would still accelerate the phase separation by pushing the fluid with the largest density (in this case, the water), to the bottom of the domain. Combination of the Cahn-Hilliard with the continuity and the Navier-Stokes equations is currently doable with the current state of the art of the applications of the phase field model e.g. (Bueno et al., 2018; Sandoval Murillo et al., 2019).

### 4.3.3. Domain geometry

The domain used in the numerical simulations is a vast simplification of the setup were the experiments have taken place in several ways: i) it is squared, which is not the case of e.g. the turbiscan sample holder, ii) it only represents part of a larger physical system and iii) it is a bi-dimensional representation of three-dimensional objects.

It can be argued that, although it would make the simulation more realistic, precise reproduction of the geometry (i) is not that important, since the experimental sampling itself (e.g. to measure the particle size distribution), disturbs the system in a way that cannot be reproduced by the model, regardless of its complexity. Accordingly, (ii) it is enough to represent the domain partially. In this work we have used zero-flux boundary conditions (equations (4) and (5)), although periodic boundary conditions would have been equally acceptable, since both are standard techniques to make computations on domains smaller than the physical system, and both enforce mass conservation. As it can be seen in Fig. 7, coalescence in the edges of the domain behaves according to equations (4) and (5). The boundary node receives contributions exclusively from the nodes that lay inside the domain (equation (4)), and the oil/water interface and the boundary are forced to form a 90° angle (equation (5)). With respect to (iii), a significant difference can be expected between a bi-dimensional and a three-dimensional domain. Had the simulation domain been three-dimensional, the droplets would have been represented as spheroids rather than cylinders. This would increase the laplace pressure, likely leading to a faster phase separation, and hence a larger value of $M$. As the Cahn-Hilliard equation includes a fourth derivative in the space, simulation of the system in 3D, even in its simplest form (extension of the squared domain into a cube), would
result in a significant increase in the computing power required to run
the simulation. As an example, in this work, a mesh conformed by 200 x
200 elements, and a τ step of 2.5·10⁻⁵ m, has been deemed necessary in
order to simulate the domain with the current interface thickness value
of 1.7·10⁻⁷ m, and obtain meaningful results. For the τ period between
0 and 0.01, this amounts to over 1.5·10⁷ equations. Extension to 3D
assuming that the same amount of elements in the additional axis would
be sufficient would result in over 3.1·10⁹ equations to solve. Finally, the
numerical sampling procedure shown in Fig. 1 would have to be adapted
to spheres rather than circles, which is non-trivial.

4.3.4. *Constant interface mobility*

As explained in §4.1.2, experimentation shows that the stabilizers used in this work lower the surface tension of the water-sunflower oil interface. It might be expected that the protein concentration difference between the water and the sunflower oil would trigger the diffusive transport of proteins through the aqueous phase, leading them to colonize the interface and to lower the local phase separation speed. This would be numerically equivalent to a non-constant mobility coefficient *M* in equations (1)–(4), dependent on the local concentration of surfactant/protein. In the current model, the phase field variable φ is the mathematical equivalent to the volumetric fraction of oil and water. Therefore, a protein dependent mobility would require i) definition of a new passive scalar variable that represents the local protein concentration, ii) an equation that describes its transport by diffusion, and iii) a function that links the mobility parameter with the local protein concentration. Literature contains phase field or VOF that couple the phase field together a scalar that is either confined to one phase, or allowed to cross the interface (Bothe and Fleckenstein, 2013; Jain and Mani, 2023). Some works use a mobility coefficient that is dependent on the phase field itself (Júnior et al., 2020) but, to the best of the author's knowledge, it has never been attempted to express it in function of the concentration of one solute within one phase. Numerical challenges aside (Jain and Mani, 2023), reliability of the numerical results with respect to experimentation would depend on a robust determination of the mobility-concentration function. As this cannot still be found in literature, the latter should be found via experimentation (e.g. determination of *M* for emulsions with different protein concentrations), or by numerical trial and error.

4.3.5. *Free energy model*

As it can be seen in equation (2), the phase separation kinetics
depend on F', which is defined as the derivative of F, with respect to the
phase field variable φ, F is a mathematical expression whose purpose is
to model the fluids mutual immiscibility (Jacqmin, 1999) and as such,
its only requirement is that it contains two minima, which correspond to
the fluids two stable phases. In this work, the function F = 0.25(φ² - 1)²
has been used, as it has two minima, at φ = -1 (pure water) and at φ = 1
(pure oil). Nevertheless, this function is also a representative of Helm-
holtz's free energy (Lee et al., 2014) which is a function of the internal
energy and the entropy of the system. Since those thermodynamic
functions are themselves a function of the composition of each phase,
which is relatively well known, their derivatives can be calculated using
thermodynamic models such as the Flory-Huggins (Júnior et al., 2020)
or a relationship between the excess properties and the activity co-
efficients (Tosun, 2021). Thanks to the current state of the art of solution
thermodynamics, the latter might be calculated via the UNIFAC model,
which in principle spares their experimental determination. The use of a
more complex free energy model could improve the results and the
reliability of the predicted mobility value.

5. Conclusions

This study has provided evidence that the phase field model can be
validated experimentally from emulsion stability data. The mobility *M*,
which is the main responsible for the phase separation speed, has been

determined by means of two different experiments, and shown to have
comparable values. Unfortunately, at the time of publication of this
manuscript, no other values can be found in literature to compare with.
Hence, until parameter sets for other systems are made available, the
model will not be generalizable to other mixtures. The procedure
described in this work to determine the model parameter set is most
likely applicable to many other systems with practical interest.

The mathematical model used in this work is relatively simple: i) the
domain is bi-dimensional ii) the only force driving phase separation is
the chemical potential gradient and iii) the free energy is computed by
means of a simple double well potential function rather than an actual
free energy model. Since this is the first time the Cahn-Hilliard model
undergoes comparison against emulsion stability data, these shortcom-
ings are acceptable. The improvements that deal with the force driving
the phase separation and with a more realistic domain/geometry are
doable under the current state of the art of and, in the author's opinion,
should be included in a future iteration of this work. Having said this,
these improvements will come at the expense of a drastic increase on the
required simulation time and computing power. Finally, but equally
important, in order to advance in the understanding of these systems,
any improved model should undergo comparison against stability data
for emulsion mixtures other than the one presented in this work.

CRediT authorship contribution statement

**Lucas Sales Queiroz:** Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Conceptualization. **Angelique Berthome:** Writing – original draft, Investigation, Formal analysis. **Hector Gomez:** Writing – review & editing, Software, Conceptualization. **Betul Yesiltas:** Methodology, Conceptualization. **Antonio Fernandes de Carvalho:** Writing – review & editing, Supervision, Funding acquisition. **Federico Casanova:** Writing – review & editing, Supervision, Project administration, Funding acquisition, Conceptualization. **Brais Martinez-Lopez:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Formal analysis, Conceptualization.

Declaration of competing interest

The authors certify that they have NO affiliations with or involve-
ment in any organization or entity with any financial interest (such as
honoraria; educational grants; participation in speakers' bureaus;
membership, employment, consultancies, stock ownership, or other
equity interest; and expert testimony or patent-licensing arrangements),
or non-financial interest (such as personal or professional relationships,
affiliations, knowledge or beliefs) in the subject matter or materials
discussed in this manuscript.

Data availability

Data will be made available on request.

Acknowledgments

The authors are grateful for the financial support of CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior), CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais). The authors acknowledge the funding ERASMUS Internship for Angelique Berthome.

**Appendix A. Supplementary data**

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jfoodeng.2024.112331.
## References

Anagnostara, I., Stubbe, P.R., Martinez-Lopez, B., 2022. Practical identifiability analysis for a multiphase diffusion-evaporation model: dough baking in a pilot scale jet impingement oven. Food Bioprod. Process. 135, 238-257. https://doi.org/10.1016/j.fbp.2022.07.005.

Berton, C., Ropers, M.-H., Viau, M., Genot, C., 2011. Contribution of the interfacial layer to the protection of emulsified lipids against oxidation. J. Agric. Food Chem. 59, 5052-5061. https://doi.org/10.1021/jf200086n.

Bessa, L.W., Pieterse, E., Marais, J., Hoffman, L.C., 2020. Why for feed and not for human consumption? The black soldier fly larvae. Compr. Rev. Food Sci. Food Saf. 19, 2747-2763. https://doi.org/10.1111/1541-4337.12609.

Bothe, D., Fleckenstein, S., 2013. A Volume-of-Fluid-based method for mass transfer processes at fluid particles. Chem. Eng. Sci. 101, 283-302. https://doi.org/10.1016/j.ces.2013.05.029.

Bueno, J., Casquero, H., Bazilevs, Y., Gomez, H., 2018. Three-dimensional dynamic simulation of elastocapillarity. Meccanica 53, 1221-1237. https://doi.org/10.1007/s11012-017-0667-4.

Cahn, J.W., Hilliard, J.E., 1958. Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258-267. https://doi.org/10.1063/1.1744102.

Chella, R., Viñals, J., 1996. Mixing of a two-phase fluid by cavity flow. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 53, 3832-3840. https://doi.org/10.1103/physreve.53.3832.

Elliott, C.M., Stuart, A.M., 1993. The global dynamics of discrete semilinear parabolic equations. SIAM J. Numer. Anal. 30, 1622-1663. https://doi.org/10.1137/0730084.

Fan, D., Chen, L.-Q., Chen, S.P., Voorhees, P.W., 1998. Phase field formulations for modeling the Ostwald ripening in two-phase systems. Comput. Mater. Sci. 9, 329-336. https://doi.org/10.1016/S0927-0256(97)00158-4.

Gould, J., Wolf, B., 2018. Interfacial and emulsifying properties of mealworm protein at the oil/water interface. Food Hydrocoll 77, 57-65. https://doi.org/10.1016/j.foodhyd.2017.09.018.

Gravel, A., Doyen, A., 2020. The use of edible insect proteins in food: challenges and issues related to their functional properties. Innov. Food Sci. Emerg. Technol. 59. https://doi.org/10.1016/j.ifset.2019.102272, 102272-102272.

Jacqmin, D., 1999. Calculation of two-phase Navier-Stokes flows using phase-field modeling. J. Comput. Phys. 155, 96-127. https://doi.org/10.1006/jcph.1999.6332.

Jain, S.S., Mani, A., 2023. A computational model for transport of immiscible scalars in two-phase flows. J. Comput. Phys. 476, 111843. https://doi.org/10.1016/j.jcp.2022.111843.

Júnior, E.J.P., Staudt, P.B., Tessaro, I.C., Cardozo, N.S.M., 2020. A new approach to phase-field model for the phase separation dynamics in polymer membrane formation by immersion precipitation method. Polymer 186, 122054. https://doi.org/10.1016/j.polymer.2019.122054.

Kim, J., Lee, S., Choi, Y., Lee, S.-M., Jeong, D., 2016. Basic principles and practical applications of the Cahn-Hilliard equation. Math. Probl Eng. 2016, 9532608. https://doi.org/10.1155/2016/9532608.

Lee, D., Huh, J.-Y., Jeong, D., Shin, J., Yun, A., Kim, J., 2014. Physical, mathematical, and numerical derivations of the Cahn-Hilliard equation. Comput. Mater. Sci. 81, 216-225. https://doi.org/10.1016/j.commatsci.2013.08.027.

Martinez-Lopez, B., Mauricio-Iglesias, M., 2022. General guidelines for a successful joint determination of the diffusion and the partition coefficients in polymeric food contact materials based on optimal experimental design. J. Appl. Polym. Sci. 139, 51691. https://doi.org/10.1002/app.51691.

McClements, D.J., Bai, L., Chung, C., 2017. Recent advances in the utilization of natural emulsifiers to form and stabilize emulsions. Annu. Rev. Food Sci. Technol. 8, 205-236. https://doi.org/10.1146/annurev-food-030216-030154.

Mshayisa, V.V., Van Wyk, J., Zozo, B., Rodríguez, S.D., 2021. Structural properties of native and conjugated black soldier fly (Hermetia illucens) larvae protein via Maillard reaction and classification by SIMCA. Heliyon 7, e07242. https://doi.org/10.1016/j.heliyon.2021.e07242, 07242.

Nongonierma, A.B., FitzGerald, R.J., 2017. Unlocking the biological potential of proteins from edible insects through enzymatic hydrolysis: a review. Innov. Food Sci. Emerg. Technol. 43, 239-252. https://doi.org/10.1016/j.ifset.2017.08.014.

Queiroz, L.S., Casanova, F., Feyissa, A.H., Jessen, F., Ajalloueian, F., Perrone, I.T., de Carvalho, A.F., Mohammadifar, M.A., Jacobsen, C., Yesiltas, B., 2021a. Physical and oxidative stability of low-fat fish oil-in-water emulsions stabilized with black soldier fly (hermetia illucens) larvae protein concentrate. Foods 10. https://doi.org/10.3390/foods10122977, 2977-2977.

Queiroz, L.S., Regnard, M., Jessen, F., Mohammadifar, M.A., Sloth, J.J., Petersen, H.O., Ajalloueian, F., Brouzes, C.M.C., Fraihi, W., Fallquist, H., de Carvalho, A.F., Casanova, F., 2021b. Physico-chemical and colloidal properties of protein extracted from black soldier fly (Hermetia illucens) larvae. Int. J. Biol. Macromol. 186, 714-723. https://doi.org/10.1016/j.ijbiomac.2021.07.081.

Sandoval Murillo, J.L., Osen, R., Hiermaier, S., Ganzenmüller, G., 2019. Towards understanding the mechanism of fibrous texture formation during high-moisture extrusion of meat substitutes. J. Food Eng. 242, 8-20. https://doi.org/10.1016/j.jfoodeng.2018.08.009.

Spellman, D., McEvoy, E., O'Cuinn, G., FitzGerald, R.J., 2003. Proteinase and exopeptidase hydrolysis of whey protein: comparison of the TNBS, OPA and pH stat methods for quantification of degree of hydrolysis. Int. Dairy J. 13, 447-453. https://doi.org/10.1016/S0958-6946(03)00053-0.

Strang, G., 2007. Computational science and engineering. In: Computational Science and Engineering. Wellesley-Cambridge Press, Wellesley, MA.

Tosun, I., 2021. Excess Mixture Properties and Activity Coefficients. Elsevier, The Netherlands.

Vetter, T., Iggland, M., Ochsenbein, D.R., Hänseler, F.S., Mazzotti, M., 2013. Modeling nucleation, growth, and Ostwald ripening in crystallization processes: a comparison between population balance and kinetic rate equation. Cryst. Growth Des. 13, 4890-4905. https://doi.org/10.1021/cg4010714.

Wang, J., Jousse, M., Jayakumar, J., Fernández-Arteaga, A., de Lamo-Castellví, S., Ferrando, M., Güell, C., 2021. Black soldier fly (hermetia illucens) protein concentrates as a sustainable source to stabilize O/W emulsions produced by a low-energy high-throughput emulsification technology. Foods 10. https://doi.org/10.3390/foods10051048, 1048-1048.

Yarranton, H.W., Masliyah, J.H., 1997. Numerical simulation of Ostwald ripening in emulsions. J. Colloid Interface Sci. 196, 157-169. https://doi.org/10.1006/jcis.1997.5186.