Deformation, mechanical breakdown, and flow of soft solid foods in the stomach

Abstract

Mechanical changes to solid foods in the stomach are crucial aspects of processes for regulating nutrient bioavailability, satiety, and glycaemic response after a meal. However, the underlying mechanisms are poorly understood. This study uses a Smoothed Particle Hydrodynamics (SPH) model to study the deformation and mechanical breakdown of solid beads in a liquid medium within a realistic three-dimensional representation of the stomach geometry. The model incorporates peristaltic contraction waves, including Terminal Antral Contractions (TACs), which are high-amplitude and high-speed travelling occlusions observed in the distal region of the stomach. The solid beads are modelled using elastic-plastic (EP) and elastic-brittle (EB) constitutive laws. Results show that the stomach wall can induce significant compression and fragmentation in the solid beads through direct contact, and the accompanying fluid flow contributes towards further mechanical change. An originally spherical EP bead closest to the pylorus is extruded into a thin cylindrical shape, generating 15 fragments with a 5% higher surface area, before being propelled away from the TAC region. A model-parameter sensitivity analysis shows that an increase in yield stress substantially reduces the fragmentation but not the elongation. An EB bead near the pylorus deforms less but fractures into 235 small fragments and a large chunk, leading to an overall 12% higher surface area. The EB bead remains near the pylorus, fracturing further over multiple peristaltic cycles. Increased fracture strength, represented by a higher threshold strain, significantly reduces the surface area change and the number of fragments generated by the wall contractions. These results show how a coupled biomechanics-fluid-elastic-plastic-fracture model can be used to investigate the mechanical breakdown of solid foods in the stomach.

---

1 Introduction

The deformation and mechanical breakdown of food in the stomach are key components of life-critical stomach functions such as nutrient release and control of satiety. These processes lead to size reduction of ingested solid food, which aids in gastric processing and is critical for nutrient liberation. For example, the degree of size reduction of starch-based foods in the stomach can control the surface area available to amylase enzyme in the small intestine, thus influencing sugar release.^1,2 The release of sugars in the small intestine, in turn, affects gastric emptying through inhibitory feedback mechanisms such as the ileal brake.^3,4 In addition, gastric emptying is regulated by the size of food fragments in the stomach.^5–7 Solid fragments must be reduced to dimensions below 3 mm to facilitate their passage through the pylorus to the duodenum.^6,8 This process of emptying of solid digesta contributes to the regulation of antral distension, a factor closely linked to feelings of hunger and satiety,⁹ which in turn can affect outcomes for metabolic conditions such as diabetes and obesity. Overall, the mechanical breakdown of food in the stomach and subsequent gastric emptying are instrumental in controlling satiety, blood glucose absorption and glycaemic index through several interconnected processes. However, characterising how a food item will deform and breakdown in the stomach based on its structure and material properties remains a challenge.

Currently, the in vivo characterisation of solid food deformation and mechanical breakdown relies on a limited set of practically and ethically difficult experimental methods. Experimental techniques such as magnetic resonance imaging^10,11 (MRI) or aspiration using nasogastric tubes^12,13 lack the ability to continuously characterise important mechanical processes such as solid deformation and fracture, or the relationship of these processes to hydrodynamic flows and contractions in the gastric wall. MRI studies are further limited by their spatial resolution, which means that small but still significant particles are usually not identified.¹⁴ Similarly, structures such as a large food mass in the form of an intact solid particle may not necessarily be distinguishable from an agglomerated semi-solid mass of similar size. The interpretation of MRI images is often complicated by noise and other acquisition artifacts.¹⁵ In general, experimental studies are also constrained by ethical considerations, such as subject recruitment from specific demographics or individuals with pre-existing conditions. Despite these limitations, MRI has been utilised to roughly estimate the size, shape, and number distribution of food fragments and ingested beads,^{10,11,14–16} using images that are collected at specific time intervals. For example, in Marciani et al.'s MRI study,¹⁶ subjects ingested liquid meals along with agar gel beads with a range of fracture strengths ranging from 0.15 N to 0.90 N, with their breakdown being monitored using MRI images collected at 15-minute intervals. The results revealed that agar beads softer than 0.65 N were rapidly broken and emptied at a nearly same rate as the liquid. In contrast, beads with fracture strength exceeding 0.65 N had a slower exit from the stomach. The study concluded that the maximum force exerted by the stomach for solid food breakdown is approximately 0.65 N.

In vitro models offer an alternative for studying the mechanical breakdown of foods^17–19 without many of the practical and ethical challenges associated with in vivo methods. However, they also present a different set of practical challenges such as accurately reproducing little understood and poorly characterised in vivo conditions, including the physical forces acting on the gastric digesta due to peristaltic movements in the stomach wall. For example, the Human Gastric Simulator (HGS) utilises a series of mechanically driven rollers to simulate peristaltic movements on a deformable vessel made of latex rubber.¹⁹ When compared with in vivo emptying in growing pig stomachs, HGS shows similar trends in gastric breakdown rate and solid emptying kinetics.²⁰ Over the last decade, several groups have also attempted to create in vitro models with more realistic descriptions of the gastric morphology and geometry.^21–24 A series of in vitro models have been developed based on silicone-based J-shaped stomach, modelled after a cadaver-stomach, including a rope-driven in vitro human stomach²¹ (RD-IV-HSM), and a subsequent dynamic in vitro human stomach²³ (DIVHS) model driven by a set of eccentric wheels and rollers. The latest generation of these models, known as dynamic in vitro human stomach system²⁴ (DHS-IV), uses a rolling extrusion device to simulate peristaltic contractions in the gastric wall including terminal antral contractions (TACs). TACs are peristaltic contraction waves in the terminal region of the antrum (Fig. 1) that accelerate while approaching the pylorus, along with a rapid increase in amplitude and pulse-width.^25–29 DHS-IV also incorporates a pneumatically controlled pyloric sieving valve, where the valve opening size and frequency can be adjusted to regulate the size of the solid fragments emptied from the stomach. The device parameters can be optimised to reproduce theoretical emptying profiles of various food samples, while particle size distribution of the retained digesta can be measured to assess food breakdown. Despite these advancements, further refinement of the in vitro models is needed. For instance, the pyloric sieving valve mechanism in DHS-IV does not accurately reflect the contraction dynamics of the human pyloric sphincter observed in endoscopic videos.³⁰ Moreover, current in vitro devices lack integrated velocimetric measurements of digesta flow or capability to track movement of solid particles, limiting their ability to identify various breakdown mechanisms and quantify their actions. Continuous monitoring of the digesta flow can also reveal spatiotemporal variations in gastric processing conditions, such as density, viscosity, and pH across the different stomach regions.

Fig. 1 Resting stomach model geometry and details of ACWs: (A) the surface mesh of the stomach is shown in its resting, non-contracted state. The dotted centreline indicates the length of the stomach along which ACWs are applied. (B–D) Plots depict velocity, pulse width, and contraction amplitude of an individual ACW as a function of its location along the stomach centreline. (E) The opening and closing dynamics of the pyloric sphincter, regulated by the location of the nearest ACW to the sphincter.

In silico models do not have the same practical difficulties encountered with in vivo and in vitro methods, although their effectiveness can be limited by the model assumptions made, the calibration data available, the capabilities of the numerical methods, the numerical resolution, and the power of the computational infrastructure. For example, in silico models provide richer quantitative data and offer the advantage of continuous and precise measurement of model variables, as well as the ability to explore variations in model parameters, which are impractical or unethical in real life experiments. In silico models are also not subject to sampling biases, such as subject recruitment from specific demographics, and they do not require large sample sizes typical in many in vitro devices.¹⁷ However, previous in silico modelling studies that focus on the flow of solid digesta in the stomach use several assumptions and approximations that limit their utility in studying mechanical breakdown. Most treat solid food fragments as continuous fluid-like phases with specified densities or mass fractions.^31–34 While computationally efficient, this approach cannot capture insights into the mechanical interaction of larger food fragments with gastric wall contractions and fluid flow, including details of their movement, deformation, and fracture. Related to this, Seo et al.³⁵ and Lee et al.³⁶ used an immersed boundary method with unstructured mesh to model the movement of an orally-ingested pill inside the stomach. However, the model uses a small, single solid body with wall contact modelled as a spring, so it does not account for the mechanical breakdown of the solid through deformation or fracturing. Only one published computational study attempts to model the deformation and fragmentation of solid food in the stomach using realistic constitutive laws and material parameters.³⁷ The model approximates the stomach geometry as a straightened axisymmetric vessel without any curvature and predicts the deformation and fragmentation of a thin soft solid food due to backward extrusion and indentation by antral contraction waves (ACWs). However, the model has limited applicability due to an absence of liquid environment and assumptions of zero gravity, as well as an unphysically idealised gastric morphology. The simplicity of the current models underscores the need for more sophisticated models integrating a broader range of aspects of the gastric biomechanics within a less idealised stomach geometry.

In this study we present an update to a previously proposed Smoothed Particle Hydrodynamics (SPH) model of gastric flow^38–40 which now includes representations of solid foods and their deformation and fracture, all coupled to the dynamics of the liquid medium and stomach walls. An improved representation of gastric wall movements has been developed, based on MRI and high-resolution electrophysiological activity data.²⁶ The deformation and fracture of elastoplastic (EP) and elastic-brittle (EB), freely moving, neutrally buoyant initially spherical solid beads are simulated. These are used to identify and analyse the biomechanical drivers underlying the deformation and fragmentation of soft solid foods. This study presents the first model of mechanical breakdown and fragmentation of solid food in a mixed solid–liquid digesta, incorporating realistic material models and gastric geometry and ACW wall deformations.

2 Computational methods

2.1 Smoothed particle hydrodynamics

SPH is a mesh-free Lagrangian approach for solving partial differential equations. The method represents material volumes as a set of particles that move with the fluid flow or solid deformation.^41–44 This approach enables efficient simulation of complex coupled systems involving fluid dynamics and large scale solid deformations, such as the breakdown of solids in the digestive tract.^42,45–48 In SPH, material properties are interpolated from the neighbouring particles using a smoothing kernel, allowing the Navier–Stokes or solid deformation equations to be expressed as a set of ordinary differential equations representing conservation of mass and momentum.^42,49 The simulations performed here were performed using the CSIRO SPH code which has been extensively used for a broad range of industrial and biophysical applications.^41,42

2.1.1 SPH formulation for fluids. Based on Monaghan,⁵⁰ the Navier–Stokes equations can be reduced to the following ordinary differential equations for conservation of mass and momentum:

Here, for particle a, ρ_a is density, m_a is mass, P_a is local pressure, and μ_a is dynamic viscosity. r_ab denotes the distance between the particle a and another particle b, while their relative velocity is given by v_ab. W_ab is a cubic-spline interpolation kernel function that is evaluated for r_ab. The kernel function is discussed in detail in Monaghan.⁴⁹η is a small number intended for avoiding singularities from small denominators. ξ is a normalisation constant for kernel function. Finally, g stands for gravitational acceleration.

A quasi-compressible SPH formulation is utilised, which relates the particle density ρ with fluid pressure P by the equation of state for a weakly compressible fluid.

where c is the speed of sound and ρ₀ is the reference density. In addition, γ is a material constant with a value equal to 7.0 for fluids with water-like properties.⁵¹ To minimise density variations from compressibility effects to the order of 1% or smaller, a match number equal to 0.1 is used.

A second-order predictor–corrector integration scheme is used,⁵⁰ and a time step is chosen to satisfy the Courant condition with a modification for the presence of viscosity:

where h is the smoothing length of the kernel. This choice of δt guarantees stability of the explicit numerical integration.

2.1.2 SPH formulation for elasticity in solids. The conservation of mass equation for solid materials is essentially the same as for fluids (1). Conservation of momentum in the SPH formulation leads to the following acceleration equation.where σ_a denotes the stress tensor for particle a, Π_ab represents the viscosity terms for both shear and bulk viscosity, and I is the identity tensor. A linear elastic model is used for the solid rheology, given by the equation of state which expresses pressure as a function of density and reference density.

P = c²(ρ − ρ₀) (6)

The speed of sound c in a solid material is calculated as a function of the bulk modulus K.

The stress tensor σ is partitioned into a pressure term and a deviatoric stress term with tensor S.

σ = −PI + S (8)

Based on Gray et al.,⁵² the evolution of S can be expressed in component form as:

where έ is the strain rate tensor, δ^ij denotes the Kronecker delta, Ω stands for the Jaumann rotation tensor, and the shear modulus is given by G. The indices i, j and k refer to the three orthogonal directions in 3D space. The strain rate and Jaumann rotation tensors are formulated in SPH as:

Similar to the integration scheme for fluids (4), a second-order predictor–corrector (explicit) integration scheme is employed with time step chosen to satisfy the Courant condition.

The SPH method can display tensile instabilities, particularly when used for modelling elastic solids.⁵³ Therefore, the principal stress correction technique proposed by Gray et al.⁵² is used to suppress such instabilities for the solids. The effectiveness and accuracy of this method have been previously validated by computing stress distributions under a uniaxial loading scenario and comparing the results with those obtained from a mesh-based Finite Element (FE) model.⁵⁴ This demonstrated that the correction technique accurately predicted the stress distributions, while fully inhibiting numerical fracture arising from tensile instabilities. A tensile instability correction coefficient of 0.3 is used in the current simulations.

2.1.3 SPH Formulation for elastic-plastic solids. EP solids behave elastically at low strains but once the compressive stress exceeds the yield strength (YS), irreversible plastic deformation occurs. In the SPH formulation, plastic deformation is tracked by the radial return plasticity model.⁵⁵ Assuming an elastic response, a trial deviatoric stress S^ij_Tr is determined at the beginning of each timestep. The resulting change in plastic strain ε^P is given by:where σ_vm denotes the Von Mises stress, σ_y the yield stress, and the hardening modulus is given by H. The plastic strain is accordingly incremented.

ε^P = ε^P + Δε^P (13)

The increment in σ_y at the end of the timestep is calculated as:

Δσ_y = HΔε^P (14)

The deviatoric stress at the end of the time step is determined based on the trial deviatoric stress as:

S^ij = r_sS^ij_Tr (15)

where r_s denotes the radial scale factor:

YS, the initial value of σ_y, is provided as an input to the model. Segments of the EP solids that experience substantial plastic deformation may be displaced far from the main body, eventually losing contact. These detached segments are considered as distinct fragments.

2.1.4 SPH Formulation for elastic-brittle solids. EB beads behave elastically until the applied strain exceeds a threshold strain (TS), which causes accumulation of damage that inhibits the transmission of tensile stress, eventually leading to macroscopic fracturing. EB solids are modelled using a continuum damage law derived from the model of Grady and Kipp.⁵⁶ Based on Benz and Asphaug,⁵⁷ each SPH particle is assigned a damage parameter, D, which has a value between 0 and 1, specifying the extent of local fracture. A non-zero value of this parameter leads to inhibition of the transmission of tensile stress across all three dimensions. The SPH formulation of the damage evolution is given by:where ε denotes the tensile strain, and α is a constant given by:

Here k is the flaw density, m denotes the Weibull modulus, and C_g is the crack growth speed. This model was originally developed for one-dimensional problems only.⁵⁷ In the current study, an equivalent 3D formulation developed by Cleary and Das⁴⁴ is applied, using an effective tensile strain from Melosh et al.⁵⁸

where α_max is the maximum principal stress. Evolution of damage only occurs if ε is larger than the TS.

2.1.5 Contact mechanics. An improvement to the previously used SPH model is the inclusion of explicit contact mechanics. Previously the SPH elastic deformation equation was solved for normal compressive behaviour across contacts between solid bodies and between solid bodies and the stomach walls. Here, a linear spring-dashpot interaction is instead used between any pair of SPH particles that are on different bodies. The forces are calculated in a canonical frame normal to the contact surfaces. This enables discontinuous behaviour such as surface friction and use of contact stiffness that are independent of the material moduli. It also means that an explicit coefficient of restitution can be used which can be calibrated to the inelasticity of the material. In contrast, a pure SPH method relies on the artificial viscosity (that is present to provide some numerical stabilisation) to provide collisional dissipation and which cannot be calibrated to specific materials. The contact force model is the same as used in the Discrete Element Method. See Thornton et al.⁵⁹ for details about various contact models that are available for inelastic contact mechanics. The application to distributed SPH contacts follows the approach given in Cummins and Cleary⁶⁰ and is similar to that of Vyas et al.⁶¹

2.1.6 SPH solid boundaries. To model solid boundaries, 3D triangular surface meshes are utilised. The nodes of the meshes are represented using boundary SPH particles whose positions are updated at each time step as a result of wall movements. These boundary particles are subjected to a penalty force along a normal direction using a Lennard-Jones style form, based on the orthogonal distance between adjacent moving SPH particle and the surface boundary particle.^41,50 For the adjacent boundary-moving SPH particle pairs, the pressure force terms in the conservation of momentum eqn (2) and (5) are replaced by the penalty force. The particle pairs are also included in the summations for the viscous stress in these equations, leading to an application of no-slip boundary conditions in directions tangential to the solid surfaces. Further details about the SPH method employed in the current model can be found in previous studies.^38,41,42,45

2.2 Stomach geometry and deformations

2.2.1 Stomach geometry. The stomach geometry is represented by a closed solid boundary (Fig. 1A) created using Autodesk Maya (Autodesk, San Francisco, USA), and based upon published imaging data. The large scale features of the 3D triangular surface mesh are based upon MRI images of the stomach in a seated posture.⁶² These images were acquired from healthy subjects who underwent scanning after ingesting a mixed solid (150 g) – liquid (150 mL) meal in a fasting state. High-resolution computed tomography gastroscopy (CTG) images^63,64 were used to design the lower oesophageal sphincter and the cardia of the model stomach.

2.2.2 Large-scale stomach wall deformations. The deformations in the 3D surface mesh are kinematically prescribed as moving boundary conditions to represent the movements of ACWs (Fig. 1). The ACWs originate at the proximal antrum at intervals of 20 s and conclude at the pyloric sphincter, travelling approximately 142 mm over a period of 54 s. The waves have temporal variation in amplitude (Fig. 1B), pulse width (Fig. 1C), and speed (Fig. 1E), which depend on a wavefront's position relative to the pyloric sphincter. Notably, after 40 s, when a wavefront is approximately 30 mm away from the sphincter, the wave enters the terminal antrum region. Subsequently, the amplitude, pulse width, and speed of the wave rise sharply, leading to the formation of a TAC. Compared to the slower peristaltic waves in rest of the antrum, TACs appear as near-simultaneous segmental contractions of the terminal antrum. The ACW parameters shown in Fig. 1B–D were derived from Ishida et al.,²⁶ who incorporated MRI data from Pal et al.²⁸ Pal et al. determined the ACW parameters by fitting a cubic spline through selected points in a 2D MRI movie obtained from a single subject following consumption of a nutrient-rich meal. Additionally, Ishida et al. integrated velocity data from an experimental study by Berry et al.,²⁵ where high-resolution electrophysiological activity data was used to quantify the ACW parameters in fasted stomachs of patients undergoing hepatobiliary or pancreatic surgery. Currently, there is not enough data to specify changes of these motility patterns in responses to stomach volume or the nature of ingested food.

To implement the travelling ACWs on the 3D surface mesh, the antral and pyloric regions of the mesh are subdivided into 55 evenly spaced radial segments. A curved centreline is defined, passing through the centres of these rings to guide the ACWs. These rings are rigged and programmed to radially contract according to the magnitude of the ACWs at their respective centres. Each node in the surface mesh undergoes deformation based on a weighted average influence of the two nearest rings. This weighted average is computed using a Gaussian smoothing operation to ensure continuous and smooth deformations across the mesh.

2.2.3 Pyloric geometry and sphincter dynamics. Given the pivotal role of TAC and pyloric sphincter dynamics in gastric breakdown of solid digesta,^37,65,66 the pylorus region of the stomach mesh was specifically developed using multiple independent imaging sources. The pylorus has two parts, the pyloric antrum, which connects to the stomach, and the pyloric canal, acting as an opening to the duodenum through the pyloric sphincter. Images from dissected cadavers⁶⁷ show a rounded appearance of the pyloric canal, where both the greater and lower curvatures of stomach converge towards the pyloric sphincter along a convex profile (Fig. 2A). This rounded appearance of the pyloric canal has been integrated in the current 3D stomach mesh. The cadavers images⁶⁷ were further used to determine the location of the pyloric sphincter relative to the rest of the model pyloric canal.

Fig. 2 Pyloric sphincter geometry: images on the left are experimental images sourced from existing literature, illustrating different degrees of pyloric sphincter opening alongside the surrounding antro-duodenal region. On the right, corresponding 3D renderings of the stomach model around the pyloric sphincter are presented. This also includes the model inputs of sphincter diameters used for generating the 3D images. (A) A photograph of dissected pyloric sphincter (PS) from a cadaver stomach, where mucosae and submucosae layers have been removed, exposing the muscle bundles (adapted from⁶⁷ with permission). Here cD denotes the duodenum, PC denotes pyloric canal, while PS indicates the location of the pyloric sphincter. (B) Image of a nearly closed pyloric sphincter and the surrounding antral region captured during endoscopic examination of a patient with epigastric pain (adapted from³⁰ with permission). (C) Fluoroscopic (X-ray) image of a closed pyloric sphincter and the surrounding antro-duodenal region, taken from a barium-filled stomach in a healthy subject (adapted from⁷¹ with permission). The model was built to closely emulate these features.

Based on data from Ishida et al.,²⁶ the pyloric sphincter was prescribed to periodically close and open in coordination with each ACW approaching the pylorus. The three timepoints when the pylorus is opening, fully opened, and closing were synchronised against the times required for an individual ACW to reach three distinct distances from the pylorus (Fig. 1E). This location-based synchronisation between pyloric behaviour and ACWs is consistent with the observation of Indireshkumar et al.⁶⁸ that transpyloric flow is blocked when ACWs reach a “zone of influence” 20–30 mm away from the pylorus. When maximally open, the resting diameter of the pyloric sphincter is 10 mm. This value is consistent with measurements from previous experimental studies assessing the internal diameter of the sphincter in healthy humans.^69,70

The model of the pyloric sphincter dynamics is evaluated using experimental images captured during various phases of sphincter closure. Fig. 2B provides a side-by-side comparison between a close-up view of the partly closed pyloric sphincter from an endoscopic video³⁰ and the corresponding representation in the stomach mesh. The relative position of the sphincter with respect to the remainder of the stomach is similar in both images. Fig. 2C shows a comparison between a side view of the completely closed pyloric sphincter from an X-ray fluoroscopy image⁷¹ and the corresponding representation in the stomach mesh. The light grey area in the X-ray image represents the barium-filled stomach in a healthy subject. The closure of the pyloric sphincter in this image creates a U-shaped gap between the stomach and the duodenum, a morphological feature which is also visible in the stomach surface mesh. Thus, the stomach mesh model integrates information from different imaging sources, such as MRI, CTG, cadaver, and endoscopy, generating a more comprehensive geometry representation of stomach morphology compared to previous computational models of human stomach in the literature.

2.3 Gastric contents

The model is to be used to simulate the gastric breakdown of a mixed solid–liquid gastric digesta, consisting of six spherical solid beads with a diameter of 12.7 mm, which are dispersed in a liquid digesta. The choice of using spherical beads as solid digesta in the model follows the experimental MRI work by Marciani et al.¹⁶ The current study is not aiming to replicate these experiments which are used to provide relevant conditions for investigating the mechanisms by which model foods, such as solid agar beads, are damaged and/or fractured inside the stomach. Two different material types are used, these being EP and EB. The bulk modulus and Poisson's ratio are kept constant across all solid beads (Table 1), with values corresponding to the elastic modulus of soft hydrogels, such as those made from gelatin⁷² or egg white.⁷³ This parameter choice ensures that the beads are sufficiently compliant to the prescribed wall contractions, without applying large reaction forces on the gastric walls which are assumed to be rigid structures. The low bulk modulus also enables observable deformations of the beads within a single peristaltic cycle, thus keeping computational costs and simulation runtime minimal.

Table 1 Material parameters used in modelling elastic-brittle and elastic-plastic deformation of solid beads along with the material parameters for the gastric liquid digesta

		Simulation cases
		Liquid- only	Elastic-plastic	Elastic-brittle
Liquid digesta	Volume	200 mL	200 mL	200 mL
	Density	1000 kg m^−3	1000 kg m^−3	1000 kg m^−3
	Viscosity	0.2 Pa s	0.2 Pa s	0.2 Pa s
Solid beads	Diameter	—	12.7 mm	12.7 mm
	Number	—	6	6
	Density	—	1005 kg m^−3	1005 kg m^−3
	Bulk modulus	—	21 000 Pa	21 000 Pa
	Poisson's ratio	—	0.45	0.45
	Yield strength	—	100, 250, and 500 Pa	—
	Threshold strain	—	—	0.01, 0.025, and 0.05

---

A sensitivity analysis is used to investigate the effect of changes to the strength of the beads on their mechanical response to gastric contractions. Three separate simulations are run for EP beads using different YS values of the solid bead material (Table 1), including a baseline YS of 100 Pa, as well as relatively stronger materials with YS of 250 Pa and 500 Pa. This broad range of YS values roughly corresponds to experimentally measured YS of spreadable food materials such as butter, margarine, and sour cream in their undigested form.^74,75 Similar to the EP beads, three separate simulations are run using different TS values for the EB material (Table 1), including a baseline TS of 0.025, as well as values of 0.01 and 0.05. These values are likely to be representative of heavily enzymatically and chemically degraded initially harder food materials.^76,77 The contact model parameters for the beads are dynamic friction coefficient of 0.25, which is relatively low to account for the effect of fluid films at contacts, and a coefficient of restitution of 0.6 to account for the dampening effect of the interstitial fluid.

The density of the liquid phase of the digesta is chosen to be 1000 kg m⁻³ equivalent to the density of water. The solid beads have a density of 1005 kg m⁻³, comparable to the density of soft solid foods such as eggs⁷⁸ and tofu.⁷⁹ The small density difference means that the beads remain submerged near the antral region, where flow dynamics is of primary interest. Lighter solids will instead rise to the top surface of the liquid digesta due to buoyancy and retrograde flows in the stomach and will not be subject to mechanical deformation.

The viscosity of the liquid phase of the digesta can influence the breakdown and emptying of solid digesta.^16,80 Previous in vivo studies typically compare pairs of test meals or beverages with differing viscosities.^16,81–84 In these studies, lower viscosity meals range from 0.042 to 0.08 Pa s, while higher viscosities range between 0.22 to 29.8 Pa s. In the current study, we aim to model the breakdown of solid beads within a liquid of representative intermediate viscosity. Viscosity variation resulting from addition of gastric secretions and chemical/enzymic dissolution of solid digesta occur over much longer timescales than those involved in this investigation, which is focussed on early-stage (60 s) digestion of large food fragments. On this basis, the liquid digesta is modelled as a Newtonian fluid with a representative constant viscosity of 0.2 Pa s.

2.4 Inlet–outlet boundary conditions

The model contains an open boundary at the proximal antrum, in addition to the closed boundary along the stomach wall. A Neumann boundary condition with zero-gradient pressure and velocity is applied at this upper surface of the digesta. This results in the surface moving freely in response to the ACWs, ensuring conservation of mass. Since ACWs in the stomach mesh are only active in antrum and pylorus, it is assumed that digesta flow in the remaining stomach would have minimal impact on flow and solid breakdown in the antrum. Therefore, the digesta in the rest of the stomach mesh is not included. In the configuration considered here, the duodenum is fully filled with digesta and assumed to have a zero-velocity closed boundary condition, resulting in a fairly strong back pressure at the pylorus. This leads to minimal duodenal flow and relatively low pyloric discharge rates through the open sphincter. The case with low duodenal pressure or partially evacuated duodenum, which will have much higher discharge rates and can be expected to materially affect the motion of the beads, will be examined in the next stage of this work.

2.5 Simulation configurations

Three sets of simulations are presented:

1. Liquid-only digesta,

2. Liquid content with EP solid beads, and

3. Liquid content with EB solid beads.

The SPH discretisation is used for all parts of the physical system. The solid and liquid digesta are discretised as SPH particles with an initial particle separation of 0.6 mm. The digesta within the antrum section of the stomach mesh (Fig. 1A) contains around 805 thousand particles, including six solid beads containing 5500 particles each. The digesta in the duodenum contains another 87 thousand particles. The overall stomach model involves approximately 100 mL of digesta, which adequately fills the pyloric and antral regions of the mesh.

Each simulation is run for a duration of 60 s, which is sufficient for at least one bead to be transported into the vicinity of the pyloric sphincter and undergo significant deformation and structure change.

2.6 Validation of the SPH model

The current SPH model for gastric breakdown of soft solid foods integrates multiple subsystems, including fluid-surface interactions at the deforming gastric walls and the resulting peristaltic flow, fluid-solid interactions, and the compression and fracturing of solid beads. The CSIRO SPH code^41,42 used in this study has been widely validated for simulating complex fluid and solid behaviours across each of these individual sub-systems. For instance, the SPH implementation has previously demonstrated high accuracy in modelling free-surface water behaviour and interactions with complex geometries, outperforming traditional grid-based methods in terms of correlation with experimental data.⁸⁵ Moreover, SPH models of interaction between flow and moving geometries in industrial mixers have excellent experimental alignment.⁸⁶ Recently, Liu et al.⁴⁰ have applied the same SPH code to model peristaltic flow in a tube, and the model predictions are consistent with analytical solutions for net flow corresponding to a range of occlusion ratios. The predicted velocity fields are also comparable to those produced by grid-based methods (finite volume method). These studies collectively validate the capability of the current SPH model to accurately simulate fluid flows involving complex gastric wall movements.

It is also essential to evaluate the accuracy of the SPH implementation in resolving fluid flow around solid objects and accurately predicting the flow dynamics, drag and solid motion. Previously, SPH have been used to simulate the motion of pellets in a large impeller-stirred mixing tank, with model predicted particulate motion showing strong agreement with experimental results.⁸⁷ To further validate SPH for the specific context of the current study, a simplified version of the stomach model is presented in the appendix. The model consists of a steady parabolic flow established through a cylinder (radius of R), which transports a spherical solid bead (radius of r) along the cylinder axis. The cylindrical domain's volume and length, along with other parameters such as flow speed, fluid density, and viscosity are adjusted to match conditions in the primary stomach simulations. Drift velocities of the bead are calculated for various ratios of r/R and compared against analytical solutions. The model predictions closely match with corresponding analytical results. Additionally, numerical convergence is assessed by varying the initial particle separation, confirming the robustness of the predictions across different spatial resolutions.

Solid foods in the antrum are subjected to multiaxial compressive loading, and the resulting compression and deformation require additional validation. Under low compressive stresses, both EP and EB beads deform elastically. Previously, Das and Cleary⁸⁸ compared stress wave attributes predicted by the SPH code to those from a commercial FE model for uniaxial compression of an elastic object. The results not only showed strong agreement with both analytical and FE solutions, but SPH also produced a smoother response than FE in early loading stages. Das and Cleary⁵⁴ further showed that SPH solutions for Von Mises stress distributions under uniaxial loading have good convergence with increasing particle resolution/decreasing initial particle separation. In another study, Das and Cleary⁸⁹ validated the fracture modelling capability of SPH by comparing its predictions for dam failure under earthquake loading with experimental scale models and FE simulations. The SPH model accurately predicted crack initiation and propagation, consistent with both experimentally observed and FE predicted failure patterns. It also avoided the need for predefined crack zones or costly remeshing required in FE models.

These sub-system validations demonstrate that the current SPH model is well-equipped to accurately simulate the complex fluid dynamics and mechanical breakdown of soft solid foods in a deforming stomach.

3 Flow behaviour of the stomach for a liquid-only digesta

First, we establish the nature of the flow of liquid-only content in the stomach subject to the stomach wall deformations and pylorus opening and closing. Fig. 3 and 4 show the flow behaviour in the stomach for this liquid-only content. Fig. 3A shows the initial setup of the liquid digesta inside the stomach mesh at 0 s. The digesta is shown as transparent light blue. The pyloric sphincter is initially closed (at 0 s), but begins to open from 2 s. Fig. 3B shows a cross-sectional view of the flow within a slice through the centre of the pyloric sphincter and terminal antral region at 3 s. The liquid is coloured by fluid speed using the spectrum given at the bottom of the figure. By 3 s, a parabolic velocity profile is established through the pyloric sphincter with flow orientated towards the duodenum. The speed of the flow varies between 10 mm s⁻¹ near the sphincter wall up to 20 mm s⁻¹ at its centre.

Fig. 3 Peristaltic mixing of liquid-only digesta: (A) initial simulation setup of a liquid-only digesta inside a deformed stomach. (B) 3D visualisation of the flow behaviour within a 0.5 mm thick cross-sectional plane passing through the centre of the pyloric sphincter, colour-coded using a rainbow spectrum. (C–F) Velocity of the liquid digesta within the 0.5 mm thick plane at various time points in the simulation is represented by constant-length vectors and coloured with the same rainbow spectrum.

Fig. 4 Retropulsive flow induced by TAC: this figure shows the velocities of the liquid digesta within the terminal antrum region, in a 0.5 mm thick cross-sectional plane passing through the centre of the pyloric sphincter, from 10.00 s to 12.00 s of simulated time (A–F). The direction of motion of the liquid digesta is indicated by constant-length vectors and flow speed by the colour.

Fig. 3C–F show vector representations of the velocities of the fluid located within the same cross-sectional plane depicted in Fig. 3B, for the period of the ACW contraction between 6 s and 10 s, during which the pyloric sphincter closes but the TAC has not yet initiated. At 6 s (Fig. 3C), the backflow through the pylorus decreases to relatively low speeds below 5 mm s⁻¹ due to pressure equalisation between the stomach and the duodenum. The liquid digesta situated ahead of the ACW flows away from the contracting antral wall, while the digesta behind the ACW flows towards the expanding antral wall. Between 7 s and 9.0 s (Fig. 3D and E), the pyloric sphincter closes and the ACW continues towards it. The backward flow intensifies, reaching speeds of up to 12 mm s⁻¹, which is shown by red arrows. Backward flow exceeding 8 mm s⁻¹ is also observed near the pylorus due to the fast final closure movements of the pyloric sphincter. At 9 s (Fig. 3E), the pylorus becomes completely closed. In the next 0.5 s, the ACW occlusion moves further towards the closed pylorus with increasingly higher occlusion amplitude, pulse width and speed. The width of the occluded antral region decreases rapidly, signalling the commencement of the TAC (at 9.5 s). The terminal antral area appears to contract rapidly similar to a segmental wave, compared to the slower peristaltic waves seen in Fig. 3B–E. Due to this rapid contraction, the backward peristaltic flow is further strengthened, as evident from the increasing number of red vector arrows displayed between 9 and 9.5 s (Fig. 3E and F). Recirculating currents are induced in both the top and bottom regions of the antrum, which could contribute strongly to mixing in these areas.

Fig. 4 shows the effect of the higher speed wall contractions during the TAC on the flow of liquid digesta in the antral region. Between 10 and 11 s (Fig. 4A–C), all of the stomach wall proximal to the pyloric sphincter contracts rapidly, leading to fast reduction in the volume available for the liquid. As a result, retrograde fluid speeds of up to 70 mm s⁻¹ are induced, which are substantially higher than the speeds associated with the ACWs or the pylorus opening and closing. Simultaneously, the pair of eddies move further left towards the pylorus. These effects are similar to predictions from previous computational fluid mechanics studies, where retropulsive jets and recirculating currents have also been reported.^25,28,35 The TAC reaches its peak strength at around 11 s (Fig. 4C), after which the contracted terminal antral region rapidly expands back to its resting configuration, leading to a reversal in the direction of the flow. A forward jet towards the pylorus is then observed in the antral region between 11.25 and 11.75 s (Fig. 4D and E), reaching speeds of up to 25 mm s⁻¹. Subsequently, by 12 s (Fig. 4F), the antral region enters a state of relative quiescence. The anticlockwise eddy at the bottom of this region expands radially and progresses towards its centre, becoming the dominant feature of the antral flow. These flow structures are likely to play a crucial role in gastric mixing, although current in vivo imaging techniques do not yet have the capability to visualise such flow at this level of detail.

4 Deformation of elastic-plastic solid beads

Fig. 5 shows the initial arrangement of six spherical EP solid beads distributed in the liquid digesta, which is the same for all three YS cases. These beads are numbered from 1 to 6, as indicated. Bead 1 is positioned at the terminal antrum near the pylorus, allowing us to study the effect of TACs on solid deformation and breakdown within the 60 s simulation duration. The remaining beads are positioned in the antrum at different distances from the pylorus, to facilitate understanding of the flow behaviours induced by ACWs at different locations within the antrum. All SPH particles in the simulations, including the solid beads, initially have zero velocity.

Fig. 5 Initial arrangement of six initially spherical solid beads and liquid digesta within the stomach geometry. The solid beads are numbered from one to six and coloured differently.

4.1 Trajectory and deformation of elastic-plastic solids

The simulation predictions of deformation and fragmentation of the EP beads are shown in Fig. 6 for the baseline YS of 100 Pa. Fig. 6A shows 3D visualisations of the deformed beads at the conclusion of each of the first three 20 s duration ACW cycles. Several large-scale deformations are prominent, which are discussed in detail in the next subsections. For instance, in the first 20 s bead 1 (purple) undergoes significant extensional deformation into a long cylindrical shape during the first ACW cycle. Between 20 s and 60 s, beads 4 (dark blue) and 6 (yellow) collide and change shape considerably. Fig. 6B shows the trajectory of each bead over one (20 s), two (40 s) and three (60 s) ACW cycles. Bead 1 follows a closed loop, returning to its initial position after each ACW cycle, as fluid flows away from and then back into the terminal antrum. Similar nearly horizontal loop-like trajectories are observed for beads 2 and 5 in the antral-pyloric region. Loop-like trajectories of beads 1, 2 and 5 align with the periodic reversal of flow direction observed in the liquid-digesta-only simulation (Fig. 4). These results are consistent with model predictions of Seo et al.³⁵ and Lee et al.,³⁶ where similar cyclic trajectories are observed for a dynamic pill in the terminal antrum. In contrast, beads 3, 4, and 6, which are located further away from the pylorus, gradually descend towards the bottom of the stomach with each ACW cycle. The effect of the retrograde and then prograde jets from and then back into the terminal antrum are weaker in the central regions of the stomach, so these beads do not exhibit the same oscillatory behaviour as the other beads.

Fig. 6 The elastic-plastic deformation of solid beads by the stomach and liquid content is shown over 60 s: (A) deformation of the beads in the antral region, driven by peristalsis and TAC. Small fragments detached from the beads are shown in opaque white. (B) Visualisation of the trajectory of the centroids of each bead up to the time indicated.

4.2 Pinching effects of the pyloric sphincter

Fig. 7 shows close-up 3D visualisations of the deformation of bead 1 (for a baseline YS of 100 Pa) induced by closing of the pyloric sphincter. The beads are coloured by plastic strain and displayed from two perspectives:

(i) from the side of the pyloric sphincter (shown on the left of each pair of images at each time), and

(ii) perpendicular to the former (shown on the right of each pair), offering a view along the direction of the pylorus.

The first view corresponds to the one used in Fig. 3 and 4 for the liquid-only digesta. Similarly, it also includes vector representations of the fluid flow in the antral-pyloric region surrounding bead 1, but colour-coded in shades of magenta from light (0 mm s⁻¹) to dark (3.5 mm s⁻¹). At 3 s (Fig. 7A), the left image shows the flow from the stomach to the duodenum driven by the pressure differential between the two regions, when the pyloric sphincter is partially open. The speed of this flow through the sphincter exceeds 3.5 mm s⁻¹, as shown by the darker arrows near the sphincter. By 4 s (Fig. 7B), the flow propels bead 1 into the open sphincter, where it becomes lodged against the sphincter wall due to the current sphincter diameter of 10 mm being smaller than the bead diameter of 12.7 mm. The jammed bead obstructs fluid flow through the pyloric sphincter, while the fluid pressure acting on the bead balances the contact forces from the sphincter wall, resulting in no net bead motion during this phase of the pylorus cycle.

Fig. 7 Visualisation of gastric fluid dynamics and deformation of an elastic-plastic solid bead due to closure of pyloric sphincter: two images are given at each time from 3.0 s to 8.0 s (A–F), showing two different perspectives. The first image at each time (on the left) shows a side view of the terminal antrum with vectors indicating the fluid velocities (coloured by fluid speed), while the second image (on the right) shows a frontal view along the direction of the pyloric sphincter. Bead colour indicates plastic strain which is generated by wall and fluid interaction during pyloric closure.

After 5 s (Fig. 7C), the pyloric sphincter gradually begins to close, radially compressing the left end of the trapped bead. The fluid pressure and wall friction prevent the bead from escaping from the pylorus, allowing the radial compression to load a thin layer of the left end of the bead that is deformed plastically with plastic strain increasing to around 0.5 (blue colours in Fig. 7E). As the net wall force acting along the direction of the pylorus exceeds the fluid pressure gradient, the closing sphincter becomes dominant, slowly pushing the bead back towards the antrum, and by 8 s, it has become dislodged from the nearly closed sphincter (Fig. 7F). This pyloric interaction results in only a small area of high plastic deformation with peak plastic strains around 0.7 at the left-end of the bead which was the main contact area between the bead and the sphincter walls. However, several tiny fragments, with sizes of 0.6 mm and high plastic strains over 1.0, are generated in the contact area (Fig. 7F). The areas of the bead surface that are in contact with the sphincter wall are restricted in their ability to move tangentially by the wall friction, while the core of the bead is forced laterally out from the sphincter. This leaves a thin layer of the retarded material that becomes pinched off, generating very small fragments. After 8 s, the bead 1 is far from its initial spherical shape due to the strong deformation it has experienced and minor pinching fragmentation occurring at the sphincter contacts.

4.3 Large scale elastoplastic deformation in the antrum

Fig. 8 shows a 3D visualisation of the effect of the TAC on bead 1 and the surrounding fluid flow, immediately after the pinching that occurs at the pyloric sphincter (Fig. 7). At 8.75 s (Fig. 8A) the occlusion in the terminal antral wall progresses towards bead 1 with increasing contraction amplitude, speed, and pulse width (as also shown in Fig. 1B–D). Concurrently, in line with observations from the liquid-only simulation (Fig. 3C), the TAC results in a strong retropulsive flow with speeds exceeding 8 mm s⁻¹ near the right end of the bead 1 (indicated by darker arrows in Fig. 8A). The force from this flow results in axial loading opposing the direction of TAC (Fig. 8B–D). Simultaneously, the contracting antral wall compresses the bead radially. Due to this combined axial and radial loading, the right half of the bead is squeezed radially inwards and is pushed axially through the narrow occlusion in the antral wall. This is similar to the manufacturing process of extrusion, where a material is pushed through a small opening of specific cross section.

Fig. 8 Visualisation of fluid dynamics and deformation/fragmentation of an elastic-plastic solid bead (yield strength of 100 Pa) due to TAC between 8.75 s and 10.00 s (A–F).The direction of the liquid digesta motion is represented by constant-length vectors, with speed shown by their colour. The colour of bead 1 indicates plastic strain.

Discussion of the complex behaviour of the bead from this point is aided by the schematic in Fig. 9. Between 9.0 and 9.5 s (Fig. 8B–D and Fig. 9A), the overall length of bead 1 increases by 40% due to the observed TAC induced extrusion of the right half of the bead. Similar backward extrusion of soft solids due to TAC has also been reported in Skamniotis et al.'s computational study³⁷ for viscoplastic materials. The contracting terminal antral wall directly applies a radial compressive force on the right half of the bead (Fig. 9A). At 9.75 s (Fig. 8E), this increasing radial and axial loading results in an axial gradient of plastic strain from the left end to the right end of the elongated bead 1. Additionally, a narrow tail-like region (Fig. 8E) is created at the left-end of the bead which is still trapped within the TAC occlusion. At 9.75 s the bead length has increased by 60% of its original diameter. Over the next 0.25 s, bead 1 is completely squeezed out of the contracted terminal antrum with a shape that is now nearly cylindrical. However, by 10.0 s (Fig. 8F), the elongation of bead 1 reduces to 44%, which is lower than the peak elongation of 60% at 9.75 s. This reduction in length can be attributed to axial compression by the strong retrograde flow at around 10.0 s, indicated by the darker velocity arrows in Fig. 8F (blue arrow in Fig. 9B). As a result, the peak plastic strain at the left-end of the bead 1 increases up to 2.1 at 10.0 s (Fig. 8F) from 1.2 at 9.5 s (Fig. 8E). This behaviour contrasts with the previous phase of bead deformation (Fig. 9A), where fluid flow led to extrusion of the bead through the contracting antral wall. These results suggest that the fluid pressure of the TAC induced retropulsive flow can strongly influence the shape and rate of deformation of different parts of an EP solid in a number of different ways (axial compression or elongation), depending on the phases of the deformation and the strength and direction of the fluid flow.

Fig. 9 Schematic representation of TAC-induced deformation of elastic-plastic bead 1 during two distinct phases: (A) 9.0–9.5 s and (B) 9.5–10.0 s. Black arrows indicate the direction of radial compressive forces exerted by the contracting terminal antral wall, blue arrows represent the hydrodynamic flow direction, and pink broken arrows show the deformation of the solid bead resulting from a combination of wall compression and flow pressure.

Fig. 10 shows the total surface area and the number of fragments generated from the six beads as a function of time. For the baseline YS of 100 Pa, total surface area increases sharply (from 39 cm² to 41 cm²) during the first TAC event between 9 s and 10 s. This is driven by the extrusion of bead 1 into a long cylindrical shape (Fig. 8 and 9). This process also results in the generation of several small fragments with sizes varying between 0.6 mm and 2 mm. This fragmentation primarily occurs at the pylorus-facing left-end of the bead 1, which has already undergone plastic deformation and minor fragmentation due to pinching during pyloric sphincter closure prior to the onset of TAC (Fig. 7F). Between 9.0 and 9.75s, a radially inward flow is generated in the narrow gap between the bead 1 and the pyloric sphincter (Fig. 8B–E), causing accumulation of additional plastic strain that results in this observed fragmentation.

Fig. 10 Quantification of temporal variation in (A) fragment generation and (B) total surface area of the six elastic-plastic solid beads, for three bead yield stresses.

4.4 Collisions between elastic-plastic solid beads

Soft food can deform not only due to fluid pressure and contact forces from the deforming stomach wall, but also through collisional interactions between different food masses. This can therefore occur in areas of the stomach other than the terminal antrum and pylorus. Fig. 11 shows a sequence of images visualising the inelastic collision of beads 4 and 6 near the bottom antral wall for the baseline YS of 100 Pa. As observed in the liquid-only simulation (Fig. 3), ACWs induce a flow towards the expanding rear end of the contracted zone of the stomach wall. This flow carries the beads downward towards the bottom of the stomach, where they collide at around 15 s (Fig. 11A). The flow pressure acting downward and the upward reaction force from the stomach wall press the beads together leading to substantial deformation. There is no adhesion in this model, consequently the beads are not bonded together and can, in principle, move independently if the fluid mechanics so dictate. By 20 s, the beads begin to deform against each other due to compression from the combination of the flow pressure and the antral wall reaction force (Fig. 11B). As each new ACW approaches (Fig. 11C and G), the increased contact force from the antral wall lifts the beads up, intensifying the compression between them, resulting in significant further deformation, for instance with peak plastic strain of 1.8 by 30 s (Fig. 11C). Soon after (Fig. 11D and H), the two beads are both moved towards the left (in the direction of the pylorus), carried by the traveling occlusion in the antral wall, akin to an escalator. At the same time, a backward flow towards the rear end of the occlusion starts developing. Over the next 5 s, this flow pushes the beads back down from the occluded region of the gastric wall towards the expanding rear end of the occlusion (Fig. 11E and I). Due to this interplay between the ACW occlusions and the backward flow, the two beads continue to be pressed together with both their deformation levels increasing. At 60 s (Fig. 11I), the peak plastic strain in the contact area between beads 4 and 6 has increased beyond 3.0. This process has also resulted in the generation of multiple small fragments that are separated from the two beads. It is slower and more continuous compared to the TAC induced rapid deformation of bead 1 (Fig. 8 and 9), which occurs over a couple of seconds. As a consequence of this deformation, the total surface area of the beads gradually rises from 41 cm² at 10 s to 53 cm² by 60 s (Fig. 10B). In the same time period, 45 small fragments are generated (Fig. 10A), which have a minimal contribution in increasing the total surface area.

Fig. 11 Visualisation of plastic deformation in two colliding elastic-plastic solid beads in the distal stomach between 15 s and 60 s (A–I). The motion of liquid digesta is represented by constant-length vectors, colour-coded based on their speed. The colour of the agar beads indicates the accumulated plastic strain.

4.5 Effect of yield strength on elastic-plastic deformation

Fig. 12 shows a 3D visualisation of fluid dynamics and TAC induced deformation of EP beads with a YS of 500 Pa. Comparison of Fig. 12 and 8 (baseline YS value of 100 Pa) shows noticeable differences in the two phases (Fig. 9) of the TAC induced deformation process of bead 1. Firstly, the bead 1 with a YS of 100 Pa shows visibly significant deformation and plastic strain before 9.0 s (Fig. 8B), unlike the stronger bead (Fig. 12B) whose first deformation is delayed. Similarly, by 9.75 s, the stronger bead undergoes substantially less radial squashing and elongation (only 38%) compared to the bead with lower YS (60% elongation). While the elongation increases during the first phase till 9.75 s, it is partially reversed in the second phase due to axial compression by strong retrograde flow. In contrast, plastic strain accumulates irreversibly during both phases. Due to the lower YS of the softest bead (100 Pa), it undergoes greater plastic deformation in both phases, resulting in a higher plastic strain. The 100 Pa bead 1 has a final peak plastic strain of 2.1 at its left end (Fig. 8F), compared to 0.8 for the 500 Pa bead 1 (Fig. 12F). A comparison of these figures further indicates that the second phase of deformation—compression of the ‘tail’-like region near the pylorus—is considerably more pronounced for the softer bead. This leads to a larger decrease in elongation for the softer bead 1. As a result, both the 500 Pa (Fig. 12) and 100 Pa YS (Fig. 8) beads ultimately have comparable final elongation (47% for YS = 500 Pa and 44% for YS = 100 Pa), despite having different shapes in the earlier and intermediate stages of deformation (reflecting the reasonably different deformation histories) and different final levels of plastic strains. Therefore, unlike plastic strain, final elongation after TAC does not linearly correlate with the YS of EP beads, as it responds differently to the two phases of deformation.

Fig. 12 Visualisation of fluid dynamics and TAC induced deformation and fracture of an elastic-plastic solid bead with a yield strength of 500 Pa between 8.75 s and 10.00 s (A–F).

At 10 s, the stronger 500 Pa beads tend to generate fewer sub-mm fragments compared to the YS of 100 Pa (Fig. 10A). Additionally, Fig. 10B shows that the sharp rise in the total surface area of bead 1 during the first TAC event (9 s–10 s) is comparable across different YS values. This observation is consistent with the similar final elongation of bead 1 across different YS values. These results suggest that the TAC can deform soft solids with a broad range of material strengths, which would appear to be biologically advantageous.

5 Deformation and fracture of elastic-brittle solid beads

Six EB spherical beads are arranged in the liquid digesta using the same configuration as for the EP beads illustrated in Fig. 5. Each bead is again numbered from one to six for reference.

5.1 Trajectory and deformation of elastic-brittle solid beads

Fig. 13A shows 3D visualisations of the EB beads at the conclusions of three consecutive ACW cycles (time period of 20 s) for the baseline model with a TS of 0.025. The same contraction patterns are present in each of the 3D visualisations, with the shape of the boundaries appearing unchanged. Bead 1 undergoes minor deformation but significant fracturing from the compression that occurs during the first TAC event (unlike the EP case, in which significant plastic deformation occurred). The bead is elongated by only 16% from its original diameter, while 7% of its volume is fragmented. The fragment products are generally smaller than 1 mm, creating moderate density clouds of fine solids (shown in golden yellow in Fig. 13A). These sub-mm sized fragments move with the fluid flow because of their very low particle Reynolds number, as a consequence of their small size.

Fig. 13 Elastic-brittle deformation of solid beads: (A) deformation and fracture of the solid beads in the antral region, driven by peristalsis and TAC. Fragments detached from the beads are shown in golden yellow. (B) Visualisation of the trajectories followed by the centroids of each solid bead over three different time intervals.

Fig. 13B shows the cumulative trajectories of each EB bead across three successive ACW cycles. The trajectories of the EB beads are very similar to those of the EP beads (Fig. 6B) up to the points of significant deformation and fracture, which should be expected as these motions are predominantly driven by the same controlling fluid flow. Beads 1, 2, and 5, which are located near the pylorus, are observed to move in loop-like trajectories. Downward trajectories are observed for beads 3, 4, and 6, leading to a collision between beads 4 and 6. However, unlike the EP beads (Fig. 6B), some of which are pressed together and slowly compressed by the flow, the EB beads simply bounce off each other after the collision and then move apart on divergent trajectories.

5.2 Elastic-brittle deformation and fragmentation in terminal antrum

Fig. 14 shows the mechanisms leading to the deformation and fragmentation of the EB bead 1 during the first ACW cycle. Fig. 15A further shows the temporal variation in average damage experienced by bead 1, for the baseline TS of 0.025. At 9 s (Fig. 14A), its left side shows low levels of damage (areas coloured blue) from earlier interaction with the closing pyloric sphincter, similar to that causing the pinching deformation behaviour observed for EP beads (in Fig. 8A and B). Over the next 0.25 s (Fig. 14B), the right end of the bead is compressed radially by the contracting terminal antral wall. As a result, the bead is elongated by approximately 5% from its initial spherical diameter. At 9.5 s (Fig. 14C), the length of the EB bead 1 has increased by 15%, which is lower than the 22% elongation found for the stronger EP bead (in Fig. 12D) and the 40% value for the softer EP bead (Fig. 8D) at the same timepoint. Moreover, unlike the EP beads, a radial crack is generated in the middle of deformed EB bead 1, where the compressive strain is maximal and exceeds the TS. Low-level damage of 0.3 can be observed in the elongated right-end of the bead which is positioned inside the occlusion.

Fig. 14 Visualisation of TAC induced fluid dynamics and deformation of a nearly intact elastic-brittle solid bead between 9.00 s and 11.50 s (A–H). The motion of the liquid digesta is represented by constant-length vectors with colouring according to the flow speed. The colour of the solid beads shows the damage accumulated during the TAC.

Fig. 15 The plots quantify the temporal variation in (A) cumulative average damage in the elastic-brittle solid bead 1, subjected to TAC deformation, (B) total number of fragments generated from the six elastic-brittle beads and (C) total surface area of the elastic-brittle beads for three threshold strains.

At 9.75 s (Fig. 14D), the EB bead 1 is subjected to a hydrodynamic pressure gradient which pushes it away from the pyloric sphincter. The left half of the bead is now subject to loading from the occlusion. At the same time, the right half of the bead, which was radially squeezed to fit the occlusion, is released with a relative increase in cross-sectional area from this elastic unloading. The radial compression of the bead between 9.5 s and 9.75 s leads to additional damage, which accumulates in the right half of the released bead. At 9.75 s, a peak damage level of 1.0 (coloured by red, indicating full damage, allowing material on either side to move apart) is observed at the radial crack. The increase in damage is also accompanied by an increase in width of the crack. Overall, these EB behaviours contrast with that of the EP beads (in Fig. 8E and 12E), where the right half of the EP bead 1 is squeezed out in an elongated shape without any recovery of the cross-sectional area and crack propagation.

By 10 s, EB bead 1 has been completely expelled from the terminal antrum (Fig. 14E). The total elongation of the bead has decreased from 26% at 9.75 s to 18% at 10 s. Visual comparison between Fig. 14D and E indicates that this length reduction is localised to the pylorus-facing left end of the bead, which is squeezed out of the occluded terminal antrum without any fracturing. This reduction in length is a consequence of the left end of the bead unloading elastically, as the radial compressive load from antral wall contraction is removed by 10 s. A similar relative decrease in elongation of the left half of the bead is also observed in the case of the EP beads (Fig. 8 and 12) over the same time interval (9.75 s–10.0 s). Fig. 14E also shows generation of several fine fragments (with sizes of 0.6 mm, equivalent to SPH spatial resolution) from the periphery of the centre-plane of the bead, orthogonal to the axis of the pylorus. These sub-mm sized fragments result from the frictional interaction of the bead surface with the occluded terminal antral wall between 9.5 s and 10 s. During this time interval, the TAC wavefront moves rapidly towards the pylorus in a direction opposite to the bead motion. Consequently, the surface of the bead in contact with the moving antral wall experiences resistance and strong shear stress, as the core of the bead is squeezed through the TAC occlusion under the pressure of the retropulsive jets. This leaves a thin peripheral layer of retarded material that loses mechanical continuity with the core of the deformed bead and is left behind in the form of small fragments (coloured as red in Fig. 14E). At 10.0 s, the fragments that have separated from bead 1 account for 0.8% of its initial volume. The average damage increases sharply from nearly zero at 9 s to 0.52 at 10s (Fig. 15A), indicating widespread structural damage to the bead. Over the next 1.5 s (Fig. 14F–H), the maximally damaged part of the bead surface is exposed to the strong TAC induced eddies in the antrum (previously discussed in section 3 and Fig. 4). Consequently, additional SPH particles in the surface regions of the bead reach the maximum damage level leading to their fragmentation and separation from the main body of the bead. The pair of eddies then disperse these clouds of sub-mm sized fragments away from the bead. At the end of the first TAC (Fig. 14H), bead 1 resembles a capsule shape with many void areas in the surface due to partial fragmentation and is surrounded by a ring-like cloud of fragmented sub-mm solids. This final fragment cloud, containing a total of 235 fragments (Fig. 15B), accounts for 7% of the initial bead volume.

Fig. 16 shows the effect of the second TAC event on the previously deformed, damaged and partially fragmented EB bead 1. At 28.25 s, bead 1 again becomes trapped within the occlusion of an approaching TAC (Fig. 16A). Over the next 0.5 s (Fig. 16B and C), the occluding terminal antral wall contracts further and radially compresses the right half of the deformed bead, leading to an axial elongation of the bead into the shape of an hourglass. Simultaneously, the retropulsive flow starts to develop, as the volume available for the liquid decreases rapidly. The elongated bead is pushed further away from the pylorus towards the occlusion due to the resulting hydrodynamic pressure gradient loading it axially (Fig. 16D). The centre of the bead, which is maximally compressed and partially fractured from the previous TAC, undergoes complete fracturing, leading to the separation of the deformed bead into two primary parts at 29.25 s (Fig. 16E). At 29.50 s (Fig. 16F), the terminal antral radially compresses the pylorus-facing left half of the bead 1, which is then squeezed through the occlusion under the pressure of the retropulsive flow. By 29.75 s (Fig. 16G), the left half has been completely released from the occlusion and moves further away, then collides with the separated right half of the bead at 30 s (Fig. 16H). However, a considerable number of small and highly damaged fragments are left behind inside the occluded terminal antrum, as the retrograde flow cannot sufficiently quickly carry all these small particles out of the occlusion. Overall, the deformation and fragmentation of bead 1 during the second TAC event occurs over a longer timeframe (∼2 s), compared to the first TAC event shown in Fig. 14 (which was around 1 s). The increment in average damage level is also minimal (Fig. 15A). This is a consequence of the more elongated and damaged nature of the bead 1 after the first TAC, compared to its initial intact spherical form. By 30 s, the volume fraction of the separated sub-mm sized fragments has increased to 10% of the initial bead 1 volume, compared to 0.8% at the same periodic timepoint during the first TAC event (10 s). The second TAC leads to the generation of an additional 124 fragments (Fig. 15B). This is followed by the generation of 350 more fragments, when further breakdown of bead 1 occurs during the third TAC (Fig. 13A and 17F), resulting from the same mechanistic behaviours as those described for the first two TACs but applying to increasingly smaller, irregularly shaped larger fragments. With each TAC, the average damage level increases progressively but in smaller increments (Fig. 15A).

Fig. 16 Visualisation of TAC induced fluid dynamics and fragmentation of a previously damaged elastic-brittle solid bead which has undergone partial breakdown from a previous TAC between 28.25 s and 30.00 s (A–H).

Fig. 17 Visualisation of TAC induced fluid dynamics, damage and fragmentation of a nearly intact elastic-brittle bead during three consecutive TAC events, corresponding to three threshold strains of the bead material: 0.05 (A–C), 0.025 (D–F), and 0.01 (G–I).

Fig. 15C shows the total surface area of the six EB beads with a baseline TS of 0.025, plotted as a function of time, which contrasts with the corresponding plot for EP beads with baseline YS of 100 Pa (Fig. 10B). For the EB beads, the total surface area increases from 38 cm² to 47 cm², which is substantially higher than the corresponding value (41 cm²) for baseline EP beads at 10 s (Fig. 10B). This large increase is driven by both deformation of the EB bead 1 and the generation of the 235 fragments during the initial TAC event (Fig. 15B). In total, more than 2000 fragments are generated in the three ACW cycles (60 s). In contrast, EP bead 1 experiences substantial elongation and reduction in cross-sectional area but limited fragmentation (Fig. 10A). Consequently, after 60 s, the total surface area for baseline EP beads reaches only 53 cm², which is considerably lower than the 121 cm² observed for baseline EB beads (Fig. 15C). This demonstrates that the rheological nature of the food material has an important effect on the deformation response under compression.

Another key difference between EP and EB bead 1 is the sequence of their deformation during successive TAC events. Deformation of the EP bead 1 occurs only during the first TAC event (Fig. 6A). This is a consequence of the elongated cylindrical shape of the EP bead 1 after the first TAC, which closely conforms to the stomach wall shape at the end of the TAC. As a result, there can only be weak further compression in subsequent TAC cycles. In contrast, the EB bead 1 forms a capsule-like shape without a large change in cross-sectional area, so it can easily be trapped and compressed by subsequent TAC waves. Therefore, during the onset of the second TAC event, the EB bead 1 is much more likely to be trapped inside the narrow occluded region of the terminal antrum with peak compressive stresses that are much higher. Moreover, due to the streamlined shape, the EP bead 1 moves faster under similar flow conditions, reducing the probability of it being trapped inside the occluded terminal antrum in future antral cycles.

5.3 Effects of threshold strain on elastic-brittle deformation

Fig. 17A–C show the deformation and damage of bead 1 at the midway points of each TAC event for the highest TS of 0.05, which is the strongest material variant. These results can be visually compared with Fig. 17D–F, which show the fragmentation of bead 1 for the baseline TS of 0.025. At 10 s, the stronger EB bead 1 (Fig. 17A) is deformed by the initial TAC to a capsule-like shape similar to that found in the baseline simulation (Fig. 17D). However, the damage levels on the surface of bead 1 in Fig. 17A are visibly lower compared to Fig. 17D. Moreover, radial crack lines are absent in the stronger bead 1. These results are consistent with the time variation of average damage in Fig. 15A, where the average damage level at 10 s is substantially smaller for TS of 0.05 than for 0.025. In the next TAC cycle at 30 s (Fig. 17B), the damage profile on the surface of the stronger bead 1 remains nearly unchanged. However, while it is squeezed out of the occlusion during the third TAC cycle at 50 s (Fig. 17C), a damaged area develops in the middle of the bead with a rough surface. At this time, only 2% of the bead volume is separated in form of fine fragments. This final form of the stronger bead 1 (TS of 0.05), after the third TAC, resembles the shape and fragmentation (0.8%) of bead 1 during the first TAC (10 s) for the baseline case with TS of 0.025 (Fig. 17D). Similarly, the total surface area of the stronger beads at 50 s is comparable to the same for the baseline case after 11 s (Fig. 15C). Therefore, higher TS reduces the degree to which mechanical breakdown can be achieved. This highlights the important role of chemical and enzymatic degradation of food materials in softening them, enabling faster and more intensive mechanical breakdown.

The deformation and fracture of bead 1 with the lowest TS of 0.01 is shown in Fig. 17G–I, corresponding to three consecutive TAC events (10, 30 and 50 s). This weakest bead 1 fractures rapidly and liberates approximately 400 small fragments during each TAC event (Fig. 15B). During the first TAC (10 s), a crack propagates through the middle of the bead 1, and the surface adjacent to the crack experiences high levels of damage (coloured as yellow to red in Fig. 17G) exceeding 0.9. The area of high damage in Fig. 17G spreads over a larger fraction of the bead surface compared to the baseline case (TS of 0.025) at the same timepoint (Fig. 17D). The first TAC event also leads to separation of several sub-mm sized fragments from this weaker bead 1. At 10 s, the 1.0% volume fraction of the fragments generated is comparable to the 0.8% value for the baseline case. The differences in fragmentation behaviour of bead 1 between TS of 0.01 and 0.025 become more prominent during the second TAC (30 s). At this timepoint (Fig. 17H), a large fraction of the weaker bead 1 is squeezed out of the occlusion in the form of mm and sub-mm sized fragments (all red since they are at the maximum level of damage). Moreover, a substantial volume of larger partially damaged fragments (shown in green) is still entrapped within the occlusion of the second TAC. In contrast, for the baseline case, only a few fragments are located within the occluded antral region (Fig. 17B). The largest chunk of the fragmented weaker bead 1 (Fig. 17H) accounts for only 50% of its initial volume, compared to 90% for the baseline case (Fig. 17B). These results are consistent with the large differences in fragment counts between TS of 0.025 and 0.01 at 30 s, as shown in Fig. 15B. For the lowest TS, many fragments are also generated between TAC events due to accumulation of damage by the forces exerted by the ACWs and the resulting fluid flow (Fig. 15A and B). Because of this continuous fracturing, the total surface area of the bead increases substantially between TAC events (Fig. 15C). During the third TAC (50 s), average damage for bead 1 approaches the maximal value of 1.0 (Fig. 17I and 15A). The largest fragment of bead 1 consists of less than 1% of the initial bead volume, signifying its disintegration. The sub-mm sized fragments collide with the surface of bead 2 located nearby, leading to its accumulation of a low-level damage below 0.5 (shown in green blue in Fig. 17I). These results highlight the role of particle-to-particle collisions in the breakdown of very soft solid digesta. Overall, the total surface area and number of fragments generated after the three ACW cycles positively correlate with the TS of the EB beads (Fig. 15B and C).

Some of the fine fragments generated from the weakest bead 1 (TS = 0.01) are transported into the duodenum through the periodically opening pyloric sphincter. At 60 s, the majority of the fragments emptied into the duodenum have sub-mm sizes, while the largest fragment measures roughly 4.2 mm × 2.0 mm. These model predictions are consistent with earlier experimental studies suggesting that solid digesta must be reduced to dimensions below 3 mm to facilitate their passage through the pylorus to the duodenum.^6,8

6 Discussion

6.1 Gastric breakdown mechanisms

The results presented for the EB and EP solid beads collectively offer some key insights into the mechanisms responsible for the gastric breakdown of solid digesta. Essentially, four distinct classes of possible breakdown mechanisms can be identified from previous experimental and modelling literature:

(i) surface erosion due to chemical reactions catalysed by gastric enzymes and acids,^90,91

(ii) breakdown due to the pressure of fluid flow,^19,91

(iii) direct interaction between the deforming gastric wall and solid digesta,³⁷ and

(iv) collision with other solid food fragments and particles.^31,92

The first mechanism has not been included in this generation of the model. The remaining three mechanisms are all found to be highly active in the current study and contribute to the breakdown behaviour of food materials. Rapid ACWs in the terminal antral wall (TAC) and the resulting retropulsive flow (mechanisms iii and ii respectively) work in coordination to apply strong compressive stresses to the food materials leading to much stronger break down of the solids for a wide range of EB and EP mechanical properties. TAC induced mechanisms, such as radial loading (squeezing) from antral wall compression and axial loading from retropulsive flow, have been shown to play a more important role in the breakdown of solids, particularly in the initial stages of digestion.

The results from the current study also demonstrate that collisions between solid materials (mechanism iv) can have a substantial effect on solid breakdown but only when the solid materials are soft or weak (either having a low TS for EB materials or a low YS for the EP materials). This demonstrates that solid-to-solid collisions can contribute to solid breakdown primarily in the later stages of digestion when the solids have become sufficiently softened due to gastric reactions and diffusion of water inside the food structure (mechanism i). The general weakening of the structure of solid foods by mechanism i is generally synergistic with all three mechanical degradation mechanisms.

6.2 Implications of food rheology on nutrient availability and satiety

Despite the relatively short time frame used in the simulation, compared to total digestion times (minutes to hours), the current study reveals some key impacts of food rheology on nutrient absorption and satiety through the model predictions for total surface area of the beads and number of bead fragments. The total bead surface area at the end of the 60 s period appears to be substantially lower for the EP beads compared to the EB beads. The increase in surface area of the EP beads is mostly a consequence of their deformation into elongated and non-spherical shapes. In contrast, EB beads are fractured into numerous sub-mm sized fragments which contributes to a larger increase in total surface area, which is inversely proportional to their TS. The higher surface area of the EB beads can lead to stronger acidification and interaction with the digestive enzymes and subsequently faster and more complete nutrient release in the intestines.² Over a longer period, smaller EB bead fragments are likely to be emptied faster from the stomach through the gastric sieving phenomenon.^5–7 On the other hand, significantly deformed yet mostly non-fragmented EP beads would be retained in the stomach over longer timeframes for additional breakdown. This delayed gastric emptying of the EP beads can result in a higher level of satiety.¹⁶

6.3 SPH model limitations

While the new SPH model of gastric breakdown provides key insights into the deformation and fracturing dynamics of soft solid foods, several limitations need to be acknowledged. Firstly, the gastric walls are assumed to be rigid structures with high stiffness and strength since their deformations are kinematically prescribed. This assumption is only valid for scenarios involving soft solids with low stiffness and material strength. For solids with stiffness and strength equivalent or higher than the gastric wall, there will be substantial mechanical compliance and feedback on the gastric wall, which could affect the wall contraction dynamics including the TAC. The wall contractions can also be influenced by factors such as gastric content volume and viscosity, leading to more irregular motility patterns. Moreover, some fragmentation of the EB beads predicted by the SPH model could include a minor amount of numerical fracture that does not correlate with the actual damage accumulation in a real material. This can potentially influence the interpretation of fragmentation patterns. The smallest size of the fragments is also limited by the spatial resolution (initial particle separation) of the SPH model. Finally, the model employs a closed (no-flux) boundary condition at the duodenal boundary, resulting in negligible duodenal flow and relatively low pyloric discharge rates through the open sphincter. A high pyloric discharge rate can potentially alter the trajectory and deformation dynamics of the beads in the terminal antrum. These limitations of the SPH model provide context for the model predictions and highlights several areas for future work.

6.4 Future work

In the current study, 60 s of gastric mechanical breakdown was simulated with the goal of understanding the impact of various mechanical breakdown mechanisms on a model solid food. In future work, the model will be extended to run over longer durations, potentially of order of an hour which would allow inclusion and coupling of chemical and enzymic weakening (mechanism 1). Such extended timeframes will facilitate comparison and validation of the model predictions with results from other in vivo and in vitro^21,91 studies investigating the breakdown of similar spherical solid beads. This computational model with a longer timeframe could also serve as an in-silico tool to assist with designing foods and oral drug delivery systems with targeted gastric retention times.^93,94

The model also needs to be enhanced with extensions relating to:

(i) Tonic contraction of stomach: the volume of gastric digesta can influence the overall size and shape of the stomach due to changes in gastric tonic contraction, which can subsequently affect the flow of digesta in the antrum. Incorporating data from MRI studies⁹⁵ that quantify total gastric volume and gastric content volume at different times after a meal will be used to model the slow tonic contraction of gastric wall and its effect on dynamics of solid foods.

(ii) Mechano–chemical breakdown of food fragments in the stomach: in our ongoing work, the secretion of gastric acid from the stomach wall is being investigated, along with its mixing with rest of the digesta through advection and diffusion. The next crucial step in the development of this model is to include coupled representations of chemical and mechanical processes that lead to food breakdown, including water absorption, key chemical reactions involving gastric acid and enzymes, surface erosion, dissolution and the softening of solid foods.

(iii) Duodenal pressure: over longer durations, duodenal pressure can significantly impact the gastric emptying rate and flow within the stomach,⁹⁶ affecting the dynamics and mechanical breakdown of solid beads. The effect of food breakdown dependence on duodenal pressure and its variation needs to be investigated.

For sufficiently accurate modelling of mechanical breakdown, it will be important to also quantify the time and spatial variation in material properties of solid digesta during gastric digestion, corresponding to different gastric conditions, degrees of water absorption, and degradation of the food matrix structure. In future work, digesta samples will be characterised using data from laboratory tests such as uniaxial compression tests. Mathematical or empirical relationships can be developed to incorporate material property variations in the food materials as water is absorbed or chemical changes occur, which can be directly coupled with mechanical effects in the SPH framework.

These model additions will enable comprehensive simulations of gastric fluid dynamics and digestive life cycle of solid foods over a broad range of gastric conditions.

7 Conclusions

This study introduces a novel computational model of gastric breakdown of soft solid foods in a liquid medium driven by realistic three-dimensional motion of the gastric walls. The model incorporates peristaltic contraction waves, including Terminal Antral Contractions (TACs), which are rapid, high-amplitude peristaltic waves travelling in the antral wall near the pylorus. Smoothed Particle Hydrodynamics (SPH) is used to simulate the movement of the stomach walls, the resulting fluid dynamics of the liquid medium, and the motion, deformation and fracturing of solid beads. Two types of solid rheology are considered – elastic-plastic (EP) and elastic-brittle (EB), which represent the majority of solid food mechanical behaviours in the stomach.

The model shows how gastric wall contractions and fluid dynamics can determine the behaviour of solid foods. When the pyloric sphincter opens, flow is established through the sphincter towards the duodenum. Nearby solids are drawn towards the sphincter, with bodies that are larger than the sphincter opening able to obstruct the transpyloric flow. These cease moving as the contact force from the sphincter wall balances the fluid pressure gradient. When the pyloric sphincter closes, these jammed solids are pushed back towards the antrum. Rapid contraction of the terminal antrum during TAC generates a high-speed retrograde jet from the pylorus towards the antrum, reversing the fluid flow direction. Solids near the pylorus therefore follow loop-like trajectories driven by the alternating direction of flow to and from the pylorus. In contrast, solids located in the middle of the antrum have generally downward trajectories, causing them to settle near the bottom of the stomach.

The results demonstrate three distinct mechanisms causing mechanical breakdown of solid foods in the stomach –

(i) During the TAC, large solids present in the terminal antrum can be subjected to strong radial compression by the occluded wall. Simultaneously, the TAC induced retropulsive flow loads the solid in the direction opposite to that of the TAC waves. Due to this combined wall and fluid mechanical loading, a solid can be extruded into an elongated shape (EP) or fractured (EB), depending on its rheological behaviour. Decreasing the yield stress of EP materials increases the amount of plastic deformation, while decreasing the threshold strain (fracture strength) of EB materials leads to more frequent fracturing into a greater number of fragments.

(ii) Fluid motion from the peristaltic waves in the bottom of the stomach can lead to collisions between solids, which can be followed by compression of the solids by the moving antral wall leading to substantial change in their shape and some size reduction, depending on their rheology (EP or EB) and strength (yield stress or threshold strain).

(iii) During closure of the open sphincter, the sphincter wall can exert a deforming force on the contact area of the solids jamming the transpyloric fluid flow. Strong shear deformation can cause a thin surface layer of the solids to be pinched off, leading to the formation of very small fragments that are suspended in the adjacent fluid.

The model predictions highlight TACs as a major driver of solid breakdown, during which solid foods of different rheologies and a wide range of material strengths can be rapidly deformed and/or fragmented. Weakening of the food structure, due to chemical processes, is likely to strongly couple with this mechanical breakdown during subsequent TAC events. The study also shows that peristaltic wall movements and solid-to-solid collisions can contribute meaningfully to soft solid breakdown when the food structures are sufficiently weakened.

Author contributions

All authors contributed to the study conception and experimental design. S. G. performed the simulation experiments, analysed the data and had the primary responsibility of writing the manuscript. All authors contributed to the interpretation of data, critically reviewed and revised the manuscript, and approved the final manuscript.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The work was funded by the Active Integrated Matter (AIM) and Microbiomes for One Systems Health (MOSH) Future Science Platforms (FSPs) in CSIRO.

References

1. A. R. Mackie, B. H. Bajka, N. M. Rigby, P. J. Wilde, F. Alves-Pereira, E. F. Mosleth, A. Rieder, B. Kirkhus and L. J. Salt, Oatmeal particle size alters glycemic index but not as a function of gastric emptying rate, Am. J. Physiol.: Gastrointest. Liver Physiol., 2017, 313, G239–G246.
2. J. Nadia, J. Bronlund, R. P. Singh, H. Singh and G. M. Bornhorst, Structural breakdown of starch–based foods during gastric digestion and its link to glycemic response: In vivo and in vitro considerations, Compr. Rev. Food Sci. Food Saf., 2021, 20, 2660–2698.
3. F. Cisse, E. A. Pletsch, D. P. Erickson, M. Chegeni, A. M. Hayes and B. R. Hamaker, Preload of slowly digestible carbohydrate microspheres decreases gastric emptying rate of subsequent meal in humans, Nutr. Res., 2017, 45, 46–51.
4. P. Maljaars, H. Peters, D. Mela and A. Masclee, Ileal brake: a sensible food target for appetite control. A review, Physiol. Behav., 2008, 95, 271–281.
5. G. M. Bornhorst, K. Kostlan and R. P. Singh, Particle size distribution of brown and white rice during gastric digestion measured by image analysis, J. Food Sci., 2013, 78, E1383–E1391.
6. J. H. Meyer, J. Elashoff, V. Porter-Fink, J. Dressman and G. L. Amidon, Human postprandial gastric emptying of 1–3-millimeter spheres, Gastroenterology, 1988, 94, 1315–1325.
7. R. Vincent, A. Roberts, M. Frier, A. Perkins, I. MacDonald and R. Spiller, Effect of bran particle size on gastric emptying and small bowel transit in humans: a scintigraphic study., Gut, 1995, 37, 216–219.
8. B. P. Brown, K. Schulze-Delrieu, J. E. Schrier and M. M. Abu-Yousef, The configuration of the human gastroduodenal junction in the separate emptying of liquids and solids, Gastroenterology, 1993, 105, 433–440.
9. A. Santangelo, M. Peracchi, D. Conte, M. Fraquelli and M. Porrini, Physical state of meal affects gastric emptying, cholecystokinin release and satiety, Br. J. Nutr., 1998, 80, 521–527.
10. S. Gallier, S. M. Rutherfurd, P. J. Moughan and H. Singh, Effect of food matrix microstructure on stomach emptying rate and apparent ileal fatty acid digestibility of almond lipids, Food Funct., 2014, 5, 2410–2419.
11. H. Parker, C. L. Hoad, E. Tucker, C. Costigan, L. Marciani, P. Gowland and M. Fox, Gastric motor and sensory function in health assessed by magnetic resonance imaging: establishment of reference intervals for the Nottingham test meal in healthy subjects, Neurogastroenterol. Motil., 2018, 30, e13463.
12. M. Armand, P. Borel, C. Dubois, M. Senft, J. Peyrot, J. Salducci, H. Lafont and D. Lairon, Characterization of emulsions and lipolysis of dietary lipids in the human stomach, Am. J. Physiol.: Gastrointest. Liver Physiol., 1994, 266, G372–G381.
13. M. Armand, B. Pasquier, M. André, P. Borel, M. Senft, J. Peyrot, J. Salducci, H. Portugal, V. Jaussan and D. Lairon, Digestion and absorption of 2 fat emulsions with different droplet sizes in the human digestive tract, Am. J. Clin. Nutr., 1999, 70, 1096–1106.
14. H. Hornby, M. Collado-González, X. Zhang, N. Abrehart, M. Alshammari, S. Bakalis, A. Mackie and L. Marciani, Size and Number of Food Boluses in the Stomach after Eating Different Meals: Magnetic Resonance Imaging Insights in Healthy Humans, Nutrients, 2021, 13, 3626.
15. C. Spann, S. Al-Maliki, F. Boué, É. Lutton and F. Vidal, Interactive Visualisation of the Food Content of a Human Stomach in MRI, Comput. Graph. Vis. Comput. CGVC Eurographics Assoc., 2022, pp. 47–54.
16. L. Marciani, P. A. Gowland, A. Fillery-Travis, P. Manoj, J. Wright, A. Smith, P. Young, R. Moore and R. C. Spiller, Assessment of antral grinding of a model solid meal with echo-planar imaging, Am. J. Physiol.: Gastrointest. Liver Physiol., 2001, 280, G844–G849.
17. D. Dupont, M. Alric, S. Blanquet-Diot, G. Bornhorst, C. Cueva, A. Deglaire, S. Denis, M. Ferrua, R. Havenaar and J. Lelieveld, Can dynamic in vitro digestion systems mimic the physiological reality?, Crit. Rev. Food Sci. Nutr., 2019, 59, 1546–1562.
18. M. Minekus, The TNO gastro-intestinal model (TIM), Springer International Publishing, 2015.
19. F. Kong and R. P. Singh, A human gastric simulator (HGS) to study food digestion in human stomach, J. Food Sci., 2010, 75, E627–E635.
20. J. Nadia, H. Singh and G. M. Bornhorst, Evaluation of the performance of the human gastric simulator using durum wheat-based foods of contrasting food structure, Food Funct., 2023, 14, 6877–6895.
21. L. Chen, Y. Xu, T. Fan, Z. Liao, P. Wu, X. Wu and X. D. Chen, Gastric emptying and morphology of a ‘near real'in vitro human stomach model (RD-IV-HSM), J. Food Eng., 2016, 183, 1–8.
22. Y. Li, L. Fortner and F. Kong, Development of a Gastric Simulation Model (GSM) incorporating gastric geometry and peristalsis for food digestion study, Food Res. Int., 2019, 125, 108598.
23. J. Wang, P. Wu, M. Liu, Z. Liao, Y. Wang, Z. Dong and X. D. Chen, An advanced near real dynamic in vitro human stomach system to study gastric digestion and emptying of beef stew and cooked rice, Food Funct., 2019, 10, 2914–2925.
24. Z. Peng, P. Wu, J. Wang, D. Dupont, O. Menard, R. Jeantet and X. D. Chen, Achieving realistic gastric emptying curve in an advanced dynamic in vitro human digestion system: experiences with cheese—a difficult to empty material, Food Funct., 2021, 12, 3965–3977.
25. R. Berry, T. Miyagawa, N. Paskaranandavadivel, P. Du, T. R. Angeli, M. L. Trew, J. A. Windsor, Y. Imai, G. O'Grady and L. K. Cheng, Functional physiology of the human terminal antrum defined by high-resolution electrical mapping and computational modeling, Am. J. Physiol.: Gastrointest. Liver Physiol., 2016, 311, G895–G902.
26. S. Ishida, T. Miyagawa, G. O'Grady, L. K. Cheng and Y. Imai, Quantification of gastric emptying caused by impaired coordination of pyloric closure with antral contraction: a simulation study, J. R. Soc., Interface, 2019, 16, 20190266.
27. D. Ramkumar and K. Schulze, The pylorus, Neurogastroenterol. Motil., 2005, 17, 22–30.
28. A. Pal, K. Indireshkumar, W. Schwizer, B. Abrahamsson, M. Fried and J. G. Brasseur, Gastric flow and mixing studied using computer simulation, Proc. R. Soc. London, Ser. B, 2004, 271, 2587–2594.
29. K. Schulze, Imaging and modelling of digestion in the stomach and the duodenum, Neurogastroenterol. Motil., 2006, 18, 172–183.
30. V. Meves and J. Pohl, Upper Gastrointestinal Endoscopy: Examination Technique and Standard Findings, Video J. Encycl. GI Endosc., 2013, 1, 202–203.
31. C. Li, J. Xiao, X. D. Chen and Y. Jin, Mixing and emptying of gastric contents in human-stomach: A numerical study, J. Biomech., 2021, 118, 110293.
32. Z. Xue, M. J. Ferrua and P. Singh, Computational fluid dynamics modeling of granular flow in human stomach, Aliment. Hoy, 2012, 21, 3–14.
33. X. Liu, D. F. Fletcher and G. M. Bornhorst, A review of the use of numerical analysis in stomach modeling, J. Food Sci., 2024, 89, 3894–3916.
34. P. Trusov, N. V. Zaitseva and M. Kamaltdinov, A multiphase flow in the antroduodenal portion of the gastrointestinal tract: a mathematical model, Comput. Math. Methods Med., 2016, 2016, 1.
35. J. H. Seo and R. Mittal, Computational Modeling of Drug Dissolution in the Human Stomach, Front. Physiol., 2021, 2080.
36. J. H. Lee, S. Kuhar, J.-H. Seo, P. J. Pasricha and R. Mittal, Computational modeling of drug dissolution in the human stomach: Effects of posture and gastroparesis on drug bioavailability, Phys. Fluids, 2022, 34, 081904.
37. C. Skamniotis, C. H. Edwards, S. Bakalis, G. Frost and M. Charalambides, Eulerian-Lagrangian finite element modelling of food flow-fracture in the stomach to engineer digestion, Innovative Food Sci. Emerging Technol., 2020, 66, 102510.
38. S. M. Harrison, P. W. Cleary and M. D. Sinnott, Investigating mixing and emptying for aqueous liquid content from the stomach using a coupled biomechanical-SPH model, Food Funct., 2018, 9, 3202–3219.
39. X. Liu, S. M. Harrison, D. F. Fletcher and P. W. Cleary, Numerical simulation of buoyancy–driven flow in a human stomach geometry: Comparison of SPH and FVM models, Appl. Math. Models, 2023, 124, 367–392.
40. X. Liu, S. M. Harrison, P. W. Cleary and D. F. Fletcher, Evaluation of SPH and FVM models of kinematically prescribed peristalsis-like flow in a tube, Fluids, 2022, 8, 6.
41. P. W. Cleary, S. M. Harrison, M. D. Sinnott, G. G. Pereira, M. Prakash, R. C. Cohen, M. Rudman and N. Stokes, Application of SPH to single and multiphase geophysical, biophysical and industrial fluid flows, Int. J. Comput. Fluid Dyn., 2021, 35, 22–78.
42. P. W. Cleary, M. Prakash, J. Ha, N. Stokes and C. Scott, Smooth particle hydrodynamics: status and future potential, Prog. Comput. Fluid Dyn. Int. J., 2007, 7, 70–90.
43. J. Monaghan, 76399: Smoothed particle hydrodynamics, Rep. Prog. Phys., 2005, 68, 1703–1759.
44. P. W. Cleary and R. Das, The potential for SPH modelling of solid deformation and fracture, in Proceedings of the IUTAM Symposium, Springer, 2008, pp. 287–296.
45. S. M. Harrison and P. W. Cleary, Towards modelling of fluid flow and food breakage by the teeth in the oral cavity using smoothed particle hydrodynamics (SPH), Eur. Food Res. Technol., 2014, 238, 185–215.
46. M. Sinnott, S. Harrison and P. Cleary, A particle-based modelling approach to food processing operations, Food Bioprod. Process., 2021, 127, 14–57.
47. P. Cleary, M. Sinnott, B. Hari, S. Bakalis and S. Harrison, in Modeling Food Processing Operations, Elsevier, 2015, pp. 255–305.
48. P. W. Cleary, S. M. Harrison and M. D. Sinnott, in Digital Human Modeling and Medicine, Elsevier, 2023, pp. 361–424.
49. J. J. Monaghan, Smoothed particle hydrodynamics, Rep. Prog. Phys., 2005, 68, 1703.
50. J. J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys., 1994, 110, 399–406.
51. G. K. Batchelor, An introduction to fluid dynamics, Cambridge university press, 2000.
52. J. P. Gray, J. J. Monaghan and R. Swift, SPH elastic dynamics, Comput. Methods Appl. Mech. Eng., 2001, 190, 6641–6662.
53. J. J. Monaghan, SPH without a tensile instability, J. Comput. Phys., 2000, 159, 290–311.
54. R. Das and P. W. Cleary, Evaluation of accuracy and stability of the classical SPH method under uniaxial compression, J. Sci. Comput., 2015, 64, 858–897.
55. M. L. Wilkins, Calculation of elastic-plastic flow, Lawrence Radiation Lab., California. Univ., Livermore, 1963.
56. D. E. Grady and M. E. Kipp, Continuum modelling of explosive fracture in oil shale, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1980, 17, 147–157.
57. W. Benz and E. Asphaug, Impact simulations with fracture. I. Method and tests, Icarus, 1994, 107, 98–116.
58. H. Melosh, E. Ryan and E. Asphaug, Dynamic fragmentation in impacts: Hydrocode simulation of laboratory impacts, J. Geophys. Res.: Planets, 1992, 97, 14735–14759.
59. C. Thornton, S. J. Cummins and P. W. Cleary, An investigation of the comparative behaviour of alternative contact force models during inelastic collisions, Powder Technol., 2013, 233, 30–46.
60. S. J. Cummins and P. W. Cleary, Using distributed contacts in DEM, Appl. Math. Models, 2011, 35, 1904–1914.
61. D. R. Vyas, S. J. Cummins, M. Rudman, P. W. Cleary, G. W. Delaney and D. V. Khakhar, Collisional SPH: A method to model frictional collisions with SPH, Appl. Math. Models, 2021, 94, 13–35.
62. R. Treier, A. Steingoetter, D. Weishaupt, O. Goetze, P. Boesiger, M. Fried and W. Schwizer, Gastric motor function and emptying in the right decubitus and seated body position as assessed by magnetic resonance imaging, J. Magn. Reson. Imaging, 2006, 23, 331–338.
63. Y.-K. Lee, H.-S. Kim and J.-Y. Park, The case study of masticatory force with food from full skull and partial model, Int. J. Precis. Eng. Manuf., 2017, 18, 1455–1462.
64. J. W. Kim, S. S. Shin, S. H. Heo, H. S. Lim, N. Y. Lim, Y. K. Park, Y. Y. Jeong and H. K. Kang, The role of three-dimensional multidetector CT gastrography in the preoperative imaging of stomach cancer: emphasis on detection and localization of the tumor, Korean J. Radiol., 2015, 16, 80–89.
65. R. K. Goyal, Y. Guo and H. Mashimo, Advances in the physiology of gastric emptying, Neurogastroenterol. Motil., 2019, 31, e13546.
66. F. Kong and R. P. Singh, Disintegration of solid foods in human stomach, J. Food Sci., 2008, 73, R67–R80.
67. M.-S. Hur, S. Lee, T. M. Kang and C.-S. Oh, The three muscle layers in the pyloric sphincter and their possible function during antropyloroduodenal motility, Sci. Rep., 2021, 11, 1–13.
68. K. Indireshkumar, J. G. Brasseur, H. Faas, G. S. Hebbard, P. Kunz, J. Dent, C. Feinle, M. Li, P. Boesiger and M. Fried, Relative contributions of “pressure pump” and “peristaltic pump” to gastric emptying, Am. J. Physiol.: Gastrointest. Liver Physiol., 2000, 278, G604–G616.
69. J. Munk, R. Gannaway, M. Hoare and A. Johnson, Direct measurement of pyloric diameter and tone in man and their response to cholecystokinin, in Gastrointestinal motility in health and disease, Springer, 1978, pp. 349–359.
70. N. Jagtap, R. Kalapala and D. N. Reddy, Assessment of pyloric sphincter physiology using functional luminal imaging probe in healthy volunteers, J. Neurogastroenterol. Motil., 2020, 26, 391.
71. A. Keet Jr, The prepyloric contractions in the normal stomach, Acta Radiol., 1957, 48, 413–424.
72. M. Czerner, L. S. Fellay, M. P. Suárez, P. M. Frontini and L. A. Fasce, Determination of elastic modulus of gelatin gels by indentation experiments, Procedia Mater. Sci., 2015, 8, 287–296.
73. Q. Chang, M. A. Darabi, Y. Liu, Y. He, W. Zhong, K. Mequanin, B. Li, F. Lu and M. M. Xing, Hydrogels from natural egg white with extraordinary stretchability, direct-writing 3D printability and self-healing for fabrication of electronic sensors and actuators, J. Mater. Chem. A, 2019, 7, 24626–24640.
74. A. Sun and S. Gunasekaran, Measuring rheological characteristics and spreadability of soft foods using a modified squeeze–flow apparatus, J. Texture Stud., 2009, 40, 275–287.
75. A. Sun and S. Gunasekaran, Yield stress in foods: measurements and applications, Int. J. Food Prop., 2009, 12, 70–101.
76. A. S. Giacomozzi, J. Benedito, A. Quiles, J. V. García-Pérez and M. E. Dalmau, Ultrasonic monitoring of softening in solid foods during in vitro gastric digestion, J. Food Eng., 2024, 374, 112033.
77. K. C. Drechsler and G. M. Bornhorst, Modeling the softening of carbohydrate-based foods during simulated gastric digestion, J. Food Eng., 2018, 222, 38–48.
78. J. Hempe, R. Lauxen and J. Savage, Rapid determination of egg weight and specific gravity using a computerized data collection system, Poult. Sci., 1988, 67, 902–907.
79. E. A. H. Akpa, H. Suwa, Y. Arakawa and K. Yasumoto, Smartphone-based food weight and calorie estimation method for effective food journaling, SICE J. Control Meas. Syst. Integr., 2017, 10, 360–369.
80. J. Meyer, Y. Gu, J. Elashoff, T. Reedy, J. Dressman and G. Amidon, Effects of viscosity and fluid outflow on postcibal gastric emptying of solids, Am. J. Physiol.: Gastrointest. Liver Physiol., 1986, 250, G161–G164.
81. Y. Jin, P. J. Wilde, Y. Hou, Y. Wang, J. Han and W. Liu, An evolving view on food viscosity regulating gastric emptying, Crit. Rev. Food Sci. Nutr., 2023, 63, 5783–5799.
82. M. A. van Nieuwenhoven, E. M. Kovacs, R.-J. M. Brummer, M. S. Westerterp-Plantenga and F. Brouns, The effect of different dosages of guar gum on gastric emptying and small intestinal transit of a consumed semisolid meal, J. Am. Coll. Nutr., 2001, 20, 87–91.
83. L. Marciani, P. A. Gowland, R. C. Spiller, P. Manoj, R. J. Moore, P. Young and A. J. Fillery-Travis, Effect of meal viscosity and nutrients on satiety, intragastric dilution, and emptying assessed by MRI, Am. J. Physiol.: Gastrointest. Liver Physiol., 2001, 280, G1227–G1233.
84. R. Gopirajah, K. P. Raichurkar, R. Wadhwa and C. Anandharamakrishnan, The glycemic response to fibre rich foods and their relationship with gastric emptying and motor functions: An MRI study, Food Funct., 2016, 7, 3964–3972.
85. J. Ha and P. W. Cleary, Comparison of SPH simulations of high pressure die casting with the experiments and VOF simulations of Schmid and Klein, Int. J. Cast Met. Res., 2000, 12, 409–418.
86. M. Robinson, P. Cleary and J. Monaghan, Analysis of mixing in a twin cam mixer using Smoothed Particle Hydrodynamics, AIChE J., 2008, 54, 1987–1998.
87. M. Prakash, P. W. Cleary, J. Ha, M. N. Noui-Mehidi, H. Blackburn and G. Brooks, Simulation of suspension of solids in a liquid in a mixing tank using SPH and comparison with physical modelling experiments, Prog. Comput. Fluid Dyn. Int. J., 2007, 7, 91–100.
88. R. Das and P. Cleary, Simulation of uniaxial compression tests and stress wave propagation using smoothed particle hydrodynamics, Int. J. Numer. Methods Eng. Submitted Publ., 2013, 238, 185–215.
89. R. Das and P. Cleary, A mesh-free approach for fracture modelling of gravity dams under earthquake, Int. J. Fract., 2013, 179, 9–33.
90. J. Nadia, A. G. Olenskyj, P. Subramanian, S. Hodgkinson, N. Stroebinger, T. G. Estevez, R. P. Singh, H. Singh and G. M. Bornhorst, Influence of food macrostructure on the kinetics of acidification in the pig stomach after the consumption of rice-and wheat-based foods: Implications for starch hydrolysis and starch emptying rate, Food Chem., 2022, 394, 133410.
91. M. Vardakou, A. Mercuri, S. A. Barker, D. Q. Craig, R. M. Faulks and M. S. Wickham, Achieving antral grinding forces in biorelevant in vitro models: comparing the USP dissolution apparatus II and the dynamic gastric model with human in vivo data, AAPS PharmSciTech, 2011, 12, 620–626.
92. C. Swackhamer, Breakdown of Solid Foods During in Vitro Gastric Digestion, UC Davis, 2022.
93. I. Locatelli, N. Nagelj Kovačič, A. Mrhar and M. Bogataj, Gastric emptying of non-disintegrating solid drug delivery systems in fasted state: relevance to drug dissolution, Expert Opin. Drug Delivery, 2010, 7, 967–976.
94. J. Tripathi, P. Thapa, R. Maharjan and S. H. Jeong, Current state and future perspectives on gastroretentive drug delivery systems, Pharmaceutics, 2019, 11, 193.
95. S. Banerjee, A. Pal and M. Fox, Volume and position change of the stomach during gastric accommodation and emptying: A detailed three–dimensional morphological analysis based on MRI, Neurogastroenterol. Motil., 2020, 32, e13865.
96. S. Ghosh, P. W. Cleary and S. M. Harrison, Calibration of a physics-based human gastric emptying model, in Proceedings of the International Congress on Modelling and Simulation, MODSIM, MSSANZ, Darwin, 2023, vol. 25, pp. 264–270.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4fo04217d

---