Analytical and numerical solutions of pore formation in elastic food materials during dehydration

In this paper, we describe a model for pore formation in food materials during drying. As a proxy for fruits and vegetables, we take a spherical hydrogel, with a stiff elastic skin, and a central cavity filled with air and water vapour. The model describes moisture transport coupled to large deformation mechanics. Both stress and chemical potential are derived from a free energy functional, following the framework developed by Suo and coworkers. We have compared Finite Volume and Finite Element implementations and analytical solutions with each other, and we show that they render similar solutions. The Finite Element solver has a larger range of numerical stability than the Finite Volume solver, and the analytical solution also has a limited range of validity. Since the Finite Element solver operates using the mathematically intricate weak form, we introduce the method in a tutorial manner for food scientists. Subsequently, we have explored the physics of the pore formation problem further with the Finite Element solver. We show that the presence of an elastic skin is a prerequisite for the growth of the central cavity. The elastic skin must have an elastic modulus of at least 10 times that of the hydrogel. An initial pore with 10% of the size of the gel can grow to 5 times its initial size. Such an increase in porosity has been reported in the literature on drying of vegetables, if a dense hard skin is formed, known as case hardening. We discuss that models as presented in this paper, where moisture transport is strongly coupled to large deformation mechanics, are required if one wants to describe pore/structure formation during drying and intensive heating (as baking and frying) of food materials from first principles.


Introduction
In the last decade, researchers are starting to understand how large mechanical deformation is interacting with (de)hydrating vegetables (Aregawi et al., 2014;Karunasena et al., 2014;Liu et al., 2015;Singh et al., 2015;Gulati and Datta, 2015) ( Van der Sman, 2015;Gulati et al., 2016;Mahbubur Rahman et al., 2018;Defraeye and Radu, 2018).For a quantitative understanding of this phenomenon, we think the multiphysics theory developed by Suo and coworkers is of seminal importance (Hong et al., 2008(Hong et al., , 2009;;Zhao et al., 2008;Marcombe et al., 2010).Their theory describes how to couple large deformation mechanics to moisture transport.Commonly, these multiphysics models are applied to the swelling of synthetic hydrogels, and the application of the theory to food materials is scarce (Van der Sman, 2015;Jin and van der Sman, 2022).The theory centers on a free energy function, which extends the Flory-Rehner theory but allows for inhomogeneous deformation.From the free energy functional the driving forces for momentum and mass transport are derived.As such the theory has been successful in describing phenomena as wrinkling (Weiss et al., 2013;Jin et al., 2015), which is especially extensive if the material has a stiff outer surface and soft inner part.
Another phenomenon due to the coupling of large deformations and dehydration is the development of pores during drying of polymeric food materials or vegetables (Maurice et al., 2017) (Jacob et al., 2016;Nguyen et al., 2018;Both et al., 2020;Siemons et al., 2020), (Zogzas et al., 1994;Karathanos et al., 1996;Khalloufi et al., 2015;Joardder et al., 2018).A requirement for this pore formation is the development of an elastic skin (Siemons et al., 2020), which can happen either that the skin forms a gel or that it enters the glassy state (Siemons et al., 2022).Further shrinkage of the drying food material, leads to buildup of elastic stresses in the skin, and underpressures in the core.Pre-existing pores, such as intercellular junctions in vegetables, will expand if the underpressure exceeds a critical value.
During vegetable drying, the formation of a rigid skin is called case hardening (Gulati and Datta, 2015).Some vegetables and fruits have a skin by nature, like cherries and plums.Especially, plums are known to form wrinkles during drying, which is taken as a prototypical example of wrinkling in soft matter sciences (Li et al., 2011) (Liu et al., 2015) (Genzer and Groenewold, 2006).
Such pore opening phenomena during drying have also been observed during drying of other soft porous materials (Nguyen et al., 2020).Such phenomena are investigated theoretically (Wang andCai, 2015a, 2015b), where a constrained hydrogel with central cavity is subject to dehydration.If the outer skin is sufficiently rigid, growth of the central cavity is observed.Comparable physics are studied in (Curatolo et al., 2018) (Curatolo et al., 2023), which also regard a hydrogel with a central cavity subject to dehydration, but where the cavity remains filled with water.
Similar phenomena also play a role in biology.At excessive stresses, the liquid pressure in the xylem of vascular plants can become negative, leading to embolism, i.e. cavitation of gas bubbles (Konrad and Roth-Nebelsick, 2003) (Stroock et al., 2014).These conditions have been investigated experimentally using gels under tension (Vincent et al., 2014) (Wang et al., 2023).Cavitation and pore growth is a means to relax the tension (Vincent et al., 2014), but the growth can happen too fast under unstable conditions (Wang et al., 2023) (Sida et al., 2023).Gas bubbles can not exist under these conditions.It will grow unboundedly.In practice these conditions will lead to bursting or fracturing, rendering access of the bubble to the environment (Wang et al., 2023).This instability of gas-filled pores under tension is used by ferns for spores dispersal (Kang et al., 2017) (Kang et al., 2021).The physical description of this instability would involve accounting for inertia, and fracture properties of the biological tissue.However, such extreme mechanics are not observed during the drying of food materials.
In this paper, we examine the transient problem of pore development during dehydration using a simplified system, capturing the essentials of drying food materials.As an idealization of drying vegetable/fruit, we assume a sphere with an elastic skin and a core made of hydrogel.In the centre of the sphere is a small cavity, representing an intercellular pore.Hydrogel material is a good approximation for vegetable/fruit tissue, as our earlier study has shown that their water-holding capacity is well described by the classical Flory-Rehner theory, which is commonly applied to swelling of hydrogels (van der Sman et al., 2013).
In a previous paper we have developed for a similar problem a Finite Volume implementation of the model, where the large deformation mechanics is coupled to moisture transport, cf.(Bertrand et al., 2016;Jin and van der Sman, 2022).With this paper, we want to validate our Finite Volume approach via comparison to a Finite Element solution and an analytical approximation to the dehydrating pore/core/shell problem.Problems with coupled mechanics and moisture transport are commonly solved with the Finite Element method using the weak formulation.However, in the field of food science, the weak formulation is rarely used because of the involved mathematical technicalities.Hence, this paper will also have a strong tutorial character.Consequently, in the appendix, we give an elaborate introduction to the model description using the weak formulation.
The Finite Volume scheme is based on a reference frame, co-moving with the solid phase (polymer network).The Finite Element solution is based on a Lagrangian description of the problem, with the reference frame equal to the initial, stress-free free-swollen state, which is solved by COMSOL.Furthermore, we give an analytical approximated solution to this problem under the assumption of large pores and negligible gradients in the moisture content.After comparison of the three solutions, we will perform an extensive parameter study, whose results we use to discuss the physics of the problem.We conclude with a discussion of the advantages of the different approaches.

State variables & kinematics
In solid mechanics it is custom to describe the large deformation with respect to a non-moving reference frame B ref (Hong et al., 2008).The position of a material point in this reference frame B ref is given by X, but due to deformation the material point has moved to a new position in the current frame (B τ ): x = f(X, t) = X + u(X, t), with u the displacement vector.Via the deformation gradient tensor F = ∇f = I + ∇u (with I the identity tensor), one can compute how geometric elements transform from the reference frame B ref to the current frame B τ .We consider the transformation of a volume element dv ref , a surface element da ref , and a normal vector m = n ref from the reference frame B ref to those (dv, da, n) in the current frame B τ : with the Jacobian J = det F, and the conjugate F* = J F − T .Below, we will use these relations for transformations of the governing equations from the current frame to the reference frame.In solid mechanics it is custom to solve the problem in the reference frame.Also, we define the stretch parameters λ i , with λ 2 i the eigenvalues of the Cauchy tensor C = F T F.
We like to note that in literature there is still a debate what is the correct reference state for polymer gels (Urayama et al., 1996) (Quesada-Pé et al., 2011) (Sakumichi et al., 2022).Often for synthetic polymers the dry state (without solvent) is taken as the reference state (Hong et al., 2009).However, in polymer physics, it is argued that the polymer fraction at preparation (with solvent present) must be taken as the reference configuration.Yet, food gels have often physical crosslinks, and consequently, there is no clearly defined preparation state.
Because we have a bilayer material (core-shell material), which is stress-free and uniformly swollen in the initial state, it can be shown that a dry configuration with zero stress and uniform deformation is physically not realizable.Hence, for our problem, it is more convenient to use the initial, free swollen configuration B 0 as the reference frame.The system attains its free swollen configuration in equilibrium with pure water, with water activity a w = 1.In the free swollen configuration the system has uniform deformation and zero stress.
Yet, (food) hydrogels the biopolymers in the gel network are elastically stretched in the initial free swollen configuration state.Consequently, in polymer physics one takes as the reference frame the condition, where the biopolymers are unstretched.In the initial state the biopolymers have an uniform, isotropic stretch λ i = λ 0 .In our earlier paper (RGM Van der Sman, 2015) we have shown that for many food materials there is a universal value J 0 = λ 3 0 = 3/2, independent on the degree of crosslinking.There we have defined the polymer volume fraction in the free-swollen configuration as φ 0 , and the polymer volume fraction in the unstretched configuration as φ ref = λ 3 0 φ 0 .As our system is a bilayer core/shell system with each layer having different crosslink density, we think it is convenient to define the deformation with the respect to take the initial free-swollen configuration B 0 , but define the elastic deformation with respect to the unstretched configuration B u .Hence, the elastic deformation F e is defined following a multiplicative decomposition of F: Thus in the initial state F --I, and F e = λ 0 I, and thus it follows that F u = λ − 1 0 I.A similar multiplicative decomposition of F is commonplace in the field of plastic deformation, where the total deformation has both inelastic (plastic) and elastic contributions, with the latter defined with respect to an intermediate stress-free configuration (Simo, 1988).In our R.G.M. van der Sman et al. problem, the unstretch configuration frame serves as our intermediate configuration, from which we compute the elastic deformation.Note, that core and shell have different unstretched reference configurations, as they differ in φ ref .This multiplicative decomposition and its connection to the three defined coordinate frames is nicely illustrated in Fig. 1.Thus, we have the following transformations.
• F maps the deformation of material points from the initial configuration B 0 to the current configuration B τ , with det F = J = φ 0 /φ, • F u maps from the initial configuration B 0 to the intermediate unstretched configuration B u , with det F u = J u = φ 0 /φ ref , and • the elastic deformation F e maps from the intermediate unstretched configuration B u to the current configuration B τ , with detF e = J e = 1/φ = φ ref /φ, Thus, instead of the immediate mapping F, one can do this also in two steps: via F u , and F e .The stress will be a function of the elastic deformation, which will be computed as: Note, φ we have used earlier in our papers about Flory-Rehner theory (RGM Van der Sman, 2015).

Spherical symmetry
Under the assumptions of spherical symmetry, a material point X is uniquely determined by the radial coordinate R, and the only non trivial component of the displacement u is the radial displacement u = u(R, t); the current position x is represented by r = R + u(R, t).The deformation gradient F is given by: where λ r , λ θ are the radial and the hoop stretches measured with respect to the free swollen configuration B o .Also the composition of deformations and their determinants simplifies as:

Balance laws & constitutive relations
In this section, we describe the model for the osmotic dehydration and deformation of a spherical hydrogel, with a gas-filled cavity, and a permeable elastic skin.The elastic skin can also absorb water, but it has a higher crosslink density than the core.Hence, the two domains of our core-shell system have the following reference values (φ ref ): where G skin and G core is the shear modulus of the two regions.The value of φ skin follows from the c*-theorem of deGennes (RGM Van der Sman, 2015).

Mass balance of water
First, we describe the balance equation of liquid in the current frame, that is using the spatial fields as φ w defined on B τ .Let D t = ∂ t + ∇ ⋅v be the material derivative, with v the spatial velocity defined by the time derivative of the deformation field u.The mass balance for the water in the core and gel is: where j w is the diffusive water flux given by the generalized Darcy's law, D s is the self-diffusivity of water, and μ w is the chemical potential (in units of [J/m 3 ]).The chemical potential is derived from a free energy functional (Hong et al., 2009).It can be decomposed into two contributions: p liq is the hydrostatic pressure in the solvent, which is a consequence of the incompressibility of the hydrogel material.μ w,mix is the so-called mixing contribution, and it can be derived from Flory-Huggins theory for polymer gels.
However, it is convenient to reformulate the mixing contribution into an osmotic pressure Π mix , which follows the scaling law of Cloizeaux (RGM Van der Sman, 2015).Earlier, we have shown that the response of food gel materials under osmotic pressure can be mapped to a single master curve (RGM Van der Sman, 2015), which is a consequence of the c* theorem of deGennes (RGM Van der Sman, 2015).From this scaling behaviour it turned out that the mixing pressure can be reformulated in terms of the elastic modulus of the material G: and α is a constant, defined in (RGM Van der Sman, 2015).

Balance of forces on B τ
We assume that inertia forces are negligible and external forces are null; the balance of forces in the current configuration writes as: where σ is the actual stress (Cauchy stress), n the normal to the boundary ∂B τ , and t the force at boundary.Under spherical symmetry hypothesis, the non zero components of the stress are σ rr , and σ θθ .
By assuming the Neo-Hookean model, and setting φ = φ/φ 0 λ 0 , the stress components are given by (Van der Sman, 2015): Using φ/φ 0 = 1/λ 2 θ λ r the stress components are: From ( 9), it follows: The initial configuration is used as the domain frame, defining the deformation F to the current (deformed) configuration.The elastic deformation F e is defined via the reference configuration.Note, that in all 3 configurations, the skin and core have matching deformations at their interface.
From the difference between radial and hoop stresses, we can eliminate the pressure: The condition for mechanical equilibrium becomes: This equation will be integrated numerically to obtain σ rr (r).Subsequently, we can compute the pressure field p liq via: which can be inserted in the definition of the water chemical potential, which enters the mass flux j w .The problem is complemented by the boundary conditions.At the outer boundary the stress balances the external pressure p 0 : Note, by convention the stress and pressure have opposite signs.At the inner boundary at the wall of the cavity holds: with p gas the gas pressure in the central cavity, which will be derived below.The boundary condition for the chemical potential at the outer boundary is Π ext is the osmotic pressure of the external fluid.
The central cavity is assumed to be filled with both air (inert nitrogen) and water vapour.We assume that the air does not dissolve in the hydrogel, but the water vapour can be exchanged with the hydrogel.
We assume local equilibrium for the gas bubble in the cavity.In the free-swollen initial state it holds that p gas = p 0 , and the vapour in the cavity is saturated (a w = 1)p vap,0 = p sat (T).Thus the amount of air is p air,0 = p 0 − p vap,0 , and N air,0 = p air,0 V gas,0 /RT, with the initial gas volume V gas,0 = 4πr 3 0 /3.This amount of air, N air,0 = N air remains constant throughout the drying.
The equilibrium condition for the cavity is: with p gas = p air + p vap , and p air = N air RT/V gas , using V gas = 4πr 3 in / 3 for the gas bubble volume.

Transient solution
With the Finite Volume we solve the problem in the current frame B τ , which moves along with the biopolymer network.The (current) domain is subdivided into control volumes being spherical shells, with same thickness smaller than the outer radius: Δr ≪ r out .The boundaries of the different domains, i.e. a) the external surface, b) the internal surface of the cavity, and c) the interface between core and shell, exactly match with the boundary between control volumes.The state variable φ w represents the volume-averaged value of the control volume and approximates the value at the centre of the control volume.The coordinate system will evolve with the dimensional changes of the control volumes, as follows from the temporal evolution of φ w .The latter will be driven by gradients in μ w , which are also computed at the centers of the control volumes.However, the fluxes J w,r are computed at the boundaries of the control volumes following central differencing: with Δr the distance between centers of the adjacent control volumes.
With the fluxes one can compute the change in the volume of water inside the control volume, using simple Euler forward time integration: with Δr + and Δr − distances to adjacent control volumes.Note, that the volume of the polymers ΔV s in the control volume always stays the same, as the coordinate system is linked to the deforming polymer network.
The polymer volume fraction then reads: Mind, that the chemical potential contains the hydrostatic pressure p liq , as μ w = − Π mix (φ s ) + p liq .p liq follows from solving the momentum balance: ∇ ⋅ σ = 0.This equation will be integrated starting at the outside r = r out .The boundary condition states σ rr (r = r out ) = − p 0 .The radial stress will be computed at the boundaries of the control volumes (with Δr the thickness of the control volume): with the gradient in the radial stress at the centre of the control volume: using this gradient we also compute σ rr (r), via which we compute the hydrostatic pressure: This pressure is substituted in μ w (r).
Via integration of the force balance follows the radial stress at the inner surface of the cavity: σ rr (r in ) = − p gas , determining the gas pressure, and the chemical potential μ w (r in ) = − R gas T log(a w,gas ) + p gas .Thus follows the vapour pressure: p vap = a w,gas p sat (T), and the air pressure p air = p gas − p vap .Subsequently, we compute the gas volume: p air V gas = p air,0 V gas,0 .This finally renders the amount vapour moles in the cavity: N vap (t + Δt) = p vap V gas /R gas T. We use this to compute the molar flux at the inner boundary: As the new gas volume is known, V gas together with volumes of water, ΔV w (r), one can compute the new positions of the boundaries (r ± Δr/2) and centers of the control volumes (r).We note that J w (r in ) represents the amount of water turned into water vapour via evaporation.As we assumed isothermal conditions, the evaporation does not lead to evaporative cooling.In Fig. 2 we have indicated schematically how the mass and momentum balances are solved, and the coupling between them.The divergence of the mass fluxes J w leads to changes in the amount of water (φ w ), which leads to deformation of the sphere -as captured by the deformation gradient F --I + ∇u.The deformation leads to stresses σ ii , which can be viewed as momentum fluxes.The incompressibility constraint leads to the hydrostatic pressure p liq , which couples back to the chemical potential μ w , driving the mass fluxes.

Steady state solution
We will also solve the steady state problem numerically, cf.(Van der Sman, 2015).As a first step, we integrate the momentum balance, given the external osmotic pressure Π ext , and an assumed stretch Λ b , which renders p gas , and Λ a (or equivalently V gas ).Local values of φ for skin and gel is calculated from the chemical equilibrium: μ w (r) = μ w,ext , or rather From the solution of the momentum balance we compute p vap , via p gas = p air + p vap , and p air = N air,0 R gas T/V gas .This value of p vap should match that of chemical equilibrium: p vap = p sat (T) exp[(μ w, ext − p gas )ν w /R gas T].This non-linear equation is solved via the secant method.

Analytical solution
For the steady state solution we can formulate an analytical approximation, cf. ( Van der Sman, 2015).Using the incompressibility condition, and defining λ = λ θ , the condition for mechanical equilibrium can be rewritten as: Following Pence and Tsai (2006) (Van der Sman, 2015) we change the integration variable using the total derivative of λ = r/R: and consequently As a first approximation, we can assume that for thin to moderately thick shells, with r a < r < r b , the swelling is uniform cf.(Pence and Tsai, 2006).
For a homogeneous shell of hydrogel material with elastic modulus G, the integration of the condition for mechanical equilibrium gives the following analytical expression for the pressure drop Δp over the shell: with Λ a = λ(r = r a ), and Λ b = λ(r = r b ) the stretch factor at respectively the inside and outside of the shell.Note that, the above expression is mathematically identical to that derived in (Curatolo et al., 2018), although the dry state is used as the reference state.
An approximation for φ follows from the boundary condition at the outer boundary: Given the swelling Λ b and the swelling ratio from the above relation, we can easily compute the inner radius r a of the shell from the incompressibility condition: with R B and R A the dimensions of the shell in the initial state.
In our problem, we have a hydrogel with a core-shell structure, with elastic moduli G core and G skin .The above expression for the pressure drop over a spherical shell can be applied to both core and shell.Hence, the pressure difference between the cavity and the environment can be decomposed in two contributions: with Λ m is the stretch factor at the interface between core and shell (r = r m ), Λ in and Λ out are the stretch factors at r in and r out respectively.φskin follows from the chemical equilibrium at r = r out , and φcore follows from the chemical equilibrium at r = r m , using that p liq (r = r m ) = p 0 + Δp skin .Note, that in steady state at every location, the chemical potential equals that of the environment: μ w (r) = μ w,ext

Finite Element solution using weak formulation
Following the common approach in solid mechanics for hydrogels (Hong et al., 2008), we solve the problem also via the Finite Element method using the weak formulation with the initial frame B o as the reference configuration.As researchers in food science are often not familiar with the weak form, we give a short explanation how to derive it in the Appendix.Before stating the weak formulations, we reformulate the mass and force balances in terms of the spatial fields defined in the initial frame B o .

Mass balance of water on B o
Above, the mass balance of water has been written on the current configuration B τ by using as state variable the volume fraction φ w , see (6).In Finite Element models one uses as the state variable for the mass balance the molar concentration c of water per volume in the initial frame B o (Hong et al., 2009).This relates to the volume change J due to deformation: where ν w is the water molar-volume.Note, that the volume fraction of water in the current frame B τ equals φ w = 1 − φ s = 1 − φ.Thus, φ w = ν w c/J.Hence, the state variable c also accounts for volume changes due to dehydration, i.e. it is linear with the pull-back of the volume fraction φ w , and it might be used for the incompressibility condition.
The mass balance (6) written on B τ , can be reformulated in the reference frame B o as follows: where h w is the pull-back of the current molar flux j w /ν w , m the normal to the boundary ∂B o , and q the liquid flux at boundary.From the formula relating the spatial gradient ∇μ w (x, t) to the material gradient ∇μ w (X, t): it follows: The mobility tensor M is given by with M the mobility in the current frame.
In spherical coordinates, the mass balance thus reads: The prefactor λ 2 θ /λ r accounts for the deformation of the exchange surface area (∼ λ 2 θ ∼ (r/r 0 ) 2 ) and the deformation of the thickness of control volumes: (~λ r ~ (dr/dr 0 )).

Balance of forces on B o
The balance equation ( 9) written in the current frame B τ can be reformulated in the reference frame B o as follows: (41) where S = σ F* is the reference stress (first Piola stress), and t o = |F* m| t the force at boundary.Under the spherical symmetry hypothesis, using Eq. ( 11), the first Piola stress S simplifies as: The balance of forces (41) rewrites as

Weak formulation
We formulate the problem in the weak form.How to derive the weak forms is shown in the Appendix.We give only the left-hand side of the weak form, as the right-hand side is handled via the boundary conditions.The weak contribution W 1 for the balance of forces reads: where ũ and p are the test functions associated to the displacement u = r − R, and to the Lagrange multiplier p, respectively, and p is the effective pressure at the boundary.Note, that different relations hold for S in core and skin, as these layers differ in their elastic modulus G.
At the boundary we have: Their weak formulation is: The initial conditions are formulated for the displacement.As we use the initial state as the domain configuration, it follows that u = 0.The weak contribution W 2 for the incompressibility condition, where the pressure is represented in terms of a Lagrangian multiplier, reads: with J = λ r λ 2 θ , and Ĵ = ν w c + φ 0 .Note, that φ 0 = φ ref /λ 3 0 is different for core and skin.No boundary conditions are required for the incompressibility condition.The initial condition for the pressure field for the freeswollen state is different for the core and shell: The weak contribution to the mass balance is split into the core contribution W c and the skin contribution W s : where R m is the position of the interface core/skin at initial conditions, q the mass flux at the boundaries of the two subdomains.c s and c c are the state variable for water content in the skin and core respectively.Finally, h w,R is the radial component of the flux h w given by: The boundary conditions are formulated as Robin type relation for the molar fluxes.These conditions are different from the above Finite Volume model, where Dirichlet conditions are assumed for μ w .As the chemical potential is not a state variable in our Finite Element formulation, we found it more convenient to formulate it as a Robin-type boundary condition for the molar fluxes.Moreover, they allow a more gradual adaptation of the core/shell system to the environmental conditions.To ensure continuity of the chemical potential at the interface of core and shell (at R m ), there a Robin-type boundary condition with the flux q linear in the difference between the chemical potentials at either side of the interface.
We assume convective boundary layers in the gas directly surrounding the skin and in the central cavity, allowing for evaporation of liquid water into vapour, and we assign different boundary conditions for the core and the skin; for the core we have: for the skin R.G.M. van der Sman et al. q The mass transfer coefficients are β cav = D vap /(R cav /5), and β air = h air / ρ air c p,air following the Lewis relation.h air ≈ 20 W/m 2 .K for air velocities about 1 m/s.RH cav = p vap /p sat (T) follows from the mass balance for the water vapour in the cavity, with p vap the vapour pressure, and p sat is the saturated vapour pressure.The vapour pressure enters the gas pressure: and p air follows Gay-Lussac: N vap is the state variable for the mass balance of the cavity, written as a 0D balance equation, implemented in COMSOL as an ordinary differential equation (ODE): Finally, the complete weak formulation for the coupled problem writes as: find u, p, c c , c s such that: The problem is supplemented with initial conditions for both core and skin:

Steady state solutions
We first validate our simulation model against the analytical solution.We assume a single domain with G skin = G core , subject the system to different external osmotic pressures Π ext .are performed for R out /R in = 0.6 and T = 308 K. Results are shown in Fig. 3.We have scaled Π ext against, (α − 1)G core which is the osmotic pressure, where the gel would be in the unstretched state, B u , with λ e,i = 1 and φ s = φ ref .The results show that both the Finite Element and the Finite Volume methods show good agreement with the analytical solution, but the Finite Volume method is less accurate in the predictions of the gas pressure, p gas .
Similar simulations are performed for the cases G skin /G core = {2, 5}.Again, we have scaled Π ext with (α − 1)G core , with G skin = G core .As above, we have R out /R in = 0.6 and T = 308 K. Results are shown in Fig. 4. Again, there is reasonably good agreement between the three solution methods.We take this result as a good validation of both the Finite Element and Finite Volume method.
We have performed a parameter study for G skin /G core = {10, 20} and 4 ≤ Π ext /G core ≤ 20, with results shown in Fig. 5.We show there the profiles of Λ θ and Λ r as a function of the position in the material frame R/ R out .Observe the continuity in Λ θ , and the discontinuity in Λ r at the core/shell interface.The solutions of the Finite Volume and Finite Element method coincide very well with each other.The results show that for G skin /G core Λ θ (R in ) = r in /R in > 1 for all values of Pi ext , meaning that in the steady state the pore is larger than in the initial state.For G 1 / G 0 = 10 the pore remains smaller than the initial pore size.
We investigate the pore growth via a parameter study for Π ext /G 0 = 20, while varying G skin /G core to indicate a critical ratio for pore opening.Results are shown in Fig. 6.We observe that for G skin /G core > 10 it occurs that p gas = − σ rr (r in ) < p 0 , meaning that the pore has expanded beyond its initial size R in .This is also observed for the actual position of the pore for G skin /G core > 10, which is right of the vertical dashed line indicating the initial pore size.
We investigate the development of gas pressure as a function of initial pore size, using Π ext /G core = 20, G skin /G core = 20, and t skin = 0.05R out .First, we analyse the principal stretches and stress profiles, which are shown in Fig. 7.We observe that for smaller pores, stronger gradients in σ rr and Λ θ develop at the interface between the core and gas Fig. 3. Steady state solution of a single domain system (G skin = G core ) regarding a) polymer volume fraction φ (left pane) and b) pressure difference between cavity and environment p gas − p 0 , as follows from the analytical solution, the Finite Element and Finite Volume method.bubble.Finite Volume and Finite Element solutions are quite similar.Additional simulations show that Finite Volume simulations converge to the Finite Element solutions with the increase of the number of grid points.However, Finite Volume simulation for the smallest pore becomes unstable.
Gas pressure as a function of pore radius is shown in Fig. 8, as obtained from Finite Element simulations.We observe that with decreasing pore radius, R in , the gas pressure approaches a limiting value, which is Fig. 6.Radial stress σ rr profiles (as function of current position r/R out ) at steady state for different ratios of G skin /G core .Skin thickness is t skin = 0.05R out , and the initial pore radius is R in = 0.6R out , which is indicated with the dashed vertical line.If − σ rr (r in )/p 0 − 1 < 0 the pore has expanded beyond its initial size.This condition is indicated by the solid horizontal line.just above zero.Hence, there is no sign of instability here.We investigate further the case of Π ext /G core = 100, which shows similar behaviour, albeit a slightly lower final pressure.The limiting pressure is about reached if R in /R out ≈ 0.05.
We investigate the performance of the transient Finite Volume solver.We compare its steady-state solution with those of the Finite Volume steady-state solver (using similar discretization), and the Finite Element solver.Simulations are performed for G skin /G core = 20, t skin /R out = 0.05, R in /R out = 0.7, and Π ext /G core = {6, 10}.These values were chosen to keep the computation time of the transient Finite Volume solver small.Results are shown in Fig. 9.We conclude that the three solvers give quite comparable solutions for σ rr , but they are in good agreement for λ θ .
Furthermore, we have compared the transient solution of the gas pressure and pore radius for the Finite Volume and Finite Element solvers, as shown in Fig. 10.Due to some differences in the boundary condition at the outer boundary, there is some difference in the dynamics, but overall there is similar behaviour.Both solutions show an equal expansion of the pore.Surprisingly, the Finite Element solution for the gas pressure follows exactly the analytical solution, Eq. ( 34), as indicated by the black lines.The mass transfer is so slow, that the gas pressure follows a quasi-steady state solution.
The (transient) Finite Element solver has a much shorter computation time than the transient Finite Volume solver, as the latter was based on Python, which is interpreter-based.We assume all solvers to be valid, but given the difference in computation time, further analysis is done with the Finite Element solver.
First, we investigated the evolution of the total and partial gas pressures for a range of external pressures Π ext /G core = {20, 100, 500, 1750}, as shown in Fig. 11.At very short times we observe compression of the central pore, shown by the increase of gas and air pressure.Subsequently, we observe that first, the partial air pressure decreases due to the instantaneous increase of pore size, because the drying front reaches the central pore.Only then the partial vapour pressure decreases, which remains in equilibrium with the hydrogel bounding the pore.At large values of the external pressure Π ext /G core ≥ 100 the air pressure goes through a minimum, before stabilizing at a larger value.Remarkably, the final air pressure exceeds the final vapour pressure at the maximal external pressures Π ext /G core = 1750.The Finite Element solver is unstable for Π ext /G core ≥ 2000, possibly due to physical instability of the central cavity.
In the right pane, one can also observe the evolution of the gel radius r out and pore radius r in with time.Similar to the air pressure, we observe a temporary extremum in the pore radius.This extremum in the radius is explained by the sudden expansion of the pore, which happens faster than the diffusion of moisture in the core.Consequently, it takes time for the vapour pressure to reach a new equilibrium, as shown in the graph for the vapour pressure p vap in the left pane.Due to this extra supplementation of total gas pressure, the radius shrinks back for a small amount.
Furthermore, the sudden increase of r in is faster for increasing external osmotic pressure Π ext .The difference between final values of r in and r out becomes smaller with increasing Π ext , but their ratio goes to an asymptote.
We perform a similar analysis for a higher ratio G skin /G core = 50.Results are shown in Fig. 12.We observe qualitatively similar behaviour as in the case of G skin /G core = 20.However, the shrinkage of the outer radius is much less, and the pore opening is larger.The instability of the code happens at similar conditions Π ext /G core ~ 2000.We have also performed similar simulations at higher temperatures, which should give higher vapour pressures.We have changed temperature to T = 338K (65 o C), with G skin /G core = 50.Results are shown in Fig. 13.Again results are qualitatively similar, albeit with higher vapour pressures.Total gas pressure is mainly determined by the vapour pressure.However, the instability of the scheme is still a similar condition: Π ext / G core ~ 2000.Investigation of how final gas pressures depend on Π ext does not show any indication of an onset of physical instability.Hence, it is likely that the instability is numerical.
Calculations show that the core material will have a w ≈ 0.07 in case of Π ext /G core ~ 2000.At this low water activity at the end of drying most Fig. 11.Left pane: Gas, and partial vapour and air pressures as a function of time, and rightpane: inner (red) and outer radius (blue) of the system as a function of time.Both graphs as from simulations with G skin /G core = 20, t skin = 0.05R out , R in = 0.1R out , Π ext /G core = {20, 100, 500, 1750} (decreasing order of transparency) following the Finite Element solution.Fig. 12. Left pane: Gas, and partial vapour and air pressures as a function of time, and right pane: inner (red) and outer radius (blue) of the system as a function of time.Both graphs as from simulations with T = 308K (35 o C), G skin /G core = 50, t skin = 0.05R out , R in = 0.1R out , Π ext /G core = {400, 800, 1800} (decreasing order of transparency) following the Finite Element solution.Fig. 13.Left pane: Gas, and partial vapour and air pressures as a function of time, and right pane: inner (red) and outer radius (blue) of the system as a function of time.Both graphs as from simulations with elevated temperature T = 338K (65 o C), and G skin /G core = 50, t skin = 0.05R out , R in = 0.1R out , Π ext /G core = {400, 800, 1800} (decreasing order of transparency) following the Finite Element solution.food biopolymers will be in the glassy state, and they will have a high elastic modulus (Van der Sman et al., 2022) (van der Sman et al., 2023), as is also indicated in our discussion section.Hence, our approximation of a constant elastic modulus of the core material is not realistic.Consequently, we think it is not worthwhile to resolve the apparent numerical instability.

Discussion
In this paper, we have shown that substantial pore formation and growth can occur during the drying of soft materials in the presence of a stiff elastic skin.The transient and steady-state solutions of the Finite Volume solver compare quite well with the Finite Element solutions.In the case of relatively large pores, the steady-state solutions also compare well with the analytical solutions.If pore formation is gradual, the gas pressure in the pore also follows the analytical solution -showing that the pore is in a quasi-steady state.
The validity of the analytical solution is restricted to cases where the polymer volume fraction remains relatively uniform in both skin and core materials.The Finite Element solver has a significantly wider range of numerical stability than the Finite Volume solver.However, a Finite Volume implementation can still have its merits, because they are much easier to implement in general-purpose programming languages like C/ C++ and Python.Also, they might be more suited in multiscale simulation, where pore formation is solved at the mesoscale, which is coupled to heat/mass transfer at the macroscale -similar to the multiscale cell model for the expansion of starchy snacks (van der Sman andBroeze, 2014a, 2014b).
We have also investigated the physics of the pore formation in more detail -bas shown in Figs.11-13.In the absence of an elastic skin, we find that the initial pore only undergoes shrinkage.Only, if the skin has sufficient stiffness, G skin /G core > 10, one observes that the pore expands.The smaller the initial pore size, the lower the final gas pressure.The relative growth of the pore r in /R in increases with decreasing initial pore size.The pore expansion problem also has interesting dynamics: after an initial pore compression stage, the pore suddenly expands -which reaches a maximum just before the drying front reaches the inner pore.This moment of drying front reaching the inner pore is indicated by the fact that the vapour pressure reaches the equilibrium value, as determined by the external osmotic pressure Π ext .
In our simulations, we have not observed any physical instability, where the pore expands without bounds.There is a regime where the Finite Element scheme becomes unstable, but extrapolations of gas pressure and pore size into this regime do not hint at the case of physical instability but numerical instability.This numerical instability occurs at relatively extreme conditions, Π ext /G core ≈ 2000, meaning that the soft material is quite dry (a w < 0.1).Food materials are commonly in the glassy state, and our approximation of the core material as a hydrogel with a constant elastic modulus does not hold anymore.Consequently, we think that the Finite Element solver is well applicable to problems of vegetable and fruit drying.The Finite Element model is also tested at elevated drying temperature, which renders elevated vapour pressures, which thus enhance the total gas pressure in the expanding pore.With increasing temperature, we observe an increase in the ratio of r in /R out , which means an increase in porosity.
For our model to approach the realistic food material/vegetable drying process, we need to account for their viscoelastic properties (Ozturk and Takhar, 2020) (Hu et al., 2023), and the change of such properties with temperature and moisture content.Fundamental theories for their viscoelastic properties still need to be developed.However, for proteins we have shown that their elastic modulus changes with T g /T, the ratio of the moisture-dependent glass transition temperature and the actual temperature (van der Sman et al., 2023), and for starch and maltodextrins we have shown that zero-shear viscosity and relaxation times also scale with T g /T (Van der Sman et al., 2022).Accounting for viscoelasticity is important for modelling food drying, as it explains hysteresis of moisture sorption (Meinders and Oliver, 2015) (van der Sman, 2023) (Hu et al., 2023), and thus that the elastic stresses need to be included in the driving force for moisture transport, as also shown in this paper.The increase of elastic modulus and relaxation times with decreasing moisture content explains the formation of a rigid skin/case hardening during drying, despite starting with a material fully in the rubbery state at the start of drying.At sufficient progression of drying the relaxation times for the core tissue become so large that elastic stresses are not relaxed away, and an increase of porosity is observed (Nguyen et al., 2018).Extension of the current model of the spherical system with pore towards viscoelastic properties is presented in a forthcoming paper (van der Sman et al., 2024).
Further extension of the model is envisaged towards expansion of food materials during intensive heating as baking or frying.This requires the coupling of the viscoelastic model towards an energy balance, and making the material properties temperature dependent.In previous papers we have shown that viscoelastic properties of starches and proteins scale with T g /T, the ratio between glass transition and actual temperature.Near and beyond the glass transition, food materials behave as elastic material with very long viscoelastic relaxation times.Prelimenary simulations show strong expansion of the pore leading to increase of the outer radius of the spherical system beyond its initial size.Preconditions for the strong expansion are: a) initial conditions are near the glass transition, and b) the product temperature needs to reach conditions near or beyond the boiling point of water -which renders gas pressures larger than atmospheric pressure.

Conclusion
In this paper, we have shown a physical model that can explain the increase of porosity in food materials during drying.A prerequisite for this phenomenon is the presence of an elastic skin, having a significantly higher elastic modulus than the core material.In our model, we included a pre-existing pore filled with air and water vapour, which gets enlarged due to tensile stresses induced by the hindered shrinkage of the elastic skin and the continued moisture removal from the core via the drying.The moisture transport and large deformation mechanics are strongly coupled via the pressure field entering both the stress and the chemical potential.
A previously developed Finite Volume implementation has been validated against a (COMSOL) Finite Element implementation and an analytical solution.The Finite Element solution is numerically stable over a wider range of parameters.The Finite Element model shows stable pore growth even under severe drying conditions resulting in low water activities a w ≈ 0.1, and elevated temperatures.Like similar models applied to dehydrating hydrogels, the Finite Element model is formulated via the weak form, which is mathematical rather intricate.Because of the rare use of the weak form in food science, we have explained in a tutorial manner how to derive this weak form.We think only this kind of modelling with a multi-physical coupling of large deformations, heat, and mass transfer can explain the microstructural changes in food material during drying, and intensive heating processes.With the latter, we think of the creation of porous foods via baking or frying.The next step in making our model more realistic is the inclusion of temperature, and moisture-dependent viscoelastic properties.

Declaration of competing interest
Authors declare there is no conflict of interest.

Fig. 1 .
Fig. 1.Definition of the various configurations for defining the deformation.The initial configuration is used as the domain frame, defining the deformation F to the current (deformed) configuration.The elastic deformation F e is defined via the reference configuration.Note, that in all 3 configurations, the skin and core have matching deformations at their interface.

Fig. 2 .
Fig. 2. Schematic representation of the computation of mass and momentum balance, as computed via the Finite Volume methods.Details about the coupling between the balances is explained in the text.
54), μ w,c = μ w (c c ) and μ w,s = μ w (c s ) represent the chemical potential at the core side and the skin side of the interface, respectively, and D μ regulates how fast the liquid mass flows through the interface; we have set it to a small number.β cav and β air are the mass transfer coefficients of the two boundary layers, c sat (T) is the saturated molar water vapour concentration (in [mol/m 3 ]) (following the Tetens relation), RH ext and RH cav are the relative humidities of environment and the cavity.The water activities relate to:

Fig. 4 .
Fig. 4.Steady state solution of a core/shell system (with either G skin = 2G core (black lines) or G skin = 5G core (green lines)) regarding a) polymer volume fraction φ (left pane) and b) pressure difference between cavity and environment p gas − p 0 , as follows from the analytical solution, the Finite Element and Finite Volume method.

Fig. 7 .Fig. 8 .
Fig. 7. Principal stretches λ θ and λ r (left pane), and stress profiles σ rr as function of current position r/R out at steady state for different ratios of R in /R out .Skin thickness is t skin = 0.05R out , and G skin /G core = 20 and Π ext /G core = 20.The Finite Volume solution is indicated with solid lines(with dots), and the Finite Element solution is indicated with crosses.

Fig. 9 .
Fig. 9. Comparison of steady-state solution (λ θ and σ rr ) according to a) Finite Element solver (translucent symbols), b) Finite Volume steady state solver (solid lines), and c) Finite Volume transient solver (dashed lines).Simulations are performed for G skin /G core = 20, t skin = 0.05R out , R in = 0.70R out , and two values of Π ext /G core = {6, 10} (indicated by yellow and cyan symbols respectively).

Fig. 10 .
Fig. 10.Comparison of transient solutions of gas pressure p gas and pore size R gas according to a) Finite Element solver (translucent symbols), b) analytical solution (black lines), and c) Finite Volume transient solver (yellow and cyan solid lines).Simulations are performed for G skin /G core = 20, t skin = 0.05R out , R in = 0.70R out , and two values of Π ext /G core = {6, 10} (indicated by yellow and cyan symbols respectively).