A particle based model to simulate microscale morphological changes of plant tissues during drying

Abstract

Fundamental understanding on microscopic physical changes of plant materials is vital to optimize product quality and processing techniques, particularly in food engineering. Although grid-based numerical modelling can assist in this regard, it becomes quite challenging to overcome the inherited complexities of these biological materials especially when such materials undergo critical processing conditions such as drying, where the cellular structure undergoes extreme deformations. In this context, a meshfree particle based model was developed which is fundamentally capable of handling extreme deformations of plant tissues during drying. The model is built by coupling a particle based meshfree technique: Smoothed Particle Hydrodynamics (SPH) and a Discrete Element Method (DEM). Plant cells were initiated as hexagons and aggregated to form a tissue which also accounts for the characteristics of the middle lamella. In each cell, SPH was used to model cell protoplasm and DEM was used to model the cell wall. Drying was incorporated by varying the moisture content, the turgor pressure, and cell wall contraction effects. Compared to the state of the art grid-based microscale plant tissue drying models, the proposed model can be used to simulate tissues under excessive moisture content reductions incorporating cell wall wrinkling. Also, compared to the state of the art SPH–DEM tissue models, the proposed model better replicates real tissues and the cell–cell interactions used ensure efficient computations. Model predictions showed good agreement both qualitatively and quantitatively with experimental findings on dried plant tissues. The proposed modelling approach is fundamentally flexible to study different cellular structures for their microscale morphological changes at dehydration.

---

1. Introduction

Drying is one of the most popular and generally used food preservation techniques which is currently used on about 20% of the world's perishable crops.¹ Since plant food materials can originally contain water even up to 90% by weight,² moisture removal is the key driving factor which causes food preservation during drying. As a result of such extensive moisture removal from the tissue structures of these food materials, they undergo excessive deformations which lead to significant changes of numerous food properties.

Studying the fundamentals which govern such property changes largely benefit the food industry in order to improve the product quality and the process performance. For this purpose, many areas are currently researched and particularly when it comes to material and process related characterization, different experimental studies have identified strong relationships between bulk scale deformations and microscale deformations of food materials during drying.^3–6 Further, such deformations are found to be mainly driven by the moisture content of plant tissue,^3–7 the drying temperature^8–10 and the cell turgor pressure.¹¹ To relate such governing parameters, micro-scale empirical^5,6 and theoretical models¹² are frequently involved. However, numerical modelling and simulation of these material deformations under drying conditions are still a considerable challenge for the research community and therefore we identified this as a research gap to be filled through a novel modelling approach.

In this context, several plant cell numerical models can be found that are mainly based on grid-based techniques such as Finite Element Methods (FEM) or Finite Difference Methods (FDM), which are primarily developed for basic cell micro-mechanical behaviour studies rather than drying.^13–16 Recently, a FEM based microscale tissue drying model is reported which couples water transport phenomenon and cellular deformations.¹⁷ In this particular model, there are several key limitations such as the limited moisture content range that the model can simulate (only 30% of normalized moisture content reduction) and inability to represent the wrinkling effects that the actual cell walls undergo during drying. Further, it appears that their main focus is on the turgor loss based cellular deformations rather than the moisture content based deformations, where the latter one is actually more critical to industrial drying applications. Also the hybrid numerical simulation approach they have involved is fairly complicated to implement. In addition to such models on drying, there are several vertex cell models reported in literature which are specifically developed to study division and reproduction of cells during tissue development.^18,19 Also, there are gel material models used to simulate bulk level deformations of plant leaves during dehydration.²⁰

However, considering the technical challenges and various limitations, these models have limited applicability in the scope of this research. Such fundamental shortcomings are: poor representation of multiphase phenomena existing in plant materials consisting of liquid, solid and gas phases, use of over-simplified continuum approximations for the tissue materials by neglecting the discrete nature of the cellular structure and its microscale mechanisms, inability to model the excessive cellular level deformations (cell wall wrinkling in particular) at excessive moisture content reductions, and limited flexibility to model the multiscale effects of plant material deformations during dying. These shortcomings need to be well addressed in order to get a clearer understanding on the underlying physics involved in real industrial drying processes. When considering these challenges of conventional grid-based modelling approaches, meshfree methods seem to be more applicable since they fundamentally do not use any interconnected grids as in case of FEM and FDM, which are entirely based on grids.²¹

When considering meshfree methods, Smoothed Particle Hydrodynamics (SPH) is a popular particle-based technique, which was originally developed for astrophysical applications.²² The technique discretises a given problem domain as a set of non-interconnected particles (without any mesh generation), where the particles own system properties which can evolve with time. SPH is also particularly adaptive such that the basic formulations can easily be combined with new physics and mechanics in order to solve different problems of interest.²¹ Using these capabilities, SPH has been successfully used for plant cell and tissue models recently by coupling with a Discrete Element Method (DEM) to simulate compression, tension, shear of fresh cells and tissues^23,24 and even cell breakage.²⁵ We have further improved the modelling technique in a way it can be used to simulate dried plant cells.^26–28 In those models we mainly focused on two-dimensional (2-D) single cells and in this work we present how a 2-D tissue structure can be developed with the use of such single cell models along with several improvements, in order to simulate morphological changes of tissues during drying. Also, we believe that, these proposed methods and findings will contribute to the development of meshfree methods itself and their applications towards developing novel modelling techniques particularly applicable to critical fields such as food drying.

This paper begins with the modelling concepts and formulations used in a single cell model. Then, we discuss how the tissue model can be developed using such individual cells. Next, the approach used to simulate tissue dryness states is presented. Then the experimental findings on dried tissues are compared with simulation results for validation of the model and several other simulations are provided along with a sensitivity analysis of several model parameters. Finally, we present key insights drawn from the study along with prospective future improvements.

2. The particle based model

2.1 2-D representation of plant cells in a tissue

A plant tissue was approximated to an aggregate of cylindrical shaped cells as shown in Fig. 1(a), where the top surfaces were taken as a 2-D model representing each cell. To facilitate this approximation, axial deformations of each cell in the Z direction were assumed to be uniform and the Z directional velocity components and plane stresses of the top and bottom surfaces (in the XY plane) were neglected.²⁴ Accordingly, as shown in Fig. 1(b), each of the individual cells were modelled assuming two basic model components: cell fluid and cell wall, such that the fluid hydrodynamic pressure is counter balanced by the cell wall tension as observed in real plant cells.^{7,14,16,29–32} Following a SPH–DEM based fresh cell modelling approach proposed by some researchers,²⁴ the cell fluid was treated as a low-Reynolds number incompressible viscous fluid and modelled with SPH. Further, the cell wall was modelled using a DEM considering the nature of the solid-dominated fibrous cell wall structure. A particle scheme was used to represent both the cell fluid and the cell wall, and governing equations were applied to these particles. In addition to those basic concepts, we further improved both the cell wall model (see Fig. 2) and cell fluid model (see Fig. 3) in our previous works,^26–28,33 so that the single cell model could be used to simulate dried conditions. The sections below summarise these schemes and concepts along with the new improvements done in this work.

Fig. 1 2-D representation of a plant tissue: (a) a plant tissue is approximated to an aggregate of individual cylindrical cells. (b) Top surface of each cylindrical cell is taken as a 2-D representation of the whole cell which consists of two basic components: cell fluid and cell wall. (c) The problem domain is discretized and represented as two sets of particles: SPH particles for the cell fluid and DEM particles for the cell wall. (d) Individual elements of the discretized cell wall.

Fig. 2 Force interactions used in the DEM-based cell wall model: (a) stiff forces. (b) Damping forces. (c) Wall–fluid repulsion forces. (d) Non-bonded wall–wall repulsion forces. (e) Wall–fluid attraction forces. (f) Wall bending stiff forces. (g) Wall contraction forces. (i: fluid particles; j, k & l: wall particles).

Fig. 3 Force interactions used in the SPH-based cell fluid model: (a) pressure forces. (b) Viscous forces. (c) Wall–fluid repulsion forces. (d) Wall–fluid attraction forces (i & i′: fluid particles; j & k: wall particles).

2.2 DEM modelling of the cell wall

The cell wall was assumed to be a visco-elastic solid material and a neo-Hookean elastic model was used with a supplementary viscous term.²⁴ Several additional features were introduced by us to improve the model performance and to account for drying related physical changes.^28,33 As presented in Fig. 2, the cell wall model accounts for seven types of force interactions: cell wall stiff forces (F^e), wall damping forces (F^d), wall–fluid repulsion forces (F^rf), non-bonded wall–wall repulsion forces (F^rw), wall–fluid attraction forces (F^a), forces due to bending stiffness of the cell wall (F^b) and forces to produce contractions of the cell wall during drying (F^c). Thereby, the total force (F_k) on any wall particle k is derived using these force interactions as:

F_k = F^e_kj + F^d_kj + F^rf_ki + F^rw_kl + F^a_ki + F^b_kj + F^c_kj, (1)

where: i is any neighbouring fluid particle, j is any bonded wall particle and l is any non-bonded wall particle. Here, the stiff force on any wall particle k due to any bonded wall particle j is defined using a spring model²⁴ as:where, G is the shear modulus (≈E/3) with E being the Young's modulus of the wall material, Z₀ is the initial cell height, T₀ is the initial cell wall thickness, λ_θ = L/L₀ is the extension ratio of any cell wall element at the current time step, L is the width of the wall element at the current time step (distance between particle k and j) and L₀ is its initial un-deformed width. Further, the parameter α is calculated with β = 0.5 for cylindrical cells as follows:²⁴

In eqn (1), viscous forces acting on any wall particle k due to the neighbouring wall particles j are defined by using a linear dashpot model as:²⁴

F^d_kj = −γv_kj, (4)

where, γ is the cell wall damping constant and v_kj is the velocity of particle k relative to particle j. In eqn (1), the F^rf, F^rw and F^a forces are used to define the wall–fluid interactions and the boundary conditions. To avoid any fluid particle penetrations through the wall, the F^rf forces are used to create repulsive forces between the cell fluid and wall particles which act through the centre of any interacting wall–fluid particle pair of interest, in an equal and opposite manner. Therefore, the repulsion forces F^rf on any wall particle k from any other fluid particle i are defined as:^21,24

F^rf_ki = f^rf_kix_ki, (5)

where, f^rf_ki is the magnitude of the repulsion force and x_ki is the position vector of particle k relative to particle i. The f^rf_ki is defined based on a Lenard-Jones (LJ) force type as:²⁴where, r₀ is the initial gap between the two particles, r_ki is the current gap between them and f^rf₀ is the strength of the LJ contact. Further, in eqn (1), self-penetration of the non-bonded wall–wall particles are prevented using a similar repulsion force: F^rw_kl with a LJ contact strength of f^rw₀. In order to limit cell wall–fluid detachments during drying, wall–fluid attraction force F^a_ki was introduced by us,²⁸ which was also defined using the above LJ interactions with a contact strength of f^a₀.

Next, in order to account for the resistance that a plant cell wall creates when it experiences local bending and wrinkling, bending stiff forces F^b_kj were used and are defined on any wall particle k within the k and j particle pair as:^28,34–36

where, k_b is the cell wall bending stiffness, L is the width of any given wall element at any given time step, θ is the external angle between the particular wall element and the adjacent wall element as shown in Fig. 2 and Δθ is the change of the θ angle during time evolution.

Next, as given in eqn (1), to incorporate experimentally observed cell perimeter reductions,^6,7,37,38 cell wall contraction forces (F^c) were used in the model and are defined as:³³

where, k_wc is the force coefficient of wall contractions, L is the current width of any particular wall element (see Fig. 1(d)), L^′₀ is the width of the wall element at fully turgid condition, a and b are empirical factors, and X/X₀ is the normalized moisture content of the dried cell to be simulated. The a and b factors were set to 0.2 and 0.9 by considering the normalized cell perimeter trends of our experiments³⁷ (Fig. 13(c) in Section 3), where the overall normalized perimeter reduction is about 0.2 when the overall normalized moisture content reduction is about 0.9. Then the appropriate value for k_wc was determined after a series of trial simulations in order to replicate the normalized perimeter observed from our experiments (see Table 1). Furthermore, we hypothesized that the cell wall itself will also undergo moisture reductions during drying and we simply considered such reductions to be proportional to the fluid mass reduction of the cell (see Section 2.5 for details). For fresh cells, the initial cell wall mass was taken as 10% of the fluid mass (see Table 1), which was calculated using the cell volume and fluid density.^24,28

Table 1 Key physical parameters used for the model

Parameter	Value	Source
Initial cell diameter (D0)	150 μm	3
Initial cell height (Z0)	100 μm	24
Wall initial thickness (T0)	6 μm	15
Initial cell fluid mass	1.77 × 10^−9 kg	Calculated^27
Wall mass (10% of cell fluid mass)	1.77 × 10^−10 kg	Calculated^24
Fluid viscosity (μ)	0.1 Pa s	23, 24
Initial fluid density (ρ0)	1000 kg m^−3	Set
Fresh cell turgor pressure (PT)	200 kPa	14, 24
Fresh cell osmotic potential (Π)	−200 kPa	(Π = −PT)^23,24
Wall permeability (LP)	2.5 × 10^−6 m^2 N^−1 s	Set
Wall shear modulus (G)	18 MPa	15, 24
Wall bending stiffness (kb)	1 × 10^−12 N m rad^−1	Set
Wall damping ratio (γ)	5 × 10^−6 N m^−1 s	Set^24
Fluid compression modulus (K)	20 MPa	24
Wall contraction force coefficient (kwc)	4 × 10^4 N m^−1	Set
LJ contact strength for wall–fluid repulsions (f^rf0)	1 × 10^−12 N m^−1	Set
LJ contact strength for wall–wall repulsions (f^rw0)	1 × 10^−12 N m^−1	Set
LJ contact strength for wall–fluid attractions (f^a0)	2 × 10^−12 N m^−1	Set
LJ contact strength for cell–cell repulsions (f^rc0)	1 × 10^−10 N m^−1	Set
Pectin layer stiffness (kpectin)	20 N m^−1	Set
Pectin layer thickness (Tp)	8 μm	Set
Initial smoothing length (h0)	1.2 × initial fluid grid spacing	Set
Time step (Δt)	2 × 10^−9 s	Set

---

2.3 SPH modelling of the cell fluid

Cell fluid was modelled using SPH by treating it as a high viscosity incompressible Newtonian fluid with low Reynolds number flow characteristics.^23,24,39 Additionally, the cell fluid–wall attraction forces were used to ensure the cell fluid is well attached to the cell wall when simulating dried cells.²⁷ As shown in Fig. 3, the cell fluid model involves four types of force interactions: fluid pressure forces (F^p), fluid viscous forces (F^v), wall–fluid repulsion forces (F^rw) and wall–fluid attraction forces (F^a). The resultant of these forces: F_i on any fluid particle i is defined as:

F_i = F^p_ii′ + F^v_ii′ + F^rw_ik + F^a_ik. (9)

Here, the momentum equation of the standard Lagrangian type SPH formulation that is used to model weakly compressible low Reynold's number fluid flows^24,40 was used to define the pressure forces (F^p) and viscous forces (F^v). Accordingly, for any given fluid particle i, F^p and F^v are defined using the properties of the neighbouring fluid particles i′ as:where at any given time, m, P, ρ, μ, Z and W are the particle mass, pressure, density, dynamic viscosity, cell height and the smoothing kernel. For the smoothing kernel W, the quartic smoothing kernel was used which is calculated for any given fluid particle i as:²¹where, h is the smoothing length at the current time step, S is the ratio of r_ii′/h and r_ii′ is the gap between particle i and any surrounding fluid particle i′ within the influence domain of the particle i (0 ≤ S ≤ 2). In order to improve computational accuracy, the h should be such that approximately 20 particles are maintained within the influence domain of any fluid particle i of this 2-D model.²⁸ Considering the cellular shrinkage, to maintain the above condition on particle number, the smoothing length h was dynamically adopted using a geometrical relationship as:²⁸where, D is the average cell ferret diameter at the current time step, D₀ is the initial cell diameter and h₀ is the initial smoothing length (see Table 1). During the time evolution of the model, fluid particle pressure was evolved using a standard SPH state equation:^21,24where, P_T is the initial cell turgor pressure, K is the fluid compression modulus that ensures the fluid behaves in a fairly incompressible manner,²⁴ρ_i is the current density of each fluid particle, and ρ₀ is its initial density assumed to be equal to the density of water. Next, the density of any fluid particle i was evolved as:²⁴

The first term in eqn (15) accounts for slight density changes of the cell fluid as the cell deforms in XY plane and ρ^*_i is the 2-D density of fluid particle i defined as ρ^*_i = Zρ_i. The ρ^*_i fluctuations are defined using the standard SPH continuity equation as:

The second term in eqn (15) ensures the density evolution compensates for any cell volume changes due to cell height variations and is defined as²⁴:

where, at any given time, Z_t+Δt and Z_t are the cell heights at the current and previous time steps, and Δt is the time step size. Cell height used in these equations were updated in each time step by considering the incompressibility property of the cell wall material as:²⁴

Z = (αλ_θ)Z₀. (18)

The third term in eqn (15) accounts for the slight density changes within the SPH scheme as a result of the cell fluid mass transfer through the semi-permeable cell wall and is defined as:^24,32

where A_c, L_p, n_f and Π represent the total surface area of the cylindrical cell at any given time, cell wall permeability, total number of fluid particles used and the osmotic potential of the cell fluid at a given dried cell state, respectively (see Section 2.5 for details). In eqn (9), the remaining two terms are the repulsion forces F^rw_ik and attraction forces F^a_ik that act on a given fluid particle which are defined using the LJ force type as:

2.4 Modelling of a tissue

In the literature, researchers initially developed a tissue model by aggregating SPH–DEM based individual cells which was used primarily to study mechanical responses of fresh tissues to external loading conditions.^23,24 The model is based on a simple concept and can be implemented fairly easily within a particle scheme. However, the model inherits several conceptual and computational shortcomings which led us to develop a novel tissue model. When considering the limitations of that model, firstly, it is based on a brick-like cell arrangement where each individual cell is initiated as a square in 2-D. As a result, the contacts between adjacent cells are not isotropic and the overall tissue does not fairly replicate an actual tissue. Further, their tissue model does not involve the actual thickness of the middle lamella (model is initiated such that the neighbouring cells have coinciding cell walls). However, they have hypothetically incorporated a middle lamella effect by using a stiff force interaction between cell–cell contacts, incorporating a spring model with a stiffness replicating the stiffness of pectic materials. This method of middle lamella treatment has eventually resulted in coincident steady state cell wall orientations in their model, which we believe are not quite realistic. Furthermore, the LJ repulsion force concept that they have used in order to avoid adjacent cell penetrations during model evolution is computationally costly. There, just to prevent cell–cell penetrations using cell–cell repulsive LJ forces, they have used the full array of fluid particles of each interacting cell, which is computationally not efficient. Considering all these limitations, we developed a new tissue model which is detailed next.

Firstly as shown in Fig. 4(a) and (b), 2-D cells in the tissue were initiated as hexagons which closely resemble the commonly observed honeycomb arrangement of real cells in tissues. This arrangement also ensures an isotropic cell–cell contact in contrast to the brick-like arrangement of the previous authors. Additionally, a positive initial gap is kept between adjacent cells representing the middle lamella of real tissues. Now, as shown in Fig. 4(c) and (d), the cell–cell contact can be defined using two force interactions: stiff forces due to the pectic materials in the middle lamella (i.e. pectin layer stiffness) and a LJ type repulsion forces in order to avoid interpenetrations of adjacent cells. The pectin layer stiff force was defined as a linear spring model²³ acting between the initially adjacent cell wall particles of any two adjacent cells. When the gap between such two interacting wall particles differ from the initially set pectin layer gap, the stiff force always act as a one-to-one contact-force only between the corresponding two particles. Accordingly, the pectin layer stiff force on any wall particle k due to its initial neighbouring particle m on the adjacent cell wall is defined as:

F^e_pectin_km = −k_pectinΔx_km, (22)

where k_pectin is the pectin layer stiffness and Δx_km is the gap difference of the two particles compared to their initial gap. This force helps to maintain the gap between the wall particle pair equal to the initially set pectin layer thickness. Further, this is the only force acting in between cells if they try to separate each other beyond the initial pectin layer gap.

Fig. 4 Tissue model and cell–cell force interactions: (a) hexagonal shaped cells are used for tissue initialization with a positive pectin layer gap. (b) Interacting wall particle pairs of adjacent cells. (c) Pectin layer stiff forces. (d) Cell–cell repulsion forces. (i: fluid particles; k & m: wall particles).

In case if the interacting cells become closer, pectin stiffness creates a repulsion force to separate the cells and tries to return them back to their initial relative positions. The intensity of this force is usually insufficient to fully prevent the cells become interpenetrated. Therefore, a LJ type force is used for this purpose.²³ However, compared to the approach of the previous researchers,²³ we propose a computationally far better LJ force field in this work. Here, the cells are repulsed only using wall particles rather than using the full array of fluid particles. Since the computational cost for SPH–DEM calculations increase nonlinearly with the particle number,²⁸ our proposed method should be computationally superior since the size of the involved wall particle array is far smaller compared to the cell fluid particle array. When implementing this LJ force field, the initial particle gap (r₀) used for LJ force calculations is set to half of the initial pectin layer gap (Fig. 4(d)). Accordingly, when the particle m approaches particle k, until their gap is equal to half of the initial pectin layer gap, only the pectin layer stiff force acts on the two particles to separate them. However, in case if the two particles become further closer (beyond half of the pectin layer thickness), the LJ forces start to act, since the r₀ is set equal to half of the initial pectin layer gap. Thereby, since the LJ force is much stronger compared to the pectin layer repulsive stiff force, when the cells are about to penetrate, the LJ forces become stronger and tends to effectively separate the two approaching particles. Accordingly, the cell–cell repulsive LJ forces are defined as:

where, f^rc_km is the strength of the LJ force field and x_km is the position vector of particle k relative to particle m. Here, the f^rc_km is defined similar to that of the cell wall LJ force field. It should be noted that, in order to effectively prevent cell–cell penetrations, f^rc_km is set comparatively higher than that of the other LJ force fields used in the cell wall model (see Table 1). Using these two force interactions, individual cells were aggregated and the tissue model was developed. The sections below describe how the model was computationally setup and time evolved.

2.5 Computational setup and time evolution of the model

Firstly, in the case of a single cell, fluid particles were placed inside a hexagonal cell wall. The hexagon size was set such that its surface area is equal to the surface area of a circle equal to the cell diameter given in Table 1. Then the model was time evolved using a Leapfrog integrator²¹ with a time step defined using the Courant–Friedrichs–Lewy (CFL)^21,41 criteria. For the wall, a set of virtual particles were also used in order to avoid any undesirable fluid particle penetrations through the wall boundary and also to limit unrealistic fluid–wall separations.^24,28 During fresh cell simulations, according to eqn (19), depending on the difference of the cell turgor pressure and the magnitude of the osmotic potential, the cell fluid mass varies in a similar way as observed in real cells that have semi-permeable cell walls. These fluid mass fluctuations result in slight fluid density variations as defined by eqn (15), which eventually cause significant turgor pressure fluctuations as defined by eqn (14). These turgor pressure fluctuations cause the cell wall to move inwards or outwards which produce cellular geometric changes. Now, as a result of such dimensional changes of the cell, the turgor pressure fluctuates again and it leads to cell fluid mass gains or losses according to eqn (19). This cycle of model evolution repeats until the cell reaches steady state condition where the cell fluid turgor pressure reaches the magnitude of the initially set osmotic potential, and thereafter any further moisture exchanges through the cell wall cease. By now, the cell shape is fairly circular which fairly replicates the top surface of a turgid cylindrical cell. This steady state particle arrangement is considered as the fresh cell state and is used to evaluate final cell fluid mass, turgor pressure and cell dimensions which are eventually used to characterise the fresh cell state in the analysis described in Section 4.

For single cell dried state simulations, instead of a time-domain based approach, we used a moisture-content-domain based approach which was introduced in our previous works.²⁸ Further, based on the experimental evidence of moisture content dependent positive turgor pressure existence in cells when they are subjected moisture deficiencies,^3,38,42–44 we hypothesised that the turgor pressure would remain positive during drying and simply assumed that it would linearly reduce with the reduction of the cell fluid mass.³³ So, for example, when 200 kPa is used as the fresh cell turgor pressure, the dried cells of: X/X₀ = 0.8, X/X₀ = 0.6, X/X₀ = 0.4 and X/X₀ = 0.3 were simulated by initially setting the turgor pressures to: 160 kPa, 120 kPa, 80 kPa and 60 kPa. At the same time, the magnitude of the osmotic potential was also set equal to the corresponding initial turgor pressure in each case and kept constant during the model evolution. As a result, when the model is time evolved as explained above, the final steady state turgor pressure of each dried cell state reaches the corresponding initial turgor pressure value. In addition to this, the cell wall drying effects were also incorporated in the model such that, for each moisture content case mentioned above, the initial cell wall mass was set to 0.8, 0.6, 0.4 and 0.3 of the cell wall mass corresponding to the fresh condition.³³ During the time evolution, the cell wall mass was kept constant at these initially set values. After time evolution, the final cell physical properties and other geometrical parameters were used to characterise each of the dried cell states and compared with the experimental results.

The above mentioned single cell model was used as a building block and the tissue model was setup by aggregating such cells according to the methodology presented in Section 2.4, with the model parameters given in Table 1. For each cell, the corresponding osmotic potential, initial turgor pressure and initial moisture content are set identically and allowed to undergo simultaneous inflation processes according to the time evolution cycle explained above. Since the cells are bonded by the pectin stiff forces, now the steady state 2-D cell shapes becomes quite closer to inflated hexagonal shapes, which reasonably replicate real cells in tissues. As a hexagonal initial shape was used for the cells, fundamental tissue models were produced by simply aggregating 7, 19 and 37 cells (see Section 4 for more details). Further, in order to ensure stability of the model, each simulation run was performed with an initial particle settling phase that runs for about a quarter of the total simulation time, in which the individual cells in the tissue are evolved by temporarily omitting cell wall bending forces and the cell wall contraction forces. The bending force omission is particularly important since the initial cell hexagonal shape has got sharp edges and straight sides which lead to higher initial bending forces that can result in model instability during initial time evolution. Additionally, the cell wall contraction forces can also restrict smooth initial evolution of the model by applying excessive force fields. Therefore, the initial settling phase was involved in order to ensure the time evaluation of the model starts smoothly. Further, an optimum particle scheme was used for each of the cells by considering the model consistency and computational cost according to the method proposed in our previous work.²⁸ Accordingly, 96 wall particles were used and the corresponding fluid particle number was 656 to ensure equal spacing between wall–wall and fluid–fluid particles, when the fluid particles are placed in a square grid. The selected particle scheme was such that the steady state model consistency errors were within 3% and density fluctuations were within 0.1%. These values are in good agreement with the state of the art SPH–DEM pant cell models reported in literature.²⁴ The model was programmed in a parallel C++ code and all the simulations were run on a high performance computer (HPC). The C++ source code was developed partly by referring to an existing FORTRAN based SPH source code.²¹ Model visualizations were performed using Open Visualization Tool (OVITO).⁴⁵

In each case, at the end of time evolution, the dry basis moisture content X (=kg_water/kg_{dry material}) was computed and related with steady state cell deformations in the tissue which were quantified using a set of cellular geometrical parameters: cell area (A), ferret diameter† (D), perimeter (P), roundness‡ (R), elongation§ (EL) and compactness¶ (C). Further, in order to facilitate analysis and comparison of the model performances, normalized parameters (X/X₀, A/A₀, D/D₀, P/P₀, R/R₀, EL/EL₀ and C/C₀) were finally used. For model validation, these results were compared with experimental data for apple cells obtained from our experiments³⁷ (see Section 3) and literature.⁶

3. Experimental studies

In order to investigate the real cellular deformations experienced by plant tissues during drying, a series of experiments were conducted by us on the Galla Apple variety (Malus ‘Domestica’).³⁷ The samples were obtained from a local supermarket in Brisbane – Australia and the initial wet basis moisture content (kg_water/kg_{wet solid}) was measured to be 0.84 ± 0.01. Firstly, ring shaped slices were cut from the middle parenchyma region of the fruit such that the inner diameter was 25 mm and outer diameter was 60 mm. Prior to the drying phase conducted in a hot air oven dryer, a 3% citric acid solution was used to pre-treat the slices for about 1 minute. The dryer was set to maintain a temperature of 70 °C with a constant air velocity of 1.5 m s⁻¹. Then, dried tissue samples corresponding to 30, 60, 90, 120, 150, 180 and 210 minutes of dying were obtained by drying independent samples individually. After each of the drying experiments, one half of obtained samples was used to quantify the dry basis moisture content and the other half was used for microstructural investigations. To determine the dry basis moisture content, weights of the oven-dried samples were recorded and those samples were further dried for about 18 hours at 100 °C in order to obtain the complete dry mass. These mass values were then used to calculate the dry basis moisture content of each sample. All the experiments were repeated once to obtain average results.

Samples used for microstructural investigations were firstly cut in to cubic specimens of 10 mm × 5 mm × 5 mm. These were initially coated with a fixative solution containing 2.5% glutaraldehyde, 4% paraformaldehyde and 0.1 M cacodylate buffer (pH 7.2) and then kept stored for about 12 hours at 4 °C. A 0.1 M cacodylate buffer (pH 7.4) was used to rinse the samples before post fixing with cacodylate-buffered 1% osmium tetroxide at room temperature for 4 hours. Then, ethanol solutions of incremental concentrations of 50%, 70% and 90% were used twice for 10 minutes in each case in order to dehydrate the samples. These were finally dehydrated with a 100% ethanol solution for 10 minutes. Next, a Critical Point Drying^46,47 apparatus was used to dry these dehydrated samples and it was conducted for 30 minutes, twice. To obtain fresh-open cross sections for microscopic examinations, these samples were fractured by freezing in liquid Nitrogen. This method of fracturing helps to prevent mechanical damages to the tissue surface which can disrupt clear viewing of the microstructure. These samples were then mounted on metal stubs and sputter coated with gold up to 10 μm using an automated sputter coater. These were then examined using a FEI Quanta 200 Environmental Scanning Electron Microscope (SEM) at 20 kV and imaged at 200× magnification with an image size of 885 × 1022 pixels. These images were processed by ImageJ software (version 1.46) and the cell-wise average area, ferret diameter, perimeter, roundness, elongation and compactness were quantified. These cellular geometrical variations were then analyzed against the variation of the dry basis moisture content of the samples. The model was validated by comparing its predictions against these findings.

As can be seen from Fig. 5, the tissue morphology strongly depends on the material moisture content. When the fresh tissue is considered as shown in Fig. 5(a), the cells are fairly circular compared to the dried states. This is mainly due to the presence of higher turgor pressure and moisture content in the cells, which apply higher forces on cell walls to make them stretched and thereby undergoing minimum wrinkling.^7,14,16,32 When the material gets gradually dried, since the moisture content and the turgor pressure reduces considerably, cells become shrinked and cell walls experience a higher degree of wrinkling or warping (Fig. 5(b) and (c)). Further, it is observed that the deformations are irregular and the cell wall structure undergoes minimum damage or breakage.^4,7,8 In summary, these findings imply that the numerical modelling of plant tissues during drying should particularly account for excessive moisture removal and cell wall wrinkling. However, according to the best of our knowledge, even the most recent FEM based plant tissue drying model¹⁷ suffers from both of these capabilities, which should be due to the fundamental limitations of the grid-based techniques involved. Section below describes how the meshfree model proposed in this work was used to simulate plant tissues by overcoming such limitations.

Fig. 5 SEM images of apple cells at different states of dryness: (a) X/X₀ = 1.0. (b) X/X₀ = 0.5. (c) X/X₀ = 0.2.

4. Results and discussion

4.1 Simulation of deformations during drying: single cell

Fig. 6 shows the fundamental model behaviour of a single cell in a tissue at fresh and dry conditions. Fig. 6(a) shows the initial hexagonal single cell model used and Fig. 6(b) corresponds to a fully turgid fresh cell. The dried cell conditions are presented in Fig. 6(c)–(f). When comparing Fig. 6(a) and (b), the inflation and size increment of the fresh single cell can be clearly observed. Further, the dried cell conditions indicate a clear reduction of the cell size which fundamentally agrees with the cellular shrinkage observed in plant cells during drying. This model behaviour is mainly influenced by three factors: lower turgor pressure, lower cell fluid and wall mass and higher cell wall contraction forces. Since the turgor pressure remains positive during drying (see Section 2.5 for details), a fairly circular cell shape is maintained throughout the drying phases. Also since the cell is not externally restricted by any other cells or intercellular influences, this circular shape characteristic is further guaranteed (see Section 4.2–4.3 for intercellular influence in bonded cells in tissues).

Fig. 6 Single cell model: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

4.2 Simulation of deformations during drying: dense tissue (7 cells)

Fig. 7 shows the initial particle placement, fresh tissue condition and dried tissue conditions of the 7-cell tissue. Fig. 8 shows enlarged views of these 7-cell tissue configurations. From Fig. 8(a), the initial pectin layer gap is clearly observed. When the fresh tissue condition is considered as seen from Fig. 7(b) and 8(b), the cell shape has differed to a hexagonal shape from the previously mentioned full-circular shape observed in the single cell fresh condition. This is mainly due to the intercellular interactions. Further, it is observed that the pectin layer is stretched particularly at the cell wall positions which correspond to initial hexagon edges. This is mainly due to the resistance created by the pectin layer for the fresh cell inflation that occurs at a higher turgor pressure. When considering the dried tissue states from both Fig. 7 and 8, a gradual cellular shrinkage pattern is clearly observed. Further, as the tissue gets further dried, the pectin layer tends to experience lower stretch conditions compared to the moist samples. This indicates that the tissue is more prone to deformations mainly due to the lower moisture content and turgor pressure.

Fig. 7 7-cell tissue model: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 8 7-cell tissue model (enlarged view): (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

4.3 Simulation of deformations during drying: dense tissue (19, 37 cells)

When considering the 19-cell and 37-cell tissues as seen from Fig. 9–12, similar shrinkage patterns can be observed during drying. But, it is clearly seen that the 37-cell tissue tend to experience some extensive shrinkage, which is due to the large number of intercellular interactions present in such models that cause cells to be more packed and attached to each other.

Fig. 9 19-cell tissue model: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Further, when referring to Fig. 10(f) and 12(f), the cell wall has experienced some local wrinkling or warping. This clearly resembles the wrinkled cell walls of realistic dried cells in tissues (see Fig. 5(b) and (c)). Also, when referring to Fig. 9 and 11, the interior and exterior cells of the tissues have undergone different and irregular deformations, which should be due to the differences of intercellular contacts and force fields that they experience. To elaborate these trends quantitatively, the geometric parameters introduced in Section 2.5 were studied and are presented in Fig. 13. It should be noted here that only the internal cells of the tissues were used for the analysis in order to calculate average cellular geometrical parameters, and therefore the outer most cell layer was omitted in each case. It is because, the internal cells are the ones that fully experience the intercellular interactions, which would fairly represent actual cells in a realistic tissue. As seen from Fig. 13(a), (b) and (d), the cell area, ferret diameter and roundness trends further confirm that the bigger tissues experience higher level of shrinkage. Further, the model predictions are fairly in a good agreement with the experimental findings. Particularly the 19-cell tissue and 37-cell tissue shows a slightly rapid shrinkage during the latter part of drying (X/X₀ < 0.6) which replicates what is observed from the experimental curves. When the cell perimeter is considered, all models seem to experience almost identical cell perimeter reduction trends. This is because, the cell perimeter variation is less influenced by the tissue scale interactions, and as defined in eqn (8), the perimeter reduction was taken to be directly governed by the moisture content of the cell. The selected parameters used in the cell wall contraction force formulations have ensured a good agreement between simulations and experimental findings. The cell elongation trends reconfirm that the 19-cell and 37-cell tissues can better mimic the real tissue drying behaviours rather than the single cell model and 7-cell tissue model. Further, as observed from experimental curves, the elongation increases during drying, which implies that the cells are subjected to non uniform deformations along major and minor axes. This is clearly observed from the experiments (Fig. 5) and from the simulated tissue models (Fig. 10 and 12). Next, in the case of cell compactness, bigger tissues generally indicate a fairly good agreement with the gradually reducing compactness trend of experimental curves. All these observations highlight the capability of the proposed tissue model to replicate realistic tissue deformations during drying.

Fig. 10 19-cell tissue model (enlarged view): (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 11 37-cell tissue model: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 12 37-cell tissue model (enlarged view): (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 13 Influence of moisture content reduction for cellular geometric property variations of dense tissues: (a) A/A₀. (b) D/D₀. (c) P/P₀. (d) R/R₀. (e) EL/EL₀. (f) C/C₀. (Error bars indicate one standard deviation.)

Further, we should mention here that in actual tissue drying, the cells that are located in the periphery of the tissue undergo intense moisture removal compared to the interior cells, leading to a case hardening phenomenon. This is particularly evident in tissues with finite size that involve a large number of cells, which undergo a rapid drying process such as in the case of forced convective drying. However, in our simulations, we simply assumed all the cells in the tissue undergo similar moisture removal, since our main objective was to demonstrate how the cells in a tissue would fundamentally deform when the moisture and turgor pressure varies. Although such a simplified approach was used in this work, the general tissue model introduced here is fully capable of simulating such case hardening effects and can easily be achieved by separately time evolving the interior and exterior cells with customized moisture content and turgor pressure values.

4.4 Simulation of deformations during drying: tissues with cellular voids (19, 37 cells)

In addition to the dense tissues simulations, another set of simulations were conducted in order to predict how the tissues with voids would deform during drying. For this purpose, only the 19-cell and 37-cell tissues were used and the voids were created by removing one cell from each tissue. The qualitative results are presented in Fig. 14–17. From Fig. 14 and 16, it is evident that the void geometry also experiences shrinkage in a fairly similar manner to that of the cells, which fundamentally agrees with experimental findings.^12,21,23 Also it is observed that the cells that are located in the periphery of the voids, tend to deform in a different manner compared to the normal cells within the dense regions of the tissues. This is due to the difference of the intercellular interactions, where cells in the periphery of the void have no intercellular contacts towards the void. Further, the enlarged views shown in Fig. 15 and 17 indicate that the cells in the 36-cell tissue tend to undergo limited irregular deformations compared to the 18-cell tissue. It should be due to the higher number of intercellular interactions present in the larger tissue which resist irregular local deformations. In addition to these qualitative analysis, the geometric parameter trends were also studies as presented in Fig. 18. The cell area and roundness trends indicate this particular deformation difference up to some extent. However, due to the averaging of the parameters, the localized changes observed in tissue images are not clearly reflected from cell elongation and compactness trends. These findings imply that the overall deformation characteristics of tissues are influenced by local voids, but in bigger tissues the effect is not so significant both qualitatively and quantitatively, provided that the void percentage is not very high (in this case it is 2–5%). Also, as described in Section 4.3, the different boundary conditions can be used in case of the voids also, by having customized moisture content and turgor pressures for the cells in the periphery of the tissue and around the void. However, considering the scope of this work, these potentially interesting developments were not further studied, but will be addressed in our future works.

Fig. 14 18-cell tissue model with voids: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 15 18-cell tissue model with voids (enlarged view): (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 16 36-cell tissue model with voids: (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 17 36-cell tissue model with voids (enlarged view): (a) initial condition before simulations. (b) Turgid condition: X/X₀ = 1.0 & P_T = 200 kPa. Dried conditions: (c) X/X₀ = 0.8 & P_T = 160 kPa. (d) X/X₀ = 0.6 & P_T = 120 kPa. (e) X/X₀ = 0.4 & P_T = 80 kPa. (f) X/X₀ = 0.3 & P_T = 60 kPa.

Fig. 18 Influence of moisture content reduction for cellular geometric property variations of tissues with voids: (a) A/A₀. (b) D/D₀. (c) P/P₀. (d) R/R₀. (e) EL/EL₀. (f) C/C₀. (Error bars indicate one standard deviation.)

4.5 Sensitivity analysis

The tissue model involves a number of physical and non-physical parameters. A sensitivity analysis was conducted for three parameters that mainly influence the model behaviour. For easy interpretation of results, only a 19 cell tissue model was used for this analysis. Also, in order to better illustrate the model sensitivity to each parameters of interest, both the fresh cell condition (i.e. X/X₀ = 1.0) and the critically dried condition (i.e. X/X₀ = 0.3) are presented as tissue images. Additionally, quantitative data on cell area, perimeter and roundness are presented for all states.

4.5.1 Model sensitivity to cell wall Young's modulus (E). Depending on the biological variability and maturity level of materials, the cell wall Young's modulus (E) can vary. Since the E mainly governs how the cell wall responds to different force interactions, it can highly influence the model behaviour. In order to study the model sensitivity for the cell wall E, a series of simulations was conducted by varying E from 1E to 2E in steps of 0.25E. From Fig. 19–22, it is observed that the influence of E is dominant in dried tissues than fresh tissues. Further, when comparing 1E and 2E cases from Fig. 21 and 22, it is clearly observed that the dried tissue deformations are resisted by stiffer cell walls and the resulting dried tissue are comparatively larger. Additionally, the 1E and 1.25E tissues have undergone considerable local wrinkling compared to 2E tissue, which also supports the above mentioned behaviour of stiffer cell walls. This is further evident from Fig. 23 where cell area and roundness of 1.75E and 2E models in particular are subjected to limited shrinkage compared to other less stiff tissue models.

Fig. 19 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for cell wall Young's modulus (E): (a) 1E (=54 MPa). (b) 1.25E. (c) 1.5E. (d) 1.75E. (e) 2E.

Fig. 20 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for cell wall Young's modulus (E) (enlarged view): (a) 1E (=54 MPa). (b) 1.25E. (c) 1.5E. (d) 1.75E. (e) 2E.

Fig. 21 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for cell wall Young's modulus (E): (a) 1E (=54 MPa). (b) 1.25E. (c) 1.5E. (d) 1.75E. (e) 2E.

Fig. 22 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for cell wall Young's modulus (E) (enlarged view): (a) 1E (=54 MPa). (b) 1.25E. (c) 1.5E. (d) 1.75E. (e) 2E.

Fig. 23 Effects of cell wall Young's modulus (E) on normalized cellular parameters of a 19-cell tissue: (a) A/A₀. (b) P/P₀. (c) R/R₀. (Error bars indicate one standard deviation.)

4.5.2 Model sensitivity to pectin layer stiffness (k_pectin). Pectin layer stiffness mainly governs the extent to which the cells are bonded to each other. It is clearly evident from Fig. 24 and 25, that the fresh cells are tightly bonded to each other in case of higher k_pectin. When comparing Fig. 25(a) and (e), it is observed that the cell shape also tend to be fairly hexagonal when higher k_pectin values are used. Further, when referring to Fig. 26 and 27, it is observed that higher k_pectin values tend to have the cells in dried tissues attached to each other with fairly consistent gaps between them. This is clearly observed from Fig. 26(e) and 27(e), where the spacing around the centre-most cell is almost equal to the initially set pectin layer gap. Fig. 26 and 27 further indicate that k_pectin particularly influences localized changes of cell shapes in dried tissue rather than global tissue shape changes (Fig. 26(a) and (e) has almost same tissue geometry even the used k_pectin values are significantly different). The qualitative results are presented in Fig. 28 and only the cell area and the roundness indicate some level of a difference between different k_pectin models. This should be mainly due to the cancellation of the effects during the averaging and normalisation process.

Fig. 24 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for pectin layer stiffness (k_pectin): (a) 0.5 k_pectin. (b) 1 k_pectin (=20 N m⁻¹). (c) 1.5 k_pectin. (d) 2.0 k_pectin. (e) 2.5 k_pectin.

Fig. 25 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for pectin layer stiffness (k_pectin) (enlarged view): (a) 0.5 k_pectin. (b) 1 k_pectin (=20 N m⁻¹). (c) 1.5 k_pectin. (d) 2.0 k_pectin. (e) 2.5 k_pectin.

Fig. 26 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for pectin layer stiffness (k_pectin): (a) 0.5 k_pectin. (b) 1 k_pectin (=20 N m⁻¹). (c) 1.5 k_pectin. (d) 2.0 k_pectin. (e) 2.5 k_pectin.

Fig. 27 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for pectin layer stiffness (k_pectin) (enlarged view): (a) 0.5 k_pectin. (b) 1 k_pectin (=20 N m⁻¹). (c) 1.5 k_pectin. (d) 2.0 k_pectin. (e) 2.5 k_pectin.

Fig. 28 Effects of pectin layer stiffness (k_pectin) on normalized cellular parameters of a 19-cell tissue: (a) A/A₀. (b) P/P₀. (c) R/R₀. (Error bars indicate one standard deviation.)

4.5.3 Model sensitivity to cell wall contraction force coefficient (k_wc). As described in Section 5.1, the cell wall contractions are directly influenced by the k_wc and four different k_wc values were tested in this study against the original model. By referring to the eqn (8), it is clearly evident that the model is not influenced by k_wc when X/X₀ = 1.0. This is reconfirmed by Fig. 29 and 30, where fresh tissues are having identical geometries even at different k_wc values. However, in case of dried tissues, as seen from Fig. 31 and 32, a clear difference is observed. Fig. 31 implies that higher k_wc values tend to increase tissue shrinkage that can even lead to bulk level tissue shape changes. This is particularly evident when comparing Fig. 31(a) and (e), where the former tissue geometry is fairly hexagonal and the latter one is quite circular. Further, when referring to Fig. 32, we can argue that the cell-wise shrinkage differences should be the main cause for these bulk level shrinkage differences of the tissue. Also in Fig. 32, cells which have higher k_wc have undergone significant local wrinkling and warping, which resemble actual cell deformations during drying. These differences are further evident from the geometric parameters as shown in Fig. 33, where the cell area, perimeter and the roundness trends indicate that, higher k_wc models tend to experience lower shrinkage. Also, particularly in cell perimeter trends as shown in Fig. 33(b), the k_wc = 0 model only subjects to very limited cell perimeter reductions and deviates significantly from experimental curves, indicating a general decreasing trend. This highlights the importance of using k_wc in the drying model, which should be carefully tuned in order to ensure better model performance.

Fig. 29 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for cell wall contraction force coefficient (k_wc): (a) 0 .k_wc (b) 0.5 .k_wc (c) 1 k_wc (=4 × 10 ⁴ N m⁻¹). (d) 1.5 k_wc. (e) 2 k_wc.

Fig. 30 Sensitivity of 19-cell tissue at X/X₀ = 1.0 for cell wall contraction force coefficient (k_wc) (enlarged view): (a) 0 k_wc. (b) 0.5 k_wc. (c) 1 k_wc (=4 × 10 ⁴ N m⁻¹). (d) 1.5 k_wc. (e) 2 k_wc.

Fig. 31 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for cell wall contraction force coefficient (k_wc): (a) 0 k_wc. (b) 0.5 k_wc. (c) 1 k_wc (=4 × 10 ⁴ N m⁻¹). (d) 1.5 k_wc. (e) 2 k_wc.

Fig. 32 Sensitivity of 19-cell tissue at X/X₀ = 0.3 for cell wall contraction force coefficient (k_wc) (enlarged view): (a) 0 k_wc. (b) 0.5 k_wc. (c) 1 k_wc (=4 × 10 ⁴ N m⁻¹). (d) 1.5 k_wc. (e) 2 k_wc.

Fig. 33 Effects of the cell wall contraction force coefficient (k_wc) on normalized cellular parameters of a 19-cell tissue: (a) A/A₀. (b) P/P₀. (c) R/R₀. (Error bars indicate one standard deviation.)

5. Conclusions and outlook

A 2-D meshfree particle based model has been developed by coupling SPH with DEM in order to simulate microscale morphological changes of plant tissues during drying. We specifically used a hexagonal cell shape and a positive pectin layer gap in the tissue model in order to represent the frequently observed honeycomb cellular structure of plant tissues, which also owns middle lamella of a finite thickness. The main findings of this work can be summarized as:

• Using a meshfree approach, this work has clearly overcome the key shortcomings of grid-based plant tissue drying models¹⁷ that have limited capability to account for excessive moisture reduction during drying as well as cell wall wrinkling. Also, the simulation approach in this work is comparatively straight forward.

• Compared to the state of the art SPH–DEM based plant tissue models,^23,24 the hexagonal initial cell arrangement and the positive gap used for the middle lamella in this work resulted in a more realistic tissue structure with improved cell–cell contacts.

• In order to ensure the stability of this novel model, a careful selection of cell wall pectin layer and LJ repulsion force parameters are recommended as is the proposed smooth time evolution approach which initially omits the cell wall bending stiffness effects and the cell wall contraction force effects.

• In order to avoid cell–cell penetrations in the tissue model, only the cell wall particles can be used rather than the large array of fluid particles which was used by the state of the art SPH–DEM based tissue models. By this means, the computational cost can be significantly reduced, as a result of lower number of particle interaction calculations.

• We introduced a means of simulating dried tissues using a SPH–DEM particle scheme by varying the cell moisture content and turgor pressure. A series of simulations were conducted on different size tissues (7-cell, 19-cell and 37-cell) using a generic dense arrangement and it was observed that the tissue model undergoes cellular and bulk level shrinkage. Particularly in critically dried states, cell wall wrinkling and warping was clearly observed. Both effects were found to be in close agreement with experimental findings of dried plant tissues.

• Qualitative studies on several key cellular geometric parameters indicated that the bigger tissue models are fairly capable of replicating most of the realistic cellular geometrical changes during drying.

• Tissue models with voids resembled a basic porous cellular structure which can be observed in real tissues and it was observed that the voids undergo similar deformation as of actual cells, which is in good agreement with the experimental findings in the literature.

• The model is highly sensitive to the cell wall Young's modulus particularly in dried conditions than fresh conditions. Stiffer walled cells tend to resist shrinkage and resulted in larger dried tissue dimensions.

• Stiffer pectin layers tend to produce strong cell–cell contacts which cause fresh cell shapes in tissues to become more hexagonal. Also, higher pectin layer stiffness tends to maintain fairly uniform gaps between cells, particularly in dried tissues.

• The model is also highly sensitive to the force coefficient of cell wall contractions and it mainly influences the dried tissue deformations by enhancing cellular shrinkage and cell wall warping.

Considering all these, we conclude that the proposed SPH–DEM tissue model provides a reliable means of numerically studying the cellular drying phenomena which involves complex physics. The model can be further improved by using advanced cell wall models, fluid–wall interaction models and incorporating multiscale techniques to model bulk tissues in a computationally efficient manner. Tissue porosity effects, bulk material deformations and shape changes can be further studied and related with product and process parameters in the drying industry. Also, 3-D models can be developed using the fundamental conceptual framework proposed in this work. This model can be regarded as a generic model valid for most of the plant materials and could be used to predict deformation characteristics during drying upon the availability of relevant model parameters. Also, temperature variations can be fundamentally incorporated into the SPH–DEM formulations in order to account for temperature influence for dried plant tissue deformations.

Acknowledgements

This research was conducted at the experimental facilities of the Queensland University of Technology (QUT) – Brisbane, Australia and the authors extend their thanks to the graduate student Ms Parva Hesami for the contributions in experiments. The research studies were financially supported by QUT, International Postgraduate Research Scholarship (IPRS) and Australian Postgraduate Award (APA) scholarship. Support from the ARC Future Fellowship Grant (FT100100172) and the High Performance Computer resources provided by the QUT are also gratefully acknowledged as is the support by University of Ruhuna – Sri Lanka.

References

1. S. Grabowski, M. Marcotte and H. S. Ramaswamy, in Handbook of postharvest technology: cereals, fruits, vegetables, tea, and spices, ed. A. Chakraverty, A. S. Mujumdar, G. S. V. Raghavan and H. Rawaswamy, Marcel Dekker, New York 2003, ch. 23, pp. 653–695.
2. S. V. Jangam, Drying Technol., 2011, 29, 1343–1357.
3. B. P. Hills and B. Remigereau, Int. J. Food Sci. Technol., 1997, 32, 51–61.
4. P. P. Lewicki and J. Drzewucka, Effect of drying on tissue structure of selected fruits and vegetables, Greece, 1998.
5. I. N. Ramos, C. L. M. Silva, A. M. Sereno and J. M. Aguilera, J. Food Eng., 2004, 62, 159–164.
6. L. Mayor, M. A. Silva and A. M. Sereno, Drying Technol., 2005, 23, 2261–2276.
7. P. P. Lewicki and G. Pawlak, Drying Technol., 2003, 21, 657–683.
8. Y. Bai, M. S. Rahman, C. O. Perera, B. Smith and L. D. Melton, J. Agric. Food Chem., 2002, 50, 3179–3185.
9. M. S. Rahman, I. Al-Zakwani and N. Guizani, J. Sci. Food Agric., 2005, 85, 979–989.
10. T. Funebo, L. l. a. Ahrné, S. Kidman, M. Langton and C. Skjöldebrand, J. Food Eng., 2000, 46, 173–182.
11. M. K. Bartlett, C. Scoffoni and L. Sack, Ecol. Lett., 2012, 15, 393–405.
12. G. H. Crapiste, S. Whitaker and E. Rotstein, Chem. Eng. Sci., 1988, 43, 2919–2928.
13. H. X. Zhu and J. R. Melrose, J. Theor. Biol., 2003, 221, 89–101.
14. C. X. Wang, L. Wang and C. R. Thomas, Ann. Bot., 2004, 93, 443–453.
15. N. Wu and M. J. Pitts, Postharvest Biol. Technol., 1999, 16, 1–8.
16. Q. Gao and R. E. Pitt, Trans. ASABE, 1991, 34, 0232–0238.
17. S. W. Fanta, M. K. Abera, W. A. Aregawi, Q. T. Ho, P. Verboven, J. Carmeliet and B. M. Nicolai, J. Food Eng., 2014, 124, 86–96.
18. T. Rudge and J. Haseloff, in Advances in Artificial Life, ed. M. Capcarrère, A. Freitas, P. Bentley, C. Johnson and J. Timmis, Springer, Berlin, Heidelberg, 2005, vol. 3630, ch. 9, pp. 78–87.
19. H. Honda, M. Tanemura and T. Nagai, J. Theor. Biol., 2004, 226, 439–453.
20. Z. Liu, W. Hong, Z. Suo, S. Swaddiwudhipong and Y. Zhang, Comput. Mater. Sci., 2010, 49, S60–S64.
21. G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific Publishing Co., Singapore, 2003.
22. R. A. Gingold and J. J. Monaghan, Mon. Not. R. Astron. Soc., 1977, 181, 375–389.
23. P. Van Liedekerke, P. Ghysels, E. Tijskens, G. Samaey, D. Roose and H. Ramon, Soft Matter, 2011, 7, 3580–3591.
24. P. V. Liedekerke, P. Ghysels, E. Tijskens, G. Samaey, B. Smeedts, D. Roose and H. Ramon, Phys. Biol., 2010, 7, 026006.
25. P. Van Liedekerke, E. Tijskens, H. Ramon, P. Ghysels, G. Samaey and D. Roose, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 061906.
26. H. C. P. Karunasena, W. Senadeera, Y. T. Gu and R. J. Brown, A Coupled SPH–DEM Model for Fluid and Solid Mechanics of Apple Parenchyma Cells During Drying, Launceston, Australia, 2012.
27. H. C. P. Karunasena, W. Senadeera, Y. T. Gu and R. J. Brown, A particle based micromechanics model to simulate drying behaviors of vegetable cells, Gold Coast, Australia, 2012.
28. H. C. P. Karunasena, W. Senadeera, Y. T. Gu and R. J. Brown, Appl. Math. Model., 2014 DOI:10.1016/j.apm.2013.12.004.
29. A. E. Smith, K. E. Moxham and A. P. J. Middelberg, Chem. Eng. Sci., 1998, 53, 3913–3922.
30. M. C. Jarvis, Plant, Cell Environ., 1998, 21, 1307–1310.
31. M. A. J. Chaplain, J. Theor. Biol., 1993, 163, 77–97.
32. L. Taiz and E. Zeiger, in Plant Physiol., Sinauer Associates, Sunderland, USA, 2010, ch. 3, pp. 73–84.
33. H. C. P. Karunasena, W. Senadeera, R. J. Brown and Y. T. Gu, Eng. Anal. Bound. Elem., 2014 DOI:10.1016/j.enganabound.2014.04.004.
34. S. M. Hosseini and J. J. Feng, Chem. Eng. Sci., 2009, 64, 4488–4497.
35. T. W. Pan and T. Wang, Int. J. Numer. Anal. Model., 2009, 6, 455–473.
36. L. Shi, T.-W. Pan and R. Glowinski, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 016307.
37. H. C. P. Karunasena, P. Hesami, W. Senadeera, Y. T. Gu, R. J. Brown and A. Oloyede, Drying Technol., 2014, 32, 455–468.
38. A. Blum, in Plant Breeding for Water-Limited Environments, Springer, New York, 2011, ch. 2, p. 26.
39. P. V. Liedekerke, P. Ghysels, E. Tijskens, G. Samaey, B. Smeedts, D. Roose and H. Ramon, presented in part at the 4th international SPHERIC workshop, Nantes, France, 26–29 May, 2009.
40. J. P. Morris, P. J. Fox and Y. Zhu, J. Comput. Phys., 1997, 136, 214–226.
41. A. Colagrossi, B. Bouscasse, M. Antuono and S. Marrone, Comput. Phys. Commun., 2012, 183, 1641–1653.
42. J. G. Marshall and E. B. Dumbroff, Plant Physiol., 1999, 119, 313–319.
43. P. M. Neumann, Crop Sci., 1995, 35, 1258–1266.
44. D. J. Barker, C. Y. Sullivan and L. E. Moser, Agron. J., 1993, 85, 270–275.
45. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
46. J. C. Araujo, F. C. Téran, R. A. Oliveira, E. A. A. Nour, M. A. P. Montenegro, J. R. Campos and R. F. Vazoller, J. Electron Microsc., 2003, 52, 429–433.
47. D. F. Bray, J. Bagu and P. Koegler, Microsc. Res. Tech., 1993, 26, 489–495.

---

Footnotes

†

‡ 4πA/P².

§

¶ Major axis length/minor axis length.

---