1 Introduction

The red blood cell (RBC, also referred to as erythrocyte) is the most common type of blood cells, occupying about 45% of the total blood volume for a man and approximately 40% for a woman. The RBCs in the healthy state resemble biconcave-shaped disks with a diameter within [6μm, 8μm] and a thickness of about 2μm^[1]. The interactions among the RBCs in plasma is a quite fundamental biomechanical topic because it plays an important role in the rheological behavior of blood flows^[2-5]. More specifically, plasma macromolecules (mainly fibrinogen) cause the phenomenon of RBC aggregation, which is mainly responsible for the non-Newtonian nature of the fluid^{[1, 3-4, 6]}. In order to investigate such mechanisms, great efforts have been put into the simulations for RBC-laden flows in recent years^[7-10], in which the RBCs are modeled as deformable solids^[11-12] or even rigid bodies^[13-16]. Such simulations are interesting and meaningful because they provide ways to deepen our insight into the physics lying behind blood flow properties from a microscale perspective^[17]. Haematocrit (Ht or HCT) is widely verified to be a crucial factor dominating blood rheodynamics^[18-19]. An important question regarding simulated RBC flows is how many RBCs shall be put into the simulation to generate blood with certain haematocrit, say for example 40% for women? Current simulations usually calculate the value of haematocrit R_Ht by dividing the RBC volume with the total volume of the mixture, i.e.,

(1)

where N is the number of the RBCs, V_RBC is the volume of a single RBC, and V₀ is the volume of the blood, i.e.,

(1)

However, in practical clinics, the packed cell volume (PCV) is usually determined as follows: firstly, the heparinized blood is centrifuged in a capillary tube (also known as a microhematocrit tube) at 10 000 r/min for five minutes. This separates the blood into layers, and the volume of the packed RBCs divided by the total volume of the blood sample gives the PCV^[1]. Therefore, the measured volume of the RBCs is not simply the sum of the geometric volume of each RBC, but with the voids (filled with plasma) among the included RBCs. Namely,

(2)

where V_void is the volume of the voids per cell among the packed RBCs. This raises another unanswered question: what is the void ratio of the packed RBCs?

The void ratio is an important geometrical and physical property for granular solids, which is also the key to calculate the packing density or packing constant. The packing constant of a geometric body is the largest average density achieved by the packing arrangements of the congruent copies of the body. For instance, if the packing constant of the uniform size spheres is 0.740 480 5, it is 0.906 9 for the bi-infinite cylinders^[20]. For most bodies, the value of the packing constant is unknown^[20]. The density of the randomly packed solids is also important for both scientific problems and many practical applications. The numerical simulations based on the discrete element method are effective tools in solving this problem, whose applications range from regular-shaped disks^[21] and spheres^[22] to polyhedral^[23] or even elastic wires^[24]. Such simulations or experiments are valuable in helping us to understand the granular or powder systems from a statistical point of view^[25-29]. Packed human RBCs, although physically deformable, fall into such systems, at least when their deformability simplicity is ignored^[13]. However, up to now, even for the simplest case (uniform size without friction, adhesion, or deformation), an accurate description of the void ratio or packed density for randomly packed human RBCs remains unavailable.

In this paper, a series of numerical simulations utilizing a finite-discrete element method (FDEM) are carried out, and a relatively accurate description on the packed volume and the void ratio of human RBCs is achieved. Of course, in practical clinics, the packing behavior of RBCs is much more complex, considering the effects from the aggregation, the gradation, the deformation, the storage methods, etc. The objective of this paper is to study the packing behavior of human RBCs under the most idealized situation from a mathematical perspective.

2 Methodology 2.1 Modeling RBC motion and interactions

To simulate the movement and collision of RBCs in the packing process, we use the in-house C code developed by Munjiza^[30] and Munjiza et al.^[31]. The method is based on the idea of the discrete element method with an implementation through the finite element framework, i.e., the combined FDEM. The code comprises a set of C libraries incorporating the latest breakthroughs in discontinua simulations, and is capable of modeling the movement and collisions of millions of particles of different shapes and sizes. During the simulations, the RBCs are treated as solids meshed with tetrahedron elements (see Fig. 1).

In the FDEM framework, the deformation and movement of the solids can be solved by a simple implementation of the second Newton's law for a node as follows:

(3)

where F is the total force acting on each node, M is the mass represented by the node, a is the acceleration to be calculated, and x is the node coordinates which will be updated at each time step.

The total force F in the present simulation is composed of the internal force F_I and the collision force F_C, i.e.,

The internal force F₁ is calculated based on the linear elastic principle according to the deformation of each element. Here, a very high shear modulus (10⁴ times of the normal value G₀=2.57× 10⁵ Pa for healthy RBCs^{[2, 14]}) is used to agree with the previously mentioned solid body assumption for RBCs. It corresponds to the RBC under low shear when its deformation is negligible^[2]. The collision force F_C is calculated based on a penalty method^[30]. The cell-cell interactions are modeled as a repulsive force to push them away when two cells have collisions, i.e., the collision force. The collision stiffness is set to be 1/10 of the shear modulus according to the recommendation of the FDEM method^[30]. The cell-wall interactions are modeled in a similar way except that the wall is non-deformable. The gravitational force is applied on the RBCs by adding a gravitational acceleration of 9.8 m·s^-2 directly to each element node. For the sake of conciseness, the details of the methodology, except the collision strategy, are not presented here, and readers may refer to previous work^[30-31] for further information. The code has a wide application, and is well validated through the transport simulations of both sediment particles with a very high modulus^[32-34] and highly deformable RBCs^[2], whose simulated deformation under the stretching forces agrees well with the experimental data in Ref. [35].

2.2 Modeling the RBC shape

The static shape of a healthy human RBC can usually be described with a biconcave discoid shape. The x, y coordinates of the cross-sectional profile of a RBC with a radius of R=1 can usually be described as follows^[36]:

(4)

where C₀, C₁, and C₂ are coefficients. Such a biconcave discoid shape has been proven by both interference microscopy observation and theoretical analysis based on the bending elasticity of the surrounding membrane. For normal cells, C₀ =0.207, C₁ =2.002, and C₂ =-1.122 are widely adopted in numerical simulations^{[12, 37-40]}.

3 Results 3.1 Geometry of typical RBCs described with a biconcave function

The biconcave discoid shape under the static state has been proven by both interference microscopy observation and theoretical analysis. For instance, by assuming that the shape of a RBC is controlled by the bending elasticity of the surrounding membrane, Deuling and Helfrich^[41] fitted theoretical shapes to the contours in Ref. [36], which were determined by interference microscopy, and obtained very good agreement for disc shapes. The biconcave discoid function presented by Evans and Fung^[36] can describe a wide range of RBC shapes considering the coefficients C₀, C₁, and C₂. Although for most cases, the coefficients are close to certain values for normal human RBCs, abnormal RBC shapes, with quite different coefficients in the biconcave discoid function, are also widely reported^[42-44]. The shape of a RBC contributes greatly to the clinical diagnosis with respect to the blood diseases^[45-46].

For the convenience of applications in RBC modeling, the effects of the coefficients of the biconcave discoid function on the RBC shapes are presented (see Fig. 2). Considering that the coordinate y in the thickness direction of a RBC profile has a relatively simple relationship with C₀, i.e., y is linearly proportional to C₀ for a given x, only the coefficients for higher order terms, i.e., C₁ and C₂, are investigated. Figure 2 shows that only the combinations of moderate values of C₁ and C₂ may represent relatively normal RBC shapes, the most common shapes under the healthy state^[45]. Moreover, small values of C₁ and C₂ represent relatively flat RBC shapes, while large values represent relatively thick RBC shapes. Here, the following coefficients C₀ =0.207, C₁ =2.002, and C₂ =-1.122 are adopted for a RBC shape of a standard cell, which are widely adopted in RBC simulations^[14]. The RBC with a deviation less than 20% from the standard cell is defined as a normal cell. The RBCs with 20% larger thickness are thick cells, and the RBCs with 20% smaller thickness are thin cells.

Two important geometric parameters are considered to describe the RBC shapes, i.e., the volume and the surface area. For the RBCs with given C₀, C₁, and C₂, the volume V can be easily obtained by the following integration:

(5)

By substituting x with the triangular functions and using a series of mathematical reduction, we obtain

(6)

The volume of the RBC is finally obtained as follows:

(7)

For the surface area of the RBCs, the following integration can be used:

(8)

Since Eq. (8) is too complex to be integrated analytically, the numerical integration using Simpson's method is used.

The effects of the coefficients C₁ and C₂ on the RBC volumes and surface areas are shown in Fig. 3. It can be seen that the RBC volume is linearly proportional to both C₁ and C₂, which agrees with the analytical solution (see Eq. (7)). For the RBC surface area, the relationship is basically linear except the regime with relatively small C₁ and C₂.

To investigate the effects of the RBC shape on the packing behavior, three typical RBC shapes are selected. The geometrical descriptions are shown in Table 1, and the cross-sectional profiles are shown in Fig. 4. It can be seen from Table 2 that the surface area-to-volume ratio of any of the three RBC shapes is much higher than that of spheres with equivalent radii, which is 3R^-1. Thick cells tend to exhibit a higher surface area-to-volume ratio than flat ones. A high surface area-to-volume ratio is usually considered to be favorable for RBCs in the exchange of oxygen and carbon dioxide^[1]. Take a RBC with a diameter of 7μm for example, the surface area-to-volume ratio for a representative normal cell is 1.585, which is close to the well measured values in experiments^[47-48].

3.2 Effects of the mesh resolution on the accuracy of the RBC shape

In most numerical simulations of RBC suspension flows, individual RBCs are usually treated as the deformable membranes of the hyper-elastic material^{[38-39, 49]}. Such membranes are usually meshed with triangular elements. Sometimes, tetrahedral elements can also be adopted when RBCs are modeled as deformable solids or shells^[2]. Considering that the numbers of RBCs in such simulations are usually exceptionally large, it is important to use elements as less as possible in meshing each RBC. However, what is the lowest limit for the RBC to yield an acceptable shape resolution from a geometric view? To answer this question, a series of tests have been carried out for the RBCs meshed with various resolutions. Typical flat, thick, and normal RBC shapes are considered with the surface element number ranging from 100 to 7 000 (see Fig. 5). The meshes are generated using the GAMBIT software based on a criterion that all nodes are attached to the surface generated by using a biconcave discoid function. It can be seen that the meshed cell volume increases with the number of the surface elements. For RBCs with surface elements less than 200, the meshed cell volume is nearly 10% less than the actual one. To achieve a geometric description for the RBC with a volume deviation less than 2.5%, the surface should be meshed with at least 700 surface triangular elements. The deviation can be further reduced to 0.5% when the surface elements exceed 3 000.

The mesh size convergence tests also show that, in order to make the RBC volume deviation lower than 2.5%, the element size shall be smaller than 0.08D (D is the cell diameter) with the surface meshed with over 700 triangular elements or volume meshed with over 3 000 tetrahedrons.

Based on the above tests on the mesh resolutions of RBCs, a uniform size 0.08D is adopted in all of the following numerical simulations in the RBC packing process. The above tests provide meaningful references on the spatial discretization for other numerical simulations with the involved RBC shapes. Recent progress also shows that the mesh size may have even more significant effects when the bending forces are focused^[50].

3.3 Packing process of the RBCs and measurement of the void ratio

In order to obtain the packed volume and the void raito of the RBCs described with the biconcave discoid function, DEM simulations are carried out for free falling RBCs in square-sectioned tubes (see Fig. 6). The RBCs are initially placed at random positions within the tubes. During the simulations, the RBCs fall freely under the gravitational force, and collide with each other and surrounding the walls. Due to the finite element approximation of the solids in the FDEM method, the RBCs are slightly compressible, which leads to oscillations of the packing height after the initial settlement. Such oscillations die out quickly under the action of the numerically introduced viscous damping, and the final packing height h can be obtained. Assume that n RBCs with a volume V_RBC are involved in the packing process and the edge size of the tube cross section is a. Then, the void ratio e of packed RBCs is

(9)

Since a very high shear modulus is adopted in the simulation for RBCs, the effect of the RBC deformation on the final packed RBC height is neglected here.

3.4 Effects of the RBC number and arrangement on the accuracy of the packed volume

In practice, the packed volume of solids is unavoidably affected by the shape of the container due to the near-wall effect^{[21, 28]}. Camenen et al. (2012)^[23] studied numerically the effects of confinement on the solid fraction and the structure of three-dimensional random close-packed granular materials subject to gravity, and showed that the solid fraction decreased linearly with the increase in the confinement, no matter the grain shape. Moreover, the number of the solids might also affect the packed volume since the uncertainty of a smaller number of solids might be higher due to the randomness. The primary objective of this paper is to work out the packed volume of RBCs as a geometrical property of the RBCs themselves. Theoretically, the simulation shall be carried out in an infinitely large domain with infinite RBC cells. In practice, we carry out a series of tests with gradually increasing container sizes and RBC numbers to check the convergence of the final results. The number of the RBCs is from 32 (2⁵) to 4 096 (2¹²) (see Fig. 7).

Each RBC packing simulation is repeated 4 times with different randomized initial setups, based on which a mean packed volume is calculated. Standard deviations are used to estimate the uncertainty of the results (see Fig. 8). It can be seen that both the packed volume and the computed volume deviation decrease with the increase in the RBC number. The estimated packed volume per RBC is 0.275D³ with a deviation around ± 3.3% when 4 096 RBCs are involved. For 32 RBCs, the computed packed volume is 23.6% larger, and the uncertainty is around 12 times higher. To obtain an estimation for the packed volume with the uncertainty lower than ± 5%, at least 1 024 RBCs shall be involved. All of the simulations exhibit a majority of near-edge RBCs parallel to the container wall. This is due to the fact that the RBCs parallel to the wall are more stable compared with those perpendicular to the wall. However, when the RBC amount is great enough, the wall effect is ignorable.

To verify that the packed volume of the RBCs based on the simulations is free from both the initial conditions, e.g., RBC arrangement, and the boundary conditions, e.g., container walls and container size, additional packing simulations for initially well-mixed RBCs are carried out with the periodic boundary conditions in the planar directions (see Fig. 9). According to the statistical results, the packed volume of a RBC is 0.276, which is very close the above packing results using 1 024 RBCs and falls well within the confidence interval shown in Fig. 8.

3.5 Effects of RBC shapes on the packed volumes

The above DEM simulations are carried out for the RBCs with three different shapes, i.e., thick cells, normal cells, and flat cells (see Fig. 4). 1 024 RBCs are involved in each simulation. The final state of the packed cells is shown in Fig. 10.

Basic on the simulation results, the packing parameters of typical human RBCs are obtained (see Table 2). The packed volumes of each RBC for flat, normal, and thick cells are 0.255, 0.275, and 0.324 (measured in D³), respectively, and the void ratios are 22.95%, 28.45%, and 23.20%, respectively. Besides, by packing the RBCs with side faces well aligned, the RBCs can be packed to the closest state, and the packing constants of the RBCs can be obtained, which are 0.798, 0.755, and 0.792 for the flat, normal, and thick cells, respectively.

3.6 Effects of the RBC deformability on the packed volumes

Human RBCs are highly deformable. Although RBCs are usually deemed as incompressible and their volume keeps constant during the deformations, the shear modulus may have effects on the void ratio of packed RBCs. In order to investigate the effects of RBC deformation on the packing void ratio, a series of DEM simulations are carried out for the RBCs with different shear moduli from G₀× 10^-4 to G₀× 10⁴, where G₀ is a reference shear modulus with a value of 2.57× 10⁵ Pa for healthy RBCs. A gravitational acceleration g=9.8 m· s^-2 is adopted. To save the computing time, only 128 RBCs are used in each simulation. The void ratio e is normalized with e₀, which is the void ratio for the solid RBCs with a shear modulus as high as G₀×10⁴.

Generally speaking, the void ratio increases with the increase in the shear modulus (see Fig. 11). The RBCs with a shear modulus 10⁴ times smaller than the normal value yield a 35.03% reduction in the packed void ratio. The RBCs near the bottom wall show a significant deformation. The biconcave shape of the RBCs is unrecognizable, and the cells closely contact each other. The RBCs with a shear modulus 10³ times smaller than the normal value yield a 11.05% reduction. The packing void ratios of RBCs with a normal and larger shear modulus show little differences. This indicates that the deformation of "real" human RBCs under normal gravitational force is negligible, while justifies the solid body assumption for the RBCs in the previous simulations, which utilizes a shear modulus of G₀×10⁴.

4 Conclusions and discussion

Numerical simulations using the combined FDEM are carried out to study the packing behavior of human RBCs. The RBCs are modeled as significantly idealized no-friction solid bodies. Based on the simulation results, we obtain the following conclusions:

(ⅰ) The packed volumes of each RBC for flat, normal, and thick cells, defined with a biconcave discoid function with different coefficients, are 0.255, 0.275, and 0.324 (measured in D³, where D is the cell diameter), respectively, and the void ratios are 22.95%, 28.45%, and 23.20%, respectively.

(ⅱ) To make the RBC volume deviation less than 2.5%, the element size shall be smaller than 0.08D with the surface meshed with over 700 triangular elements or volume meshed with over 3 000 tetrahedrons.

(ⅲ) To obtain an estimation for the packed volume with uncertainty lower than 5%, at least 1 024 RBCs shall be involved.

(ⅳ) The analytical solution for the RBC volume is

The numerical solutions for the RBC surface areas are presented by the integration of the RBC profiles. The above information adds to our new knowledge on the geometric feature of human RBCs described by the biconcave discoid function, and can be directly used to work out the number of RBCs required to generate the blood flow with certain haematocrits in a more accurate way.

Note that the analysis on the RBC packing in this paper is limited in quite idealized conditions. In practical clinics, RBC packing is much more complex. The packing behavior can be significantly affected by the RBC deformation, aggregation, red cell distribution width, etc. The discussion in this paper focuses on the RBC packing from a purely geometrical perspective. The findings here are particularly informative for the community of RBC simulations using numerical methods.