Application of a coupled smoothed particle hydrodynamics (SPH) and coarse-grained (CG) numerical modelling approach to study three-dimensional (3-D) deformations of single cells of diﬀerent food-plant materials during drying

Numerical modelling has gained popularity in many science and engineering streams due to the economic feasibility and advanced analytical features compared to conventional experimental and theoretical models. Food drying is one of the areas where numerical modelling is increasingly applied to improve drying process performance and product quality. This investigation applies a three dimensional (3-D) Smoothed Particle Hydrodynamics (SPH) and Coarse-Grained (CG) numerical approach to predict the morphological changes of diﬀerent categories of food-plant cells such as apple, grape, potato and carrot during drying. To validate the model predictions, experimental findings from in-house experimental procedures (for apple) and sources of literature (for grape, potato and carrot) have been utilised. The subsequent comaprison indicate that the model predictions demonstrate a reasonable agreement with the experimental findings, both qualitatively and quantitatively. In this numerical model, a higher computational accuracy has been maintained by limiting the consistency error below 1% for all four cell types. The proposed meshfree-based approach is well-equipped to predict the morphological changes of plant cellular structure over a wide range of moisture contents (10% to 100% dry basis). Compared to the previous 2-D meshfree-based models developed for plant cell drying, the proposed model can draw more useful insights on the morphological behaviour due to the 3-D nature of the model. In addition, the proposed computational modelling approach has a high potential to be used as a comprehensive tool in many other tissue morphology related investigations.


Introduction
Drying is one of the most common and cost effective techniques for the preservation of food and for the production of traditional as well as innovative processed products. 1 According to statistical information, drying is generally employed to preserve approximately 20% of the entire world's perishable crops annually. 2 During drying, moisture content of the food cellular structure is reduced to increase the shelf life by obstructing the microbiological activities. The moisture content reductions lead to volumetric deformations in the cellular structure. These changes directly influence the drying process performance and dried food quality. Therefore, a sound understanding on the underlying mechanisms is necessary to optimally control such characteristics as influenced by microstructural deformations. In addition, moisture content [3][4][5][6][7][8] and drying temperature 9 also act as driving forces for the microstructural deformations during drying. Moisture content has a strong correlation with cell turgor pressure, 10 and drying temperature is directly correlated with the rate of the moisture removal from the drying environment. To analyse such relationships, numerous microscale theoretical 11,12 and empirical models 3,5,13,14 have been reported in literature.
Most of the food drying modelling applications in the literature study the relationship between macro-scale characteristics and drying process parameters. 15 Usually, these models are based on empirical or theoretical investigations. In general, they lack common applicability due to the specific experimental setup usage. 14 Accordingly, the established relationships tend to have limitations in the in validity for the specific plant food categories under pre-specified drying conditions. 16 The empirical coefficients used in those models could deviate from realistic behaviours. 17 The theoretical models are generally based on heat and mass transfer principles (e.g. Fick's law). As they use sophisticated theoretical relationships along with oversimplified boundary conditions, there are limitations in the applicability. 15 In addition, as nonlinear diffusion equations are often used in these theoretical drying models, the computational cost of the numerical implementation have the tendency to be excessive. 14,18,19 Numerical modelling is an efficient and effective tool in morphological analysis of various types of materials. Until the recent times, it had not been used for comprehensively studying morphological changes of food-plant cellular structure during drying. However recently, a meshfree-based novel numerical modelling approach has attracted attention as a feasible technique for serving this purpose. [20][21][22] Through such comprehensive numerical models, benefits could be achieved in food drying engineering to improve the drying process performance and product quality. 14 Further, there has been progress in the field of computational modelling of morphological instability and related surface wrinkling phenomena in soft materials. 23 Surface wrinkling phenomena present in dried plant cellular structure is an interesting area that has not been explored with a computational modelling perspective until the recent times. Microstructural morphological changes in plant cellular structure during drying often involves large deformations, nonlinear constitutive relations and complicated multiphase phenomena. This makes it difficult to be modelled and simulated with conventional grid-based methods such as Finite Element Method (FEM) and Finite Volume Method (FVM). Hence, numerically modelling morphological characteristics of plant cellular structure during drying using a meshfree approach has gained popularity due to the ability to handle such complicated physics in a versatile manner. 14 Accordingly, there are a number of recently reported efforts to numerically model the macro and micro mechanics of food-plant cells/tissues using a coupled Smoothed Particle Hydrodynamics (SPH) and Discrete Element Method (DEM) approach. 14 Some of these studies focus only on the fresh plant cellular structure and its behaviour under external mechanical loading. 24 There are some studies considering the morphological changes of both fresh and dried plant cellular structure in two dimensions (2-D). 22 The predictions of this model agree well with the experimental findings both quantitatively and qualitatively. 25 Drying phenomena have been incorporated into the models through changes of the moisture content, turgor pressure and cell wall contractions. The higher capability of this meshfree-based approach is evident, particularly in handling large deformations and a higher degree of moisture reduction. 22 The computational approach in these investigations has been numerically validated through computational consistency comparisons. 20 In addition, the predictions of these 2-D models have been validated by comparing the model predictions with experimental findings. 26 However, there is a major limitation in these plant cellular drying models. It is the 2-D nature inhibiting the opportunity to study the true cellular level deformations which are 3-D by nature. This highlights the necessity for a 3-D meshfree-based cellular drying model, which is a key aim of our investigation.
To develop such a 3-D cellular drying model, there are several conceptual constraints in using the Discrete Element Method (DEM) to represent the cell wall membrane. Originally DEM had been introduced to the numerical modelling field in order to solve problems in soil mechanics. According to DEM fundamentals, particles can have different geometries and physical properties. The interactions between DEM particles are indirect ones. On the other hand, plant cell wall can be treated as a continuous thin membrane consisting of various biopolymers. The literature [15][16][17] suggests that conceptually, a Coarse-Grained method (CG) could be more suitable for this application. In CG, a network of representative particles or molecules is used to define the entire system where the characteristics of the entire system are concentrated into these particles. 27 The interactions between CG particles can be specifically interpreted to reflect the real physical behaviour of the systems. The literature suggests that CG has been widely used to model and simulate various kinds of biopolymers including the ones that make up the plant cell walls. [28][29][30] Within this background, the aim of this investigation was to develop a more realistic three-dimensional (3-D) numerical model using a SPH-CG coupled approach. Further, it aimed to simulate the morphological behaviour of apple, potato, grape and carrot plant cells during drying. In this article, firstly the numerical modelling methodology will be discussed followed by the computational implementation methodology. Next, 3-D model predictions will be qualitatively and quantitatively analysed through comprehensive comparisons with relevant experimental findings. Finally, the key conclusions will be discussed along with potential for future work.

3-D particle representation of a single cell
Plant cells are the fundamental building blocks which make up the whole plant structure. There are different physiological and biochemical functions implemented by different types of cells and the agglomeration of all these functions establishes a biological unity. Out of the main kinds of the tissues in plants, ground tissue system contributes to the largest part and the parenchyma cells are the major component involved there. 31 Parenchyma refers to the tissue consisting of living cells. Usually parenchyma cells have thin cell walls and they usually occur as a continuous medium in the primary plant body. 31,32 They also play crucial roles of regeneration and wound healing. In this background, the main focus of this study was on parenchyma cells. The thinness of the parenchyma cell walls was also a reason behind the selection as it significantly reduces the complexity of the cell wall. Many state-of-the-art theoretical and numerical studies have focused mainly on the single cell systems, which assist understanding the fundamental characteristics of the cellular dynamics during drying. 20,25 Accordingly, in recent meshfree-based single cell modelling approaches introduced by Van Liedekerke et al. 24,33,34 and Karunasena et al., [20][21][22]25,26,[35][36][37] the cell fluid and cell wall have been considered as the two main components of a single cell. 14 There, the cell fluid has been approximated to a viscous homogeneous Newtonian liquid and modelled with Smoothed Particle Hydrodynamics (SPH). In this study, we also used the same approach to treat the cell fluid. The cell wall was considered as a semi-permeable solid membrane with viscoelastic properties, and modelled as an incompressible Neo-Hookean material using a Coarse-Grained approach (CG). At the same time, the morphological descriptions of plant cells which come from relevant experimental literature were taken into consideration 8,38,39 in the design and development of this 3-D model. Accordingly, the basic structure of a plant parenchyma cell was considered as a flexible solid wall enclosing a fluid mass. 8,20,31,38,[40][41][42][43][44][45] A physical balance was hypothesised between the cell fluid hydrodynamic pressure and the cell wall tension. In other words, the flexible cell wall holds the cell fluid mass inside by balancing the forces exerted by the fluid and in return, there is a stress development in the cell wall itself. Depending on the nature of the dynamics that the cell is subjected, the cell wall can either stretch or contract (i.e. inflation or shrinkage). The basic geometry of a single cell was considered to be spherical for all the food-plant categories. 38 The cell fluid geometry was approximated to a solid sphere and the cell wall was taken as a hollow threedimensional (3-D) spherical shell enclosing the fluid content ( Fig. 1). The cell fluid was assumed to be incompressible and the whole system is regarded as isothermal. 20,25 With these fundamental approximations, the physical nature of the initial cell model was established. Next, the cell fluid and cell wall were discretised using particle schemes in order to computationally implement the model (see Fig. 1). This is a basic principle that comes under SPH, [46][47][48] where the interactions among the particles are described using a number of governing equations representing different types of interacting force fields. Due to the flexibility of the SPH-CG particle framework used in this model, there is a significant potential to be upgraded to a tissue system. 24,34 When considering the similar work in this line of research, some researchers have only addressed the 2-D behaviour of plant cellular systems, [20][21][22]25,26,35,36 whereas some others have studied the mechanical response of fresh cells and tissues under external mechanical loading. 24,33,34 Therefore, this investigation addresses the research gap of modelling three dimensional (3-D) morphological characteristics of single plant cells during drying.

Cell fluid model
The water content of the cell protoplasm (fluid) can be as high as 80-90% by volume of the entire cell. 34 Therefore, for numerical modelling purposes, the cell fluid can be approximated to an incompressible homogeneous Newtonian fluid equivalent to water. 14 Accordingly, the Navier-Stokes equations can be used for modelling. However, an elevated viscosity value has to be incorporated in order to comply with low Reynolds number viscous characteristics. 20 Accordingly, as depicted in Fig. 2, the cell fluid can be modelled with four different types of force interactions: pressure forces (F p ), viscous forces (F v ), wall-fluid repulsion forces (F rw ) and wall-fluid attraction forces (F a ). 49 The summation of all four forces defines the total force F i on any fluid particle i as (eqn (1)), where i 0 represents the neighbouring fluid particles and k, the interacting wall particles. According to the standard Lagrangian type SPH equations, which are used to model weakly compressible low Reynold's number fluid flows, 24,50 the momentum equation estimates the pressure forces (F p ii 0 ) and viscous forces (F v ii 0 ) for any given fluid particle i as a summed influence from its neighbouring fluid particles i 0 . These are formulated as (eqn (2) and (3)), Here, m, P, r, m, v and W are the fluid particle mass, cell turgor pressure, density, dynamic viscosity, velocity and the smoothing kernel. The cubic spline has been selected as the smoothing kernel as given in eqn (4) which is a widely used kernel function in the recent SPH studies due to its stability as well as computational efficiency. 48 Here, h is the smoothing length at the current time step, s is the ratio of r ii 0 /h and r ii 0 is the distance between particle i and any neighbouring fluid particle i 0 within the influence domain of the particle i (0 r s r 2). From eqn (2)-(4), it is evident that the smoothing kernel (W) and the value of s is influential in determining pressure and viscous forces (F p ii 0 and F v ii 0 ). It is because the values of field properties are taken as smoothed and summed influences across the SPH influence domain of each particle. 48 With the subsequent deformations of the cellular model, the smoothing length h has to evolve. This is a computational requirement for maintaining an optimum number of fluid particles within the influence domain of each particle. 25 For this purpose, a simple geometrical relationship has been employed as, Here, D is the average cell diameter at the current time step, D 0 is the initial cell diameter and h 0 is the initial smoothing length. With the numerical evolution of the system with time, an Equation of State (EOS) (eqn (6)) is used to maintain the relationship between the density and the pressure.
Here, P T is the initial cell turgor pressure, K the fluid compression modulus, r i the density of each fluid particle at the current time step, and r 0 the initial density of the cell fluid. The value of fluid compression modulus (K) should be so that the fluid behaves in a sufficiently incompressible manner. 24 Accordingly, the density r of each particle is defined as below Here v is the volume of the fluid subdomain represented by a given particle. When this equation is differentiated with respect to time (eqn (8)), To update the density, the standard SPH continuity equation is used as given in eqn (9), 48 Here, r i * denotes the density of a given particle assuming a constant fluid particle mass. 34 The first term in the right hand side of eqn (8) is the change of density due to deformations of the cell, and the second term corresponds to the change in water content of the cell. Since plant cells have semipermeable walls, as long as the turgor pressure of the cell does not equal to the osmotic potential P (P o 0), there will be a net water transport across the cell wall. 31 If the cell fluid mass loss or gain is not significant, P could be assumed constant. This fluid mass transport through the cell wall results in changes of the mass of individual fluid particles according to the constitutive relationship given in eqn (10). 31 Here, L P is the cell wall permeability (hydraulic conductivity) which is assumed to be isotropic over the cell surface, n f the number of fluid particles and A C the total cell surface area.
If the cell absorbs water, the density will initially increase augmenting the pressure. This is counterbalanced by the push from the cell fluid in the outward direction, lowering the density. The final density, which should vary only slightly from the initial density, will be obtained when the fluid particles cease to move further, i.e. when there is a physical balance between the fluid pressure and the cell wall tension. 34 The repulsion forces F rw ik and attraction forces F a ik act on a given fluid particle i due to the influence of the surrounding wall particles. These forces are described similar to their counterpart LJ forces in the cell wall domain and are defined as given in eqn (11) and (12) (see Section 2.3 for more details).

Cell wall model
Wall is the system boundary of a cell. It plays a critical role in defining the morphology of a plant. Plant cell walls are mainly made up of biopolymers including cellulose, hemicellulose and pectin. [28][29][30] The combination of these biopolymers provides the mechanical strength to a cell and eventually to the whole plant.
The bio-polymeric cell wall membrane could not be analysed with a simple linear elasticity theory. Mechanically, the cell wall material exhibits both elastic and plastic behaviour while energy dissipation could be attributed to viscous and structural damping. This behaviour strongly depends on the time scale of the analysis. 34,51 In this particular work, a Coarse-Grained (CG) approach was used where the representative CG particles were distributed on a spherical surface, having a local connectivity while interacting through a number of force fields. Accordingly, each CG element carries properties of the corresponding cell wall segment. The deformations are represented by the displacement of respective particles using six types of force interactions: Stiff forces (harmonic energy) (F e ), damping forces (F d ), wall-fluid repulsion forces (F rf ), wall-fluid attraction forces (F a ), bending stiffness forces (F b ) and cell wall contraction forces (F c ) as presented in Fig. 3. 20,22,25 Accordingly, the total force (F k ) on any wall particle k is derived as in eqn (13), Here, for each wall particle k, i are the neighbouring fluid particles and j bonded wall particles. Stiff forces are simply defined using a spring network model to represent the cell wall resistance to any extensions or contractions using the derivation from the harmonic energy. 49 Accordingly, the stiff force F e kj on any wall particle k due to any other bonded wall particle j is calculated individually for each wall element as given in eqn (14). 34 Here, G is the shear modulus (EE/3), E the Young's modulus of the cell wall, l = l/l 0 the stretch ratio of any cell wall element, l the current length of the wall element (distance between particle k and j), l 0 its initial undeformed length, t 0 the initial cell wall thickness and n the unit normal vector. Damping forces (F d ) were incorporated to the model in order to account for viscous characteristics of the fibrous-polymeric cell wall material and have been defined using a linear dashpot model. Accordingly, the viscous forces F d kj acting on any wall particle k were calculated as given in eqn (15). 20 Here, g is the wall damping constant and v kj the velocity of the particle k relative to particle j. Wall-fluid interactions and boundary conditions were defined using wall-fluid repulsion forces (F rf ) and wall-fluid attraction forces (F a ). Both these force fields were defined as Lennard-Jones (LJ) forces. Wallfluid repulsion forces assure that the fluid particles are constrained within the cell wall. It helps to avoid any unrealistic fluid particle penetrations (slip conditions) through the cell wall which can lead to computational instability. These F rf forces act through the centre of any interacting wall-fluid particle pair of interest, in an equal and opposite manner (Fig. 3). The repulsion force F rf ki on a wall particle k due to another fluid particle i is defined as given in eqn (16). 20,24,48 Here, f rf ki is the magnitude of the force and x ki is the position vector of k relative to i. The f rf ki is defined as given in eqn (17). 24 Here, r 0 is the initial gap between two particles, r ki the current gap between them and f rf 0 the LJ contact strength. wall-fluid attraction forces have been defined with the intention of preventing the fluid particles from unrealistically detaching from the cell wall at different cell dryness states. For this purpose, the attraction force F a ki on a wall particle k due to a neighbouring fluid particle i is defined similar to that of F rf ki using a LJ force type with a LJ contact strength of f a 0 as given in eqn (18). These forces apply only when the distance between fluid and wall particles increases relative to the initially defined value.
Karunasena et al. have used wall bending stiffness forces (F b ) and wall contraction forces (F c ) in the physical interaction description of their 2-D coupled SPH-DEM cellular drying models. 20,25 The bending stiffness forces represent the bending resistance present in plant cell wall microstructures. They resist any unrealistic deformations which could occur in the cell wall during simulations. 52 In alignment with a number of recent studies on red blood cell models, 3-D wall bending forces (F b ) were incorporated through wall bending stiffness and bending energy (E b ). [53][54][55][56] The bending energy (E b ) could be determined using the formulation shown below (eqn (19)).
Here k b is the bending stiffness and y, the external angle between the adjacent wall elements (see Fig. 3), Dy, the difference in the angle compared to the previous time step and L n , a geometry specific parameter. 20 According to the geometry of the cell wall particle network, the equation for E b was rewritten as, where n ijk and n ilj are the unit normal vectors specific to the geometry used in this computational procedure. Next, this equation was partially differentiated to obtain the bending force (F b ) acting on the corresponding particle (eqn (21)) Next, wall contraction forces (F c ) were included to represent the cell perimeter reductions observed in drying experiments. 25 An empirical-analytical formulation which was defined by Karunasena et al. in their SPH-DEM cell modelling studies 25 was used to calculate the 3-D wall contraction force term (F c ) as given in eqn (22).
Here, k wc is the force coefficient for wall contractions, L, the current width of a given wall element, L 0 0 the width at fully turgid condition. a and b are empirical factors which were selected depending on the plant cell category (i.e. apple, potato etc.) 25 and the normalised moisture content of the dry cell to be simulated (X/X 0 ).

Modelling different categories of food-plant materials
Selection of model parameters. The above discussed model formulation was used for modelling different plant cell categories by customising the model through appropriate model parameters (see Table 1). These parameters were obtained through microscopic experiments and other numerical models reported in literature. 35 In addition, some model parameter values were 'set' at their optimum values via trial and error methods. Furthermore, Table 2 summarises all the other model parameters which were commonly used for all cell types.
Setting the particle scheme and the smoothing length. Corresponding to different cell types, different particle schemes were used during the simulations to accommodate different sizes (i.e. cell radius) of the cells. Accordingly, separate particle schemes were used for modelling carrot and potato cells while same particle scheme was used for apple and grape cells as they have similar cell radii. In addition. As a result, the smoothing length (h) of the SPH calculations had to be varied to maximise the computational accuracy of the developed 3-D SPH-CG models. This was implemented through determining the percentage model consistency error for each particle scheme (see Section 3 for more details).

Computational implementation and numerical evolution of the model during simulations
COMSOL Multiphysics software was used to create the initial particle geometries. Here, the fluid particle scheme was placed without any interconnections among particles, adhering to the meshfree fundamentals used in SPH. In the cell wall, the particles were placed adhering to the CG fundamentals, mainly with equal spacing. 20,22,24 During the numerical evolution of the model, the mass of the cell fluid tends to change leading to minor density variations. 20 The key driving force behind this is the difference between the cell turgor pressure and the osmotic potential. For the time integration of the particle scheme, Leapfrog method 22,48 was used where the magnitude of the time step was determined through the Courant-Friedrichs-Lewy (CFL) stability criterion. 48,64 In alignment with SPH basics, the variation of the fluid density leads to fluctuation of the cell turgor pressure as governed by an equation of state (EOS) as given in eqn (6). Such turgor pressure variations drive the cell wall inwards (during shrinkage) or outwards (during inflation), causing variations of the cell volume. Based on such cell volume changes, the cell turgor pressure varies again since it has to be counterbalanced by the cell wall tension. The changes in cell turgor pressure leads to the cell fluid mass increments or losses as governed by a mass transfer equation represented in eqn (10). 20 This cycle repeats during the computational evolution of the cell and eventually reaches the steady state condition.
For the implementation of boundary conditions, a method that has been successfully used in a number of recent meshfreebased numerical modelling studies on plant cells 20,22,24,25,[33][34][35] has been employed. 14 It mainly involves a Lennard-Jones (LJ) type approach, where meshfree particles are used to represent the fluid and wall particles along with repulsive LJ force fields. Under this method, cell wall particles repel the fluid via these LJ type interactive forces. 24,25,48 In order to make the repulsion forces more effective, a set of virtual particles are also incorporated, which are mass-less and artificially placed among the cell wall particles. 20,22,34 It is noteworthy that a moisture-domain-based approach has been employed to represent the drying mechanism rather than adopting a time-domain-based approach. The main reason behind this was to reduce the computational cost associated with the time-domain-based approach. 20 As an outcome of this, the whole drying mechanism for a single apple cell has been broken down into a number of separate steps. This approach has been successfully used before by Karunasena et al. in a numerical modelling study for 2-D cellular scale plant cell drying. 25 A set of representative dried cell conditions were selected having separate moisture content and turgor pressure values for characterising the entire drying operation. The developed single cell numerical models were then employed to simulate those selected dryness states. Independent simulations were conducted via initiating the model with relevant  25 The turgor pressure of a fresh apple cell (X/X 0 = 1) was taken as 200 kPa. 3 Accordingly, the dryness states of X/X 0 = 0.9, X/X 0 = 0.7, X/X 0 = 0.5, X/X 0 = 0.3 and X/X 0 = 0.1 were simulated at turgor pressure values of 180, 140, 100, 60 and 20 kPa, respectively. In each simulation, the value of the osmotic potential (P) was maintained equal, but at the negative value of the corresponding turgor pressure value (P i ). After the achievement of steady state, the eventual cell physical properties and other geometrical parameters were used to represent the model predictions. Then, these model predictions were compared in detail with relevant experimental findings. 25 For the experiments, Royal Gala apple (Malus Domestica) procured from Brisbane (Australia) were used with appropriate sample preparation techniques. 65 For drying, a convective air dryer (Excalibur's five-tray dehydrator, USA) was used which has an electric heater and a fan to produce a controlled hot air flow across the samples throughout the experiments. The air temperature can be conveniently adjusted using a thermostat. For this particular series of experiments, an air temperature of 70 1C with a constant hot air velocity was maintained. The samples were introduced only after the dryer reaches the steady state condition following the initial warming-up cycle. It should be noted here that in order to capture the gradual morphological changes of the samples in each drying experiment, the same apple sample was used and intermittently observed using the microscope.
To quantify the cellular deformations in the model predictions, four geometrical parameters were employed: cell area (A), diameter † (D), perimeter (P) and roundness ‡ (R). The variation of these parameters were analysed with the dry basis moisture content X (mass water /mass dry solid ). To enable a better comparison, these parameters were normalised (A/A 0 , D/D 0 , P/P 0 and R/R 0 ) by dividing the current value of the parameter by the initial value at the fresh cell state (X 0 , A 0 , D 0 , P 0 and R 0 ). The model was established as a C++ source code and executed in a High Performance Computer (HPC). Algorithms of existing SPH source codes based on FORTRAN 48,53 and C++ 20 were referred to in the development of the C++ source code for the 3-D cell.
For the visualisations, Open Visualization Tool (OVITO) 66 was used. To make these quantifications in the model predictions, inbuilt image processing techniques available in the ImageJ software was used. 67,68

Computational accuracy of the model
To determine the computational accuracy of the developed SPH-CG single cell model, a method adopted by several previous researchers 20,24,25,34 in their SPH-DEM cell models was employed. In this method, the model prediction for the average hoop directional inter-particle force in the cell wall was compared with the theoretically expected value. § 34 The percentage model consistency error ¶ was calculated based on that difference. It assumes an equilibrium between the cell fluid turgor pressure and the tension in the cell wall according to Newton's third law and Young-Laplace law. For all the cell categories simulated in this study, the percentage model consistency error was maintained below 1%. It should be mentioned here that, this computational accuracy is significantly higher than the recently reported SPH-DEM 3-D fresh cell models. 34

Simulation of deformations during drying
As it has been described in the computational implementation section (Section 2.5), for each plant cell category, a number of cell dryness states were simulated starting from a fresh state (X/X 0 = 1.0 and P T = 200 kPa) to an extremely dried state (X/X 0 = 0.1 and P T = 20 kPa). At this extremely dried state, 90% of the total moisture content gets removed from the cell model. Even the most recent grid-based numerical models which study the deformations of fruit tissues during dehydration 69 are not capable of simulating a moisture content reduction of this degree. 14 This highlights the significant ability of meshfree-based methods to model multiphase phenomena and large deformations more effectively.
In Fig. 4, the model predictions are visualised for all the plant cell categories at all the dryness states. To make the model predictions further illustrative, sectioned side views corresponding to each dryness state are presented in Fig. 5. Particles have been colour-coded to easily distinguish between fluid (blue) and wall (green) particles. As it is visible in Fig. 4  and 5, the sizes of the cells gradually decrease with the increasing degree of dryness. When these predictions were compared with the results produced by 2-D SPH-DEM-based cellular drying simulations, [20][21][22]25,35,36 it could be noted that the 3-D model depicts the true drying scenario in a more detailed and realistic manner. Furthermore, it could be observed that the apple cells and grape cells have similar sizes while carrot cells are smaller in size. The potato cells are the largest. In all simulations, when comparing the initial particle configuration of the cells (i.e. Fig. 4(a) and 5(a)) and the fresh states (i.e. Fig. 4(b) and 5(b)), it could be observed that the cells have inflated (mainly due to the turgor pressure), resembling the turgid nature that is frequently observed in fresh cells and tissues. The dried cells have shrunk as seen in Fig. 4(c)-(g) and 5(c)-(g), approximating the shrinkage behaviour of real plant cells when subjected to drying.  It could be qualitatively observed that there are differences in shrinkage among different categories of food-plant cells. The degree of shrinkage demonstrated by the potato cell is relatively lower compared to that of the carrot cell. A major reason behind this could possibly be the difference of 'a' value used in the wall contraction force field (F c kj ) (eqn (22)). It should be noted here that the cell fluid particle number remains constant in all simulations and only the particle mass is reduced during drying following the fundamental 2-D meshfree based cellular drying models. 20,25 It is also noteworthy that as this is a basic single cell scale model, a number of properties across the model have been considered to be isotropic for the sake of simplicity. For example, cell wall properties like permeability, thickness, shear modulus and damping coefficient could be given. However, the computational modelling approach discussed here have the capability to incorporate any local variations of such characteristics due to the particle-based methodology. In such occasions, different sections of the cell model can be assigned different values of properties such as wall permeability. The heterogeneity of the cellular structure could be addressed in a similar manner.
Next, the post processed model predictions were qualitatively compared with the results of experimental investigations on real plant cellular structures of apple, grape, carrot and potato. For this purpose, experimental findings from this study 65 and from investigations reported in literature 3,26 were used. In Fig. 6, microscopic images from an experimental investigation on the apple cellular structure during drying (involving 3-D digital light microscopy and image analysis techniques) are shown. 65 It should be noted that these images exhibit the behaviour of a single apple parenchyma cell during the drying process. During the microscopic investigations, the same cell has been imaged at all different dryness states. This makes those images and subsequent results highly comparable with the developed SPH-CG single cell model. When Fig. 4-6 are compared, it is evident that there are similarities in the morphological changes during drying.

Quantifying the morphological behaviour of the models through geometrical parameters
The simulated deformations of the dried cell models for all food-plant categories were quantified using normalised cell area (A), feret diameter (D), perimeter (P) and roundness (R) as discussed in Section 2.5. These quantitative results were then compared with experimental findings from literature for apple, grape, carrot and potato cellular structures. 3,26 At the same time, they were compared with the model predictions presented by Karunasena et al. 25,35 in their SPH-DEM 2-D meshfree-based single cell drying models.
3.3.1 Morphological behaviour of the apple single cell model during drying. Fig. 7 shows the simulated apple cells at different dryness states. Corresponding results of the geometrical parameter comparisons are presented in the graphs of Fig. 8. The variation of each geometrical parameter has been plotted against the normalised moisture content (X/X 0 ) during drying. The variation of normalised area is shown in Fig. 8(a). Overall, the model predictions of this study demonstrate an agreement with the experimental findings. 3  tissues (experimental findings), interacts closely with its neighbouring cells by default. However, in these single cell computational models, such inter-cell interactions are not present. This affects a 3-D single cell model in a more severe manner than a 2-D single cell model because of the increased number of degrees of freedom in 3-D.
In case of diameter variations (D/D 0 ), the agreement between the model outcomes and the experimental findings show a similar trend (Fig. 8(b)) to the area variation. The model predictions from this study show an agreement with the experimental values from literature as well as with the model predictions from the previous study from Karunasena et al. The agreement is stronger in the high moisture contents. At low moisture levels, the model predictions from this study tends to slightly deviate from the experimental findings. For normalised perimeter, the comparison follows a similar trend to those of the area and the diameter (see Fig. 8(c)). Experimental findings and the previous 2-D model predictions show a reduction in the cell perimeter with decreasing moisture content. The results from this study demonstrate a favourable agreement with that behaviour. The model predictions show a decreasing cell perimeter which is in close coincidence with those values from the experiments and 2-D modelling attempts.
When it comes to cell roundness (R/R 0 ), experimental findings show an unchanged value which is close to 1 throughout the drying process. As seen in Fig. 8(d), the predictions for the roundness in the 3-D models of this study follow a similar behaviour in the higher moisture levels before exhibiting a  considerable discrepancy in the low moisture contents (X/X 0 r 0.5). That is, the roundness predicted by the model in this study becomes lower with the decreasing moisture content. This deviating behaviour is evident in the recently reported SPH-DEM 2-D plant cell drying models in a similar or higher degree. 3.3.2 Morphological behaviour of potato cell model during drying. Fig. 9 shows the simulated potato cells at different dryness states while corresponding results of the geometrical parameter comparisons are presented in the graphs of Fig. 10. It is evident that the potato cell undergoes a relatively lower degree of morphological variations during drying. This is reflected through very small variations occurring in the geometric parameters (i.e. normalised area, diameter, perimeter and roundness). A similar trend could be observed in the previously reported investigations in literature on cellular scale morphological variations of potato during drying. 35,70 As seen in Fig. 4, potato cells are significantly larger than the rest of the cell categories considered in this study. This is in correspondence with the parameters used in determining the strength of the cell wall contraction force field (F c kj ). The SEM images reported in literature 35 agree with these model predictions. The relatively larger size of the cells as well as the relative isotropic/anisotropic nature of the overall cellular structure could have an effect on this comparatively lower degree of morphological variations.
It could be seen from Fig. 10(a) that there is a favourable agreement between the model predictions of this study and the experimental findings on potato cellular structure in terms of normalised area variation. 62 Similar to the apple cell model behaviour discussed in the Section 3.3.1, there is a small deviation between the model predictions and experimental findings at the extremely dried states (X/X 0 r 0.3). The agreement between the model predictions and experimental findings demonstrate a similar behaviour for the cell diameter and the perimeter as seen from Fig. 10(b) and (c). There is a matching behaviour of the results apart from a slight deviation towards the lower end of the moisture domain. The normalised roundness variation of the cell model is also in close agreement with the experimental cell roundness variation. 70 The value of roundness stays close to a value of 1.0 throughout the process. This is also evident from the post-processed model predictions seen in Fig. 9 where the cell exhibits a circular shape at all dryness states. In addition, the results of the 2-D SPH-DEM cellular drying model 25 agree with the findings from this 3-D SPH-CG study.

3.3.3
Morphological behaviour of grape cell model during drying. Fig. 11 shows the simulated grape single cells at different dryness states. Corresponding geometrical parameter variations are presented in the graphs of Fig. 12. As seen there, the 3-D SPH-CG model predictions from this study have been compared to experimental findings obtained through a stereo microscopy investigation. 5 However, these microscopic investigations have considered drying only up to a normalised moisture content (X/X 0 ) of 0.6 as seen from the graphs in Fig. 12. In the 3-D model predictions, there is a gradual shrinkage behaviour which lead to a favourable agreement with the experimental findings. It could be observed that the grape cells go through a slightly higher degree of shrinkage compared to apple and potato cells. For an example, the image processing results show that the values of normalised area (A/A 0 ), diameter (D/D 0 ) and perimeter (P/P 0 ) for grape cells at an extremely dried state of X/X 0 = 0.1 are respectively 0.77, 0.88 and 0.88. The corresponding values of those parameters for apple cells are respectively 0.81, 0.90 and 0.90. The degree of shrinkage of potato cell models are lower than that of both apple and grape cells. One possible reason causing this difference is the variation of physical characteristics.
Another key conclusion derived through the observation of Fig. 12 is that the 3-D SPH-CG numerical cell model developed under this study has the ability to model the cells throughout a larger moisture domain compared to recently reported 2-D SPH-DEM cellular drying models. 35 2-D models have considered a moisture removal from the cellular system only up to approximately X/X 0 = 0.3 while the 3-D model of this study has been able to model moisture removal up to X/X 0 = 0.1. This provides a hint about the stability of the adopted computational coupling between SPH and CG. There is a slight deviation between the model predictions and experimental findings 5,63 especially in the range of X/X 0 = 0.8 to X/X 0 = 0.6. The 2-D SPH-DEM model predictions also show this mismatch with the experimental findings to a certain degree as depicted in Fig. 12(a)-(c) respectively. One of the reasons behind this could possibly be the significantly smaller moisture content range addressed in the experimental studies. In addition, the lower magnification used in the microscopy could potentially be a contributing a reason for the discrepancy. 5,35,63 As seen from Fig. 12(d), the roundness of the grape cell model stays almost constant throughout the entire drying process. This agrees with the conclusions drawn from the visualised numerical results in Fig. 11 where the cell maintains its circular shape even at very low moisture contents. This agrees with the findings from relevant experimental investigations 5,63 as well as 2-D model predictions. 35 However, 3-D model predictions do not demonstrate a significant cell wall wrinkling behaviour during drying. The absence of the cell-cell interactions could be one probable cause behind this. This again highlights the necessity to extend these single cell models into the multiple-cell levels to improve the efficacy of approximation.
3.3.4 Morphological behaviour of carrot cell model during drying. Fig. 13 shows the simulated carrot cells at different dryness states while geometrical parameter variations are presented in the graphs of Fig. 14. There is a significantly higher degree of shrinkage in the carrot cell model compared to the other three types of food-plant cells. To validate this observation, model predictions at the extremely dried state of X/X 0 = 0.1 has been tabulated in Table 3 for all four food-plant cell categories. This provides conclusive evidence that the carrot cell model goes through a relatively higher degree of shrinkage compared to the other three types of cells. Potato cell model corresponds to the lowest degree of shrinkage while apple and grape cells demonstrate moderate shrinkage characteristics. In addition, information in Table 3 suggests that there could be a correlation between the degree of shrinkage  and the relative size of the cell. This phenomenon has the potential to be further investigated in detail during future investigations. 3-D model predictions for carrot cells during drying have been compared with experimental findings on real carrot microstructure 35,59 (see Fig. 14). It is evident that there is an agreement between the 3-D model predictions and the experimental findings. The agreement is favourable even at low moisture content values (X/X 0 r 0.3). This is a unique positive characteristic observed in the carrot cell model because the other three types of cell models demonstrated a deviating behaviour from the experimental findings at the low end of the moisture domain. This 3-D SPH-CG numerical cell model has the ability to model the carrot cells in a stable manner throughout a larger moisture domain compared to the 2-D SPH-DEM cellular drying models 35 (see Fig. 14(a)-(c)). The 3-D cell model from this study has considered drying process up to a moisture content value of X/X 0 = 0.1 while the 2-D models have considered up to X/X 0 = 0.3. This occurrence could be observed in potato and grape cell models as well. The added stability of the models in this study could be arising due to the 3-D nature of the developed computational model.
In all the food-plant cellular drying models considered in this study, there are notable agreements with the  View Article Online experimental findings. However, there is further room for improving the performance of this computational modelling approach. As mentioned previously, one of the key limitations in this investigation is not taking cell-cell interactions into account in the physical description of the plant cellular structure. This could possibly be one of the major reasons behind the disagreements between model predictions and experimental findings. Also, the morphological behaviour of the cell wall plays a critical role towards the overall morphological behaviour of the cell model. Incorporating an enhanced cellular turgor pressure relationship could improve the morphological description of the cell wall model (compared to the realistic morphological behaviour). This could lead to overall performance improvements of this meshfree-based cellular drying model.

Conclusions
In this investigation, a three dimensional (   potato and grape cells. In carrot cells, there is a favourable agreement even at very low moisture content values. This is a positive characteristic observed in this novel computational model relative to the previous 2-D SPH-DEM plant cell drying models. In addition, single cell-based drying experiments have been conducted to obtain 3-D details on the morphological changes of the plant cellular structure during drying. This experimental information has been a valuable source of validation data for the 3-D SPH-CG model predictions. Improved stability of the computations is another key characteristic that was exemplified by this 3-D meshfree-based numerical model. In addition, it could be concluded that employing a CG approach to model the plant cell wall membrane has established a compatible numerical coupling with SPH. To improve the validity and applicability of the developed numerical model, extending this single cell approach to tissue level will be an important future work. In doing so, using enhanced computational algorithms might be necessary to cope with the increased computations in such complicated systems. 71 In addition, the heterogeneity and anisotropic nature in the plant cellular structure could be addressed in future attempts through introducing localised property variations (e.g. cell wall permeability etc.). Incorporating temperature variation related effects is a further area for improvement. The proposed 3-D SPH-CG computational modelling approach could potentially be applied in modelling the morphological changes of animal cells during dehydration processes.

Conflicts of interest
There are no conflicts of interest to declare.