1. Introduction
The assessment of cell deformability as a mechanical biomarker has significant value in the study of various cellular processes (Starodubtseva Reference Starodubtseva2011; Khani et al. Reference Khani, Tafazzoli-Shadpour, Rostami, Peirovi and Janmaleki2014) and disease diagnosis (Costa Reference Costa2003; Suresh et al. Reference Suresh, Spatz, Mills, Micoulet, Dao, Lim, Beil and Seufferlein2005). Human red blood cells (RBCs) are responsible for transporting oxygen and need to pass repeatedly through capillaries, during which they undergo a high degree of deformation. However, a number of pathophysiological states can significantly reduce the deformability of RBCs, including conditions such as sickle cell anaemia (Dao, Lim & Suresh Reference Dao, Lim and Suresh2003) and malaria infection (Dondorp et al. Reference Dondorp, Angus, Chotivanich, Silamut, Ruangveerayuth, Hardeman, Kager, Vreeken and White1999), which in turn have a detrimental effect on blood circulation. In addition, ATP (adenosine triphosphate) is an important signalling molecule in the blood microcirculation and is released when RBCs are subjected to mechanical stress and cause sufficient deformation (Zhang et al. Reference Zhang, Shen, Hogan, Barakat and Misbah2018). Therefore, the study of erythrocyte deformation has potential clinical significance in the diagnosis of disease.
Over the years, the development of microelectromechanical systems has given rise to various methods for quantifying cell deformability. These methods can be divided into two categories, namely in situ and flow-based approaches (Xie et al. Reference Xie, Nama, Li, Mao, Huang, Zhao, Costanzo and Huang2016). In situ techniques, such as atomic force microscopy (Haase & Pelling Reference Haase and Pelling2015), micropipette aspiration (Pierrat, Brochard-Wyart & Nassoy Reference Pierrat, Brochard-Wyart and Nassoy2004), laser cavitation (Quinto-Su et al. Reference Quinto-Su, Kuss, Preiser and Ohl2011; Tandiono et al. Reference Tandiono, Klaseboer, Ohl, Ow, Choo, Li and Ohl2013) and optical tweezers (Dao et al. Reference Dao, Lim and Suresh2003) are usually tested at low speeds. In contrast, flow-based methods (Bow et al. Reference Bow, Pivkin, Diez-Silva, Goldfless, Dao, Niles, Suresh and Han2011; Gossett et al. Reference Gossett, Tse, Lee, Ying, Lindgren, Yang, Rao, Clark and Di Carlo2012; Zheng et al. Reference Zheng, Nguyen, Wang and Sun2013) usually only capture the deformation state of cells as they pass through a detection point. The combination of ultrasound and microfluidics has emerged as an attractive solution that combines high throughput, increased data recording, non-contact characteristics and various other advantages (Aubert et al. Reference Aubert, Wunenburger, Valier-Brasier, Rabaud, Kleman and Poulain2016; Xie, Bachman & Huang Reference Xie, Bachman and Huang2019).
The propagation of ultrasound waves through a medium not only triggers the time-harmonic oscillation, but also induces a slow time-averaged response due to acoustic nonlinear effects, including acoustic radiation force and acoustic streaming. Once suspended cells are introduced into the standing waves, they are subjected to the acoustic radiation force directed towards the acoustic pressure node. Cells trapped at the acoustic pressure node are deformed due to unevenly distributed time-averaged stress on the membrane surface. Based on this, ultrasonic standing waves have been successfully applied in biological systems for trapping (Jakobsson, Antfolk & Laurell Reference Jakobsson, Antfolk and Laurell2014; Rocha et al. Reference Rocha, Silva, Sales, D’amato, Leite and Silva2023; Li et al. Reference Li, Nunn, Brumley, Sader and Collis2024; Sujith, Malik & Sen Reference Sujith, Malik and Sen2024) and deforming (Mishra, Hill & Glynne-Jones Reference Mishra, Hill and Glynne-Jones2014; Wijaya et al. Reference Wijaya, Mohapatra, Sepehrirahnama and Lim2016; Silva et al. Reference Silva, Tian, Franklin, Wang, Han, Mann and Drinkwater2019) biological microparticles. The acoustic deformation of cells is quantified by deformation index and relaxation time, which reveals the viscous and elastic properties of cells (Link & Franke Reference Link and Franke2020). While microfluidic parameters induced by cells, such as flow patterns and velocity magnitudes, can serve as label-free proxies for measuring the mechanical properties of cells (Salari et al. Reference Salari, Appak-Baskoy, Ezzo, Hinz, Kolios and Tsai2020).
To explain the experimental phenomenon of acoustic-induced deformation, several theoretical and numerical studies have made progress in recent years. Mishra et al. (Reference Mishra, Hill and Glynne-Jones2014) established a finite element model for the deformation of swollen RBCs based on linear elastic thin-shell theory. Wijaya et al. (Reference Wijaya, Mohapatra, Sepehrirahnama and Lim2016) improved the numerical scheme and established a more efficient boundary element model since they only discretized the cell surface and omitted the acoustic domain discretization. Silva et al. (Reference Silva, Tian, Franklin, Wang, Han, Mann and Drinkwater2019) theoretically analysed the acoustic deformation of lipid vesicles in the regime of long wavelengths and small deformations. However, the acoustic radiation force formulations used in these studies are still based on the assumption of an ideal inviscid fluid.
To address the issue of viscous fluid, direct numerical simulation can theoretically be used to calculate the transient Navier–Stokes equations. However, it is essential to take into account two different time scales, i.e. acoustic oscillations and fluid dynamics, which are very different (Chini, Malecha & Dreeben Reference Chini, Malecha and Dreeben2014; Karlsen, Augustsson & Bruus Reference Karlsen, Augustsson and Bruus2016; Guglietta et al. Reference Guglietta, Behr, Biferale, Falcucci and Sbragaglia2020). This requires extremely high time resolution, resulting in expensive computational costs (Dubrovski, Friend & Manor Reference Dubrovski, Friend and Manor2023). Therefore, direct numerical simulation methods are not commonly used in practice. Typically, research in this field is based on acoustic perturbation theory (Bruus Reference Bruus2012; Rajendran et al. Reference Rajendran, Jayakumar, Azharudeen and Subramani2022), which decomposes the fluid variables into first-order and second-order components, achieving an explicit separation of time scales. This method simplifies the problem into two smaller sets of equations, where the fast time scale acoustic wave propagation equations are solved in the frequency domain, and the slow time scale fluid dynamics equations are solved in the time domain. It provides a good approximation of the flow field behaviour while maintaining reasonable computational cost. Currently, most of the studies applying acoustic perturbation theory are conducted within a fully Eulerian framework (Hahn, Lamprecht & Dual Reference Hahn, Lamprecht and Dual2016; Cui, Bhosale & Gazzola Reference Cui, Bhosale and Gazzola2024). Nama, Huang & Costanzo (Reference Nama, Huang and Costanzo2017) employed a generalized Lagrangian framework to analyse acoustic streaming. They directly express second-order problems using Lagrangian averaged variables, eliminating the need for postprocessing of second-order particle velocities related to the Eulerian velocity field. This framework is more suitable for handling flow problems involving large deformations or motion. Applying the Lagrangian acoustic perturbation method, we have established a numerical model to study the deformation of swollen spherical RBCs in acoustic streaming (Liu & Xin Reference Liu and Xin2023a ). The calculated steady-state deformation results were compared with Mishra et al.’s experimental data (Mishra et al. Reference Mishra, Hill and Glynne-Jones2014) to verify the computational model. Most of the related research objects focus on initially spherical or ellipsoidal cells (Liu & Xin Reference Liu and Xin2023b , Reference Liu and Xin2024). However, the initial shape of cells significantly affects the propagation of acoustic waves, the deformation patterns of cells and the distribution of the flow field. Under normal physiological conditions, RBCs exhibit an initial biconcave shape. It can be said that the acoustic deformation of biconcave RBCs differs significantly from that of spherical RBCs. To address the challenge of cell deformation with more complex shapes, this paper investigates the large deformation of biconcave RBCs under the action of acoustic streaming.
Unlike our previous work, which only considered swollen spherical cells, this work establishes for the first time a fully dynamic deformation model of biconcave RBCs in a single standing wave acoustic field, and specifically develops a curvature adaptive mesh refinement algorithm to address the error issues in fluid–structure interaction calculations caused by the nonlinear response of biconcave boundaries (such as curvature transitions). Through a high-precision fluid–structure interaction algorithm, this work has predicted and discovered the snap-through instability phenomenon of biconcave RBCs driven by ultrasound. Moreover, this study deeply analyses the effects of acoustic pressure, fluid viscosity, membrane shear modulus and bending stiffness on RBC deformation, revealing the transition from a biconcave shape to an oblate disk under high acoustic pressure and the occurrence of snap-through instability. These findings provide theoretical support for the application of acoustofluidic technology in the detection of RBC-related diseases. The rest of this paper is organized as follows. In § 2, the acoustic governing equations and the mechanical constitutive equations of the biconcave RBC membrane are provided. In § 3, the finite element model and numerical schemes used to solve the governing equations are described. In § 4, mesh-convergence analysis is performed as the first step to ensure the accuracy of the numerical model. Then, the numerical results of acoustic streaming pattern around the deformed cells and the distribution of fluid stress on the membrane surface are presented. Finally, the physical factors affecting the deformation of biconcave RBCs in acoustic streaming are discussed are discussed. In § 5, the conclusions are summarized.
2. Theoretical model
As shown in figure 1, the mechanical deformation of biconcave RBCs subjected to ultrasonic standing waves in a microfluidic channel is investigated. In figure 1(a), the channel, which is filled with a water-based buffer solution, features an inlet and an outlet for fluid flow. A piezoelectric transducer is installed at the bottom of the microchannel, which excites a half-wavelength ultrasonic standing wave with reflection layer, establishing a pressure node at the centre of the channel. The RBCs are effectively trapped at the acoustic pressure node by acoustic radiation forces. Figure 1(b) depicts the shape transition of the RBCs. Initially, the RBCs maintain their characteristic biconcave shape, illustrated in figure 1(bi). The stationary shape of a normal RBC without any external force is an axisymmetric biconcave shape, which is represented in the cylindrical coordinate system as
\begin{equation} z=\sqrt{1-\left(\frac{r}{R_{0}}\right)^{2}}\left[0.207+2.003\left(\frac{r}{R_{0}}\right)^{2}-1.123\left(\frac{r}{R_{0}}\right)^{4}\right] \end{equation}
where
$R_{0}$
denotes the maximum transverse radius of the cell. As the ultrasonic standing wave acts on the RBCs, a clear transformation occurs, as shown in figure 1(b ii), where the RBCs exhibit deformation. The RBC structure consists of a viscoelastic membrane enclosing a cytoplasm formed from a viscous haemoglobin solution. The membrane consists of a lipid bilayer, which has bending stiffness and area incompressibility to resist mechanical deformation, and a protein cytoskeleton, which provides in-plane elasticity. Both the cytoplasm and the water-based buffer solution are considered as homogeneous viscous fluids.
(a) Schematic of biconcave RBCs subjected to ultrasonic standing waves in a microfluidic channel. (b) Initially biconcave RBC captured at the acoustic pressure node and elliptically deformed RBC captured at the acoustic pressure node.
In order to theoretically describe the mechanical deformation of the biconcave RBC subjected to ultrasonic standing waves, a time scale separation method is employed. Here, the fluid response under ultrasonic excitation can be typically considered as a superposition of two components: fast time scale acoustic oscillations (typically
$0.1\,\unicode{x03BC} \text{s}$
) and slow time scale mean dynamics (typically 0.1 s). Typically, there are two coordinate frames that can be applied in the field of fluid dynamics, namely the Eulerian frame and the Lagrangian frame. The problem of deformation dynamics of biconcave RBCs in a viscous fluid driven by an ultrasonic standing wave is investigated here, which involves the interaction of the mean flow with the cell membrane and the deformation dynamics process of the cell membrane, which is a dynamic boundary acoustic problem. If described using the Eulerian framework as a coordinate system, it is difficult to deal with the nonlinear and time-varying nature of these interactions, although the mathematical form of the equations is simple. The Lagrangian framework can be used to track dynamic boundaries and fluid particles, and is a more effective coordinate system description method to deal with dynamic boundary problems. Therefore, for the coupled acoustic–fluid-structural dynamics problem in this study, an acoustic perturbation method based on the generalized Lagrangian framework is adopted here.
Specifically, the current configuration
$B_{t}$
is perturbed by acoustic oscillations, and the Eulerian framework requires a Taylor series expansion to describe the boundary conditions of the mean dynamics, whereas the generalized Lagrangian framework is adopted in this study, which can effectively avoid the truncated approximation of the Taylor series. As shown in figure 2, the configuration before acoustic excitation is represented by the reference configuration
$B_{0}$
. The actual motion of the particles in the system is represented by the current configuration
$B_{t}$
. The introduced mean configuration
$B$
is an intermediate configuration between the reference configuration and the current configuration. The actual motion of particles in the current configuration can be regarded as the superposition of mean motion in the mean configuration and the acoustic oscillation. In the Eulerian configuration, all fluid quantities are originally understood as functions of time
$t$
and position
${y}$
in the current configuration, which does not explicitly include the mean motion. Now, in the Lagrangian configuration, all flow variables are remapped to the mean configuration
$B$
, and these variables are understood as functions of time t and position
${x}$
in the mean configuration, independent of acoustic oscillations. The treatment of second-order time-averaged motion in the Lagrangian framework is more efficient compared with the Eulerian framework.
Generalized Lagrangian configurations in the theoretical model. Here
$B_{0}$
indicates the reference configuration before acoustic actuation,
$B$
indicates the mean configuration of particle mean motion and
$B_{t}$
indicates the current configuration of particle actual motion. The actual motion of fluid and membrane particles can be considered as a combination of the mean motion on the slow time scale and the acoustic oscillation on the fast time scale.
The mapping relationship between the current configuration
$B_{t}$
and the mean configuration
$B$
is
$\boldsymbol{y}=\boldsymbol{x}+\boldsymbol{\xi }$
. According to this relationship, the Navier–Stokes equations in the current configuration
$B_{t}$
can be mapped to the mean configuration
$B$
. The Navier–Stokes equations in the mean configuration
$B$
can be reformulated as
\begin{align} \partial _{t}\rho +\boldsymbol{F}_{\xi }^{-T}\boldsymbol{\cdot }\boldsymbol{\nabla }\rho \boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{v}_{\xi })+\rho \boldsymbol{F}_{\xi }^{-T}\colon \boldsymbol{\nabla }\boldsymbol{v}&=0 , \end{align}
\begin{align} J_{\xi }\rho \big[\partial _{t}\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{F}_{\xi }^{-1}\boldsymbol{\cdot }(\boldsymbol{v}-\boldsymbol{v}_{\xi })\big]-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}&=\boldsymbol{0} , \end{align}
where
$\rho$
is the fluid density,
$\boldsymbol{v}$
is the fluid velocity,
$\boldsymbol{F}_{\xi }=\boldsymbol{\nabla }\boldsymbol{y}=\boldsymbol{i}+\boldsymbol{\nabla }\boldsymbol{\xi }$
represents the gradient of acoustic oscillation displacement,
$J_{\xi }=\det (\boldsymbol{F}_{\xi })$
is the Jacobian determinant of the acoustic displacement gradient,
$\boldsymbol{v}_{\xi }=\partial _{t}\boldsymbol{\xi }$
is the time derivative of the acoustic displacement and
$\boldsymbol{P}=J_{\xi }\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{F}_{\xi }^{-T}$
is the Piola–Kirchhoff stress defined in the mean configuration, where
$\boldsymbol{\sigma }$
is the Cauchy stress defined in the current configuration. The Cauchy stress of a viscous compressible fluid is given by
\begin{equation} \boldsymbol{\sigma }=-p\boldsymbol{i}+\mu \big(\boldsymbol{\nabla} _{y}\boldsymbol{v}+\boldsymbol{\nabla} _{y}\boldsymbol{v}^{T}\big)+\left(-\frac{2}{3}\mu +\mu _{b}\right)\left(\boldsymbol{\nabla} _{y}\boldsymbol{\cdot }\boldsymbol{v}\right)\boldsymbol{i} ,\end{equation}
where
$\boldsymbol{\nabla} _{y}=\boldsymbol{F}_{\xi }^{-1}\boldsymbol{\cdot }\boldsymbol{\nabla }$
represents the gradient operator in the current configuration,
$\boldsymbol{i}$
is the unit tensor,
$\mu$
and
$\mu _{b}$
are the shear viscosity and the bulk viscosity, respectively. The fluid pressure
$p$
has the relationship with the density as
\begin{equation} p-p_{0}=c_{0}^{2}\left(\rho -\rho _{0}\right)\! ,\end{equation}
where
$p_{0} , \rho _{0}$
and
$c_{0}$
are pressure, density and sound speed of the stationary fluid, respectively. The acoustic perturbation method is applied to obtain the governing equations for the two time scales, in which the pressure and velocity of the fluid are expanded to the second order:
$p=p_{0}+p_{1}+p_{2}+ {\cdots} , \boldsymbol{v}=\boldsymbol{v}_{1}+\boldsymbol{v}_{2}+{\cdots}$
, where the subscripts denote the respective orders. Correspondingly, the governing equations of (2.2) and (2.3) can be divided into the first-order equations governing the ultrasonic propagation on the fast time scale and the second-order equations governing the time-averaged dynamics on the slow time scale.
According to the acoustic perturbation method, substituting the perturbed fields into the governing equation, the first-order equations governing the propagation of acoustic waves can be expressed as
\begin{align}&\quad\,\, \partial _{t}p_{1}+c_{0}^{2}\rho _{0}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{1}=0 \quad {\rm and}\quad \rho _{0}\partial _{t}\boldsymbol{v}_{1}-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_{\boldsymbol{1}}=\boldsymbol{0} , \end{align}
\begin{align}& \boldsymbol{P}_{1}=-p_{1}\boldsymbol{i}+\mu \big(\boldsymbol{\nabla }\boldsymbol{v}_{1}+\boldsymbol{\nabla }\boldsymbol{v}_{1}^{T}\big)+\left(-\frac{2}{3}\mu +\mu _{b}\right)\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{1}\right)\boldsymbol{i} , \end{align}
where all variables with subscript 1 are the first-order variables,
$p_{1}$
is the first-order acoustic pressure,
$\boldsymbol{v}_{1}$
is the first-order acoustic particle velocity and
$\boldsymbol{P}_{1}$
the first-order Piola–Kirchhoff stress. Utilizing the relationship
$\boldsymbol{v}_{1}=\partial _{t}\boldsymbol{\xi }=i\omega \boldsymbol{\xi }$
with
$\boldsymbol{\xi }$
being the acoustic particle displacement, the first-order equation (2.6) can be expressed in terms of the acoustic particle displacement field
$\boldsymbol{\xi }$
as
\begin{equation} \rho _{0}\left(\mathrm{i}\omega \right)^{2}\boldsymbol{\xi }-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_{1}\left(\boldsymbol{\xi }\right)=\boldsymbol{0} \end{equation}
with the first-order Piola–Kirchhoff stress as
\begin{equation} \boldsymbol{P}_{1}\left(\boldsymbol{\xi }\right)=\left[c_{0}^{2}\rho _{0}+\mathrm{i}\omega \left(-\frac{2}{3}\mu +\mu _{b}\right)\right]\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi }\right)\boldsymbol{i}+\mathrm{i}\omega \mu \big(\boldsymbol{\nabla }\boldsymbol{\xi }+\boldsymbol{\nabla }\boldsymbol{\xi }^{T}\big). \end{equation}
To describe the mean deformation dynamics of the fluid–cell system on the slow time scale, the time-averaging operator is introduced for time averaging processing. The time averaging operation is performed to filter out periodic variations in the fluid flow on the fast time scale, while preserving the nonlinear mean behaviour on the slow time scale over one oscillation period of the acoustic wave
$T=2\pi /\omega$
. For any acoustic quantity
$\phi (\boldsymbol{x},t)$
, the time-averaging operator is
\begin{equation} \left\langle \phi \left(\boldsymbol{x},t\right)\right\rangle =\frac{1}{T}\int _{t}^{t+T}\phi \left(\boldsymbol{x},\tau \right){\rm d}\tau .\end{equation}
The time-averaged second-order equations governing the mean fluid dynamics can be given by
\begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{v}_{2}\rangle=0\quad{\rm and}\quad\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{P}_{2}\rangle=0 \,{\rm in}\,\mathrm{\mathcal{B}}, \end{equation}
with the time-averaged second-order Piola–Kirchhoff stress
$\langle \boldsymbol{P}_{2}\rangle$
given by
\begin{align} \left\langle \boldsymbol{P}_{2}\right\rangle &=-\left\langle p_{2}\right\rangle \boldsymbol{i}+\mu \big(\langle \boldsymbol{\nabla }\boldsymbol{v}_{2}\rangle +\big\langle \boldsymbol{\nabla }\boldsymbol{v}_{2}^{T}\big\rangle \big)-\mu \big(\langle \boldsymbol{\nabla }\boldsymbol{v}_{1}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\xi }\rangle +\big\langle \boldsymbol{\nabla }\boldsymbol{\xi }^{T}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{1}^{T}\big\rangle \big)\nonumber \\&\quad -\left(-\frac{2}{3}\mu +\mu _{b}\right)\big\langle \boldsymbol{\nabla }\boldsymbol{\xi }^{T}\colon \boldsymbol{\nabla }\boldsymbol{v}_{1}\big\rangle \boldsymbol{i}+\big\langle \boldsymbol{P}_{1}\boldsymbol{\cdot }\big[(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi })\boldsymbol{i}-\boldsymbol{\nabla }\boldsymbol{\xi }^{T}\big]\big\rangle . \end{align}
In (2.11), the divergence-free velocity field implies that the volume remains constant in the mean deformation dynamics of the fluid–cell system. The governing equations (2.11) do not contain time derivatives, indicating that this second-order problem is time-independent, where the second-order velocity field is an acoustic streaming field driven by acoustic viscous dissipation, and where the acoustic streaming field is intrinsically stable. However, when the acoustic streaming is coupled to the motion and deformation of the cell, the boundary conditions at the fluid–membrane interface change with time, making it a time-dependent problem. In other words, the components of the acoustic streaming field
$\langle \boldsymbol{v}_{2}\rangle$
include the flow due to the transient deformation of the cell membrane, which is time-dependent.
Since the considered biconcave RBCs lack nuclei and organelles, but contain a homogeneous haemoglobin solution, that is the cytoplasm. From mechanical perspective, RBC can be modelled as a viscous droplet encapsulated by a viscoelastic membrane. The internal cytoplasm and external solution are modelled as a Newtonian fluid with viscosities
$\mu _{{\textit{in}}}$
and
$\mu _{{out}}$
, respectively. The cell membrane is considered as a two-dimensional (2-D) surface embedded in three-dimensional (3-D) fluids because its thickness is extremely thin relative to the cell diameter. To describe the mechanical model of the cell membrane, the Kirchhoff–Love theory of hyperelastic thin shells (Duong, Roohbakhshan & Sauer Reference Duong, Roohbakhshan and Sauer2017) is applied here, which can describe the large deformation and rotation of the cell membrane. To apply the hyperelastic thin-shell theory, the curvature tensor of the mean deformed surface
$\mathrm{\mathcal{S}}$
and the mean curvature are introduced,
\begin{align} \boldsymbol{b}&=-\boldsymbol{\nabla} _{\!s}\boldsymbol{n} , \end{align}
\begin{align} h&=\frac{1}{2}\boldsymbol{b}\colon \boldsymbol{i}_{\!s} , \end{align}
where
$\boldsymbol{\nabla} _{\!s}=\boldsymbol{i}_{\!s}\boldsymbol{\cdot }\boldsymbol{\nabla }$
is the surface gradient operator, and
$\boldsymbol{i}_{\!s}=\boldsymbol{i}-\boldsymbol{n}\otimes \boldsymbol{n}$
is the surface identity tensor with
$\boldsymbol{n}$
is the outward unit normal of the mean deformed surface
$\mathrm{\mathcal{S}}$
, and
$\boldsymbol{\nabla }$
is the gradient operator. As shown in figure 2, based on nonlinear continuum mechanics, let
$\boldsymbol{X}_{\!s}$
represent the particle position in the reference configuration (i.e. the initial shape of the RBC) and
$\boldsymbol{x}_{\!s}(t)$
represent the particle positions in the mean configuration of the cell membrane at time
$t$
. The deformation map between
$\mathrm{\mathcal{S}}_{0}$
and
$\mathrm{\mathcal{S}}$
is characterized by the surface deformation gradient
$\boldsymbol{F}_{\!s}$
obtained from
$\boldsymbol{F}_{\!s}=\boldsymbol{\nabla} _{\!{s_{0}}}\boldsymbol{\varphi }=\mathrm{d}\boldsymbol{x}_{\!s}/\mathrm{d}\boldsymbol{X}_{\!s}\boldsymbol{\cdot }\boldsymbol{i}_{\!s}$
. The surface left Cauchy–Green deformation tensor
$\boldsymbol{B}_{\!s}=\boldsymbol{F}_{\!s}\boldsymbol{\cdot }\boldsymbol{F}_{\!s}^{T}$
measures the in-plane deformation of the cell membrane, and its two invariants are
\begin{equation} I_{1}=\mathrm{tr}(\boldsymbol{B}_{\!s})-2\,\,{\rm{and}}\,\,I_{2}=\det \left(\boldsymbol{B}_{\!s}\right)-1=J_{\!s}^{2}-1 .\end{equation}
Here,
$J_{\!s}=\sqrt{\det (\boldsymbol{B}_{\!s})}$
represents the ratio of the mean deformed surface area to the undeformed surface area. The equation governing the temporal evolution of
$\boldsymbol{B}_{\!s}$
can be obtained by taking the material derivative for the defining equation of
$\boldsymbol{B}_{\!s}$
,
\begin{equation} \frac{\mathrm{d}\boldsymbol{B}_{\!s}}{\mathrm{d}t}=\boldsymbol{B}_{\!s}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!s}\boldsymbol{v}_{\!s}+\boldsymbol{\nabla} _{\!s}\boldsymbol{v}_{\!s}^{T}\boldsymbol{\cdot }\boldsymbol{B}_{\!s} ,\end{equation}
where
$\boldsymbol{v}_{\!s}=\mathrm{d}\boldsymbol{x}_{\!s}/\mathrm{d}t$
is the velocity of the particle on the membrane.
Assuming that the viscoelastic behaviour of the membrane follows the Kelvin–Voigt model, the total membrane tension is the sum of the elastic and viscous tensions,
\begin{equation} \boldsymbol{\varSigma }_{\!s}=\boldsymbol{\varSigma }_{\!s}^{e}+\boldsymbol{\varSigma }_{\!s}^{b}+\boldsymbol{\varSigma }_{\!s}^{\nu } .\end{equation}
From a mechanical point of view, the structure of the RBC membrane, which consists mainly of a lipid bilayer and an embedded protein cytoskeleton, is resistant to shear deformation, area expansion and bending. In the case of shear and area expansion deformations, the elastic tension
$\boldsymbol{\varSigma }_{\!s}^{e}$
can be calculated using the strain-hardening Skalak’s law (Skalak et al. Reference Skalak, Tozeren, Zarda and Chien1973), and the in-plane strain energy function is given by
\begin{equation} W_{e}=\frac{G_{\!s}}{4}\big(I_{1}^{2}+2I_{1}-2I_{2}+CI_{2}^{2}\big),\end{equation}
where
$G_{\!s}$
is the surface shear modulus, and C represents the elastic constant related to membrane area expansion. The larger the value of C, the more restricted the area expansion. Here, C is set to 100 to ensure surface incompressibility. In this way, the total change in surface area is limited to a negligible amount, ensuring the stability of numerical calculations. The elastic tension
$\boldsymbol{\varSigma }_{\!s}^{e}$
associated with shear deformation and area expansion can be derived from the strain energy function
\begin{equation} \boldsymbol{\varSigma }_{\!s}^{e}=\frac{G_{\!s}}{J_{\!s}}\big[\left(I_{1}+1\right)\boldsymbol{B}_{\!s}+J_{\!s}^{2}\left(-1+CI_{2}\right)\boldsymbol{i}_{\!s}\big].\end{equation}
The bending elasticity is modelled according to the Helfrich’s bending energy formulation (Helfrich Reference Helfrich1973)
\begin{equation} W_{b}=\frac{1}{2}E_{b}J_{\!s}\left(h-h_{0}\right)^{2} ,\end{equation}
where
$E_{b}$
is bending stiffness,
$h$
is the mean curvature and
$h_{0}$
is the spontaneous curvature in the Helfrich curvature elastic energy. The membrane spontaneously bends to one side, and the spontaneous curvature is used to characterize the asymmetry between the two sides of the membrane. Given the bending energy density function
$W_{b}$
, the constitutive equations can be given by
$\boldsymbol{\varSigma }_{\!s}^{b}=({2}/{J_{\!s}})\boldsymbol{F}_{\!s}\boldsymbol{\cdot }({\partial W_{b}}/{\partial \boldsymbol{C}_{\!s}})\boldsymbol{\cdot }\boldsymbol{F}_{\!s}^{T}$
. The tension
$\sum _{s}^{b}$
resisting bending deformation is calculated as
\begin{equation} \boldsymbol{\varSigma }_{\!s}^{b}=2E_{b}\big[(h-h_{0})^{2}\boldsymbol{i}_{\!s}-(h-h_{0})\boldsymbol{b}+\boldsymbol{n}\otimes \boldsymbol{\nabla} _{\!s}(h-h_{0})\big].\end{equation}
Furthermore, the membrane viscous effect can be modelled by considering the linear viscous constitutive relation, and the membrane tension due to the membrane viscosity is (Barthès-Biesel Reference Barthès-Biesel2016)
\begin{equation} \boldsymbol{\varSigma }_{\!s}^{\nu }=\mu _{s}\big[\boldsymbol{i}_{\!s}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!s}\boldsymbol{v}_{\!s}+\boldsymbol{\nabla} _{\!s}\boldsymbol{v}_{\!s}^{T}\boldsymbol{\cdot }\boldsymbol{i}_{\!s}-\left(\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{v}_{\!s}\right)\boldsymbol{i}_{\!s}\big]+\mu _{\!s}^{b}\left(\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{v}_{\!s}\right)\boldsymbol{i}_{\!s} ,\end{equation}
where
$\mu _{\!s}$
and
$\mu _{\!s}^{b}$
are the membrane shear and expansion viscosities, respectively.
Considering the mechanical coupling relationship between the RBC membrane and the surrounding fluid, the first is quasistatic mechanical equilibrium: due to the negligible inertia of membrane, the membrane viscoelastic tension generated by the membrane deformation is in equilibrium with fluid force
\begin{align}&\quad\,\, \boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{\varSigma }_{\!s}+\Delta \boldsymbol{f}=\boldsymbol{0} , \end{align}
\begin{align}& \Delta \boldsymbol{f}=\big(\big\langle \boldsymbol{P}_{2}^{\textit{ex}}\big\rangle -\big\langle \boldsymbol{P}_{2}^{\textit{in}}\big\rangle \big)\boldsymbol{\cdot }\boldsymbol{n} , \end{align}
where
$\Delta \boldsymbol{f}$
is the jump in fluid stress at the interface (force per unit area of the mean deformed surface),
$\boldsymbol{P}_{2}^{\textit{in}}$
and
$\boldsymbol{P}_{2}^{\textit{ex}}$
are the fluid stress on the inner and outer sides of the cell membrane and
$\boldsymbol{\varSigma }_{\!s}$
is the membrane tension in (2.17). Second, in the mean dynamics, there is a no-slip boundary condition at the interface: there is no relative motion between the second-order mean fluid flow and the cell membrane at the interface
\begin{equation} \boldsymbol{v}_{\!s}=\big\langle \boldsymbol{v}_{2}^{\textit{ex}}\big\rangle =\big\langle \boldsymbol{v}_{2}^{\textit{in}}\big\rangle ,\end{equation}
where
$\boldsymbol{v}_{\!s}$
is membrane velocity,
$\langle \boldsymbol{v}_{2}^{\textit{in}}\rangle$
and
$\langle \boldsymbol{v}_{2}^{\textit{ex}}\rangle$
are the mean fluid velocity on the inner and outer sides of the cell membrane. Based on the above mechanical coupling relationship, the second-order fluid mean dynamics and the membrane deformation dynamics are effectively combined.
3. Finite element model
The governing equations describing the acoustically induced deformation of the RBC are solved numerically using the finite element method. Since both the RBC and the incident acoustic wave are symmetric with respect to the polar axis of cell, a 2-D axisymmetric finite element model is established to reduce the computational cost. As shown in figure 3, the geometry of the model in the mean configuration consists of the inner and outer fluids of the cell, which are separated by the cell membrane. The red region represents the cytoplasm of the RBC, and the green region represents the buffer solution. Since the cell diameter is much smaller than the microchannel size, the modelling of the microchannel wall structure is neglected and only the physical field at the acoustic pressure node at the centre of the channel is focused on, as shown in figure 3(a). To establish the standing wave acoustic field around the cell, the acoustic particle displacement
$\boldsymbol{\xi }$
is applied on the boundary of the entire computational domain
\begin{equation} \boldsymbol{\xi }=\frac{p_{0}}{2\omega \rho _{0}^{\textit{ex}}c_{0}^{\textit{ex}}}\left[e^{i{k^{\textit{ex}}}z}+e^{-i{k^{\textit{ex}}}z}\right]e^{i\omega t}\boldsymbol{e}_{z} \end{equation}
where
$p_{0}$
is the acoustic pressure amplitude,
$k=\omega /c_{0}^{\textit{ex}}$
is the wavenumber in the outer fluid. The parameters used in the calculations are detailed in table 1.
Finite element model of a biconcave RBC in ultrasonic standing waves: (a) 2-D axisymmetric computational domain; (b) typical mesh; (c) a zoom into RBC.
Parameters used in the calculation. The bulk viscosities of the fluids and cell membrane are set to zero. The choice of these parameters refers to the parameters of RBCs (Cordasco, Yazdani & Bagchi Reference Cordasco, Yazdani and Bagchi2014; Mishra et al. Reference Mishra, Hill and Glynne-Jones2014; Guglietta et al. Reference Guglietta, Behr, Biferale, Falcucci and Sbragaglia2020).
Here, the Boltzmann constant
$k_{b}=1.38\times 10^{-23}\,\mathrm{J}\,\mathrm{K}^{-1}$
and the temperature
$T = 300\,\rm K$
.
Numerical simulations are performed using the finite element method to convert the theoretical model equations into the weak form of the finite element model. To model the acoustic field, the first-order acoustic governing equation is multiplied by the test function
$\delta \boldsymbol{\xi }$
of acoustic particle displacement
$\boldsymbol{\xi }$
and integrated over the entire computational domain, and using the divergence theorem, it is obtained that
\begin{equation} \int _{\mathcal{B}}\big[\rho _{0}\left(\mathrm{i}\omega \right)^{2}\boldsymbol{\xi }\boldsymbol{\cdot }\delta \boldsymbol{\xi }+\boldsymbol{P}_{1}\colon \boldsymbol{\nabla }\delta \boldsymbol{\xi }\big]\mathrm{d}v-\int _{\mathrm{\mathcal{S}}}\big(\boldsymbol{P}_{1}^{\textit{ex}}-\boldsymbol{P}_{1}^{\textit{in}}\big)\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\delta \boldsymbol{\xi }\mathrm{d}a=0. \end{equation}
The integral contribution at the interface is zero due to the first-order stress continuity condition, so (3.2) becomes
\begin{equation} \int _{\mathcal{B}}\big[\rho _{0}\left(\mathrm{i}\omega \right)^{2}\boldsymbol{\xi }\boldsymbol{\cdot }\delta \boldsymbol{\xi }+\boldsymbol{P}_{1}\colon \boldsymbol{\nabla }\delta \boldsymbol{\xi }\big]\mathrm{d}v=0 .\end{equation}
The acoustic particle displacement continuity condition is ensured by the continuity between elements in the finite element method.
To simulate the mean fluid dynamics, the second-order mean fluid dynamics (2.11) is multiplied by the test function pair
$(\delta \langle p_{2}\rangle ,\delta \langle \boldsymbol{v}_{2}\rangle )$
of mean flow pressure and velocity
$(\langle p_{2}\rangle ,\langle \boldsymbol{v}_{2}\rangle )$
and integrated over the entire fluid domain, and using the divergence theorem
\begin{align}&\qquad\qquad\qquad\qquad \int _{\mathcal{B}}\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\left\langle \boldsymbol{v}_{2}\right\rangle \right)\boldsymbol{\cdot }\delta \left\langle p_{2}\right\rangle \mathrm{d}v=0, \end{align}
\begin{align}& \int _{\mathcal{B}}-\left\langle \boldsymbol{P}_{2}\right\rangle \boldsymbol{\cdot }\boldsymbol{\nabla }\delta \left\langle \boldsymbol{v}_{2}\right\rangle \mathrm{d}v-\int _{\mathrm{\mathcal{S}}}\big(\big\langle \boldsymbol{P}_{2}^{\textit{ex}}\big\rangle -\big\langle \boldsymbol{P}_{2}^{\textit{in}}\big\rangle \big)\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\delta \left\langle \boldsymbol{v}_{2}\right\rangle \mathrm{d}a=0 . \end{align}
In the numerical modelling of membrane deformation mechanics problem, the calculation of the membrane tension
$\boldsymbol{\varSigma }_{\!s}$
involves the curvature tensor
$\boldsymbol{b}$
and the surface left Cauchy–Green deformation tensor
$\boldsymbol{B}_{\!s}$
. In this part,
$\boldsymbol{b}$
and
$\boldsymbol{B}_{\!s}$
are calculated according to the weak form of the relevant (2.13) and (2.16). Similarly, the weak formulations can be obtained by multiplying (2.13) and (2.16) with the respective test function
$\delta \boldsymbol{b}$
and
$\delta \boldsymbol{B}_{\!s}$
, integrating over the interface
$\mathcal{S}$
, and using the surface divergence theorem,
\begin{align} \int _{\!s}\left[\boldsymbol{b}\colon \delta \boldsymbol{b}-\left(\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\delta \boldsymbol{b}\right)\boldsymbol{\cdot }\boldsymbol{n}-\left(\boldsymbol{\nabla} _{\!s}\boldsymbol{\cdot }\boldsymbol{b}\right)\boldsymbol{n}\otimes \boldsymbol{n}\colon \delta \boldsymbol{b}\right]\mathrm{d}a&=0 , \end{align}
\begin{align} \int _{{s}}\left(\frac{{\rm d}\boldsymbol{B}_{\!s}}{{\rm d}{t}}-\boldsymbol{B}_{\!s}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!s}\left\langle \boldsymbol{v}_{2}\right\rangle -\boldsymbol{\nabla} _{\!s}\left\langle \boldsymbol{v}_{2}\right\rangle ^{T}\boldsymbol{\cdot }\boldsymbol{B}_{\!s}\right)\colon \delta \boldsymbol{B}_{\!s}\mathrm{d}a&=0 . \end{align}
Following the no-slip condition of (2.25), the membrane velocity
$\boldsymbol{v}_{\!s}$
can be replaced by
$\langle \boldsymbol{v}_{2}\rangle$
in (3.7). According to (2.23) and (2.24), (30) is transformed into a weak form that directly represents the fluid–membrane coupling problem
\begin{equation} \int _{\beta }-\left\langle \boldsymbol{P}_{2}\right\rangle \colon \boldsymbol{\nabla }\delta \left\langle \boldsymbol{v}_{2}\right\rangle \mathrm{d}v+\int _{\!s}\boldsymbol{\varSigma }_{\!s}\colon \boldsymbol{\nabla} _{\!s}\delta \left\langle \boldsymbol{v}_{2}\right\rangle \mathrm{d}a=0 .\end{equation}
Based on the moving mesh method, large membrane deformations can be tracked in the deformation domain. In the numerical discretization process, the solutions of the acoustic particle displacement
$\boldsymbol{\xi }$
are approximated with fourth-order Lagrange elements, and the solutions of the curvature tensor
$\boldsymbol{b}$
and the surface left Cauchy–Green deformation tensor
$\boldsymbol{B}_{\!s}$
are approximated with third-order Lagrange elements. The solutions of the pressure
$\langle p_{2}\rangle$
and the velocity
$\langle \boldsymbol{v}_{2}\rangle$
are approximated with third-order and second-order Lagrange elements, respectively, to ensure full conservation properties and high accuracy. In the 2-D axisymmetric model, the entire fluid domain is meshed using the ‘free triangle’ mesh, which is suitable for discretizing complex geometric shapes. To address the error issues caused by the nonlinear response of biconcave boundaries (such as curvature transitions) in fluid–structure interaction calculations, a curvature-adaptive mesh refinement algorithm is specifically developed (see the Supplementary materials for details). For the time stepping strategy, an implicit time stepping method, the backward difference formulation, is used and
$5\times 10^{-4}$
s is chosen as a typical time step to ensure numerical stability.
4. Results and discussion
To validate the reliability of the fluid–membrane coupling algorithm based on the finite element method in this study, the present finite element model of cell deformation in the acoustic field is degraded to calculate the transient deformation of an RBC in an extensional flow and compared with previous results (Sudip Das, Deshmukh & Thaokar Reference Das, Deshmukh and Thaokar2019). The validation model only requires setting the acoustic input to zero and applying an extensional flow field, while keeping all other simulation settings unchanged. As shown in figure 4, the deformation of a biconcave RBC in a polar extensional flow is simulated. In the 2-D axisymmetric coordinate system used in our modelling, the extensional flow is given by
\begin{equation} \boldsymbol{v}=e\left[\begin{array}{lll} -1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 2 \end{array}\right]\left[\begin{array}{l} r\\ 0\\ z \end{array}\right] \end{equation}
where e is the axisymmetric uniaxial extensional strain rate. The symmetry axis of the RBC coincides with the z-axis, so its centre of mass is located at the origin (see figure 4). The flow capillary number
$Ca_{f}$
is a dimensionless parameter used to determine the deformation dynamics of the RBCs,
\begin{equation} Ca_{f}=\frac{\mu er_{0}}{E_{\!s}} ,\end{equation}
where
$\mu$
is the fluid viscosity,
$r_{0}$
is the maximum transverse radius of the cell and
$E_{s}=2G_{\!s}(1+2C)/(1+C)$
is the surface Young’s modulus, which measures the resistance of the RBC to deformation (Barthès-Biesel et al. Reference Barthès-Biesel, Diaz and Dhenin2002). In the verification simulation, the area expansion related parameter C is set to 10 to facilitate comparison with previous studies.
Schematic diagram of a biconcave RBC under extensional flow.
Figure 5 shows the deformation of biconcave RBCs under extensional flow along the polar direction at different fluid capillary numbers
$Ca_{f}$
, with figure 5(a–d) showing the steady-state deformation shapes of biconcave RBCs in the extensional flow. When the capillary number is small (e.g,
$Ca_{f}=0.01$
), it means that the strength of the extensional flow is small relative to the deformation resistance of the RBC. As shown in figure 5(a), the cell deformation is small and maintains a biconcave shape. Figure 5(b) shows that as the capillary number increases, the biconcave RBC is stretched along the polar direction, and the steady-state shape reaches a biconcave cylinder. As shown in figure 5(c), when the capillary number reaches the critical value (
$Ca_{f}=0.04$
), the biconcave cavity of the RBC suddenly opens outward, and the cell shape suddenly changes from biconcave to biconvex, i.e. a snap-through instability phenomenon occurs. When the capillary number is even higher (
$Ca_{f}=0.08$
), as shown in figure 5(d), the cell eventually evolves into a highly stretched prolate spheroid shape.
Steady shapes of deformed RBCs under extensional flow along the polar direction at different capillary numbers
$Ca_{f}$
: (a)
$Ca_{f}=0.01$
; (b)
$Ca_{f}=0.038$
; (c)
$Ca_{f}=0.04$
; (d)
$Ca_{f}=0.08$
. (e) Polar displacement of deformed RBCs at different capillary numbers. Red circles connected by dashed lines represent the numerical results of Sudip Das et al. (Reference Das, Deshmukh and Thaokar2019) based on the boundary integral method. Blue symbols connected by solid lines represent the present numerical results. Supplementary movies of our computational results are available in the Supplementary material. (Since the snap-through instability occurs near the capillary number
$Ca_{f}=0.04$
, the deformation process of
$Ca_{f}=0.038 , Ca_{f}=0.04 , Ca_{f}=0.042$
are specially selected to provide three animation demonstrations.)
In figure 5(e), the polar displacement of the deformed RBC is plotted as a function of the flow capillary number. The polar displacement, defined as
$(z-z_{0})/z_{0}$
, is a parameter to measure the deformation of the cell in the extensional flow, which represents the equilibrium displacement of the north pole, where
$z$
and
$z_{0}$
are the position of polar steady deformation and the position at the initial moment, respectively. Figure 5(e) shows that the biconcave RBC undergoes snap-through instability through the polar displacement, transforming from the biconcave disk to the prolate-spheroid. The shape evolution trend and polar displacement are consistent with the results obtained by Sudip Das et al. (Reference Das, Deshmukh and Thaokar2019) using the boundary integral method. This confirms the reliability of our RBC model based on the finite element method, indicating that the RBC model in this study can accurately simulate RBC deformation.
In addition to verifying the reliability of the fluid–membrane coupling algorithm based on the finite element method, to ensure the accuracy and efficiency of the calculations, the convergence verification of the spatial mesh density and time step of the computational model is conducted in the Appendix. It is found that when the number of cell membrane line elements is selected as
$N=100$
, the number of fluid elements inside and outside the cell membrane is selected as
$M=12\,157$
, and the time step is selected as
$t=5\times 10^{-4}\,\mathrm{s}$
, the results of the computational model has converged. Therefore, these spatial mesh densities and time steps are used in subsequent computational studies.
To investigate the effect of acoustic pressure input on cell deformation, a quantitative comparison of different acoustic pressures is performed. In previous studies, most research subjects were elliptical RBCs, whose deformation pattern could be measured using a single aspect ratio parameter. In contrast, this study investigates the deformation of biconcave RBCs, whose biconcave shape is more complex than that of elliptical cells, thus requiring the use of two parameters, transverse diameter and polar diameter, to measure the deformation of biconcave RBCs. Figures 6(a) and 6(b) show the time evolution of the transverse diameter and polar diameter of biconcave RBCs at different acoustic pressure inputs. As shown in figures 6(a) and 6(b), when the acoustic pressure input is minimum, i.e.
$p_{0}=0.58\,\mathrm{MPa}$
, the RBC deformation increases slowly and reaches a steady state in more than 0.5 s. When the acoustic pressure input is
$p_{0}=0.80\,\mathrm{MPa}$
, the final deformation values of RBCs increase. The deformation time decreases and reaches equilibrium in approximately 0.4 s. When the acoustic pressure input is maximum
$p_{0}=1.83\,\mathrm{MPa}$
, the transverse and polar diameters of the RBC elongate the fastest and largest, reaching equilibrium immediately in less than 0.05 s. Therefore, as the acoustic pressure input increases, the deformation of the RBC accelerates and the degree of deformation increases.
Time evolution of RBCs at different acoustic pressure inputs
$p_{0}=0.58\,\mathrm{MPa} , 0.80\,\mathrm{MPa}$
and
$1.83\,\mathrm{MPa}$
: (a) transverse diameter; (b) polar diameter.
The acoustic-induced deformation of the cell is driven by the jump of time-averaged second-order mean stress
$\boldsymbol{P}_{2}$
across the cell membrane. According to (2.12), the mean stress
$\langle \boldsymbol{P}_{2}\rangle$
driving the cell membrane deformation is closely related to the first-order acoustic velocity
$\boldsymbol{v}_{1}$
and second-order mean flow velocity
$\langle \boldsymbol{v}_{2}\rangle$
. To understand the transient dynamics of the RBC deformation behaviour and the surrounding flow, the 3-D shape of biconcave RBCs, the second-order mean flow pattern and the first-order acoustic velocity distribution during the deformation process under different acoustic pressure inputs
$p_{0}=0.58\,\mathrm{MPa}$
and
$p_{0}=1.83\,\mathrm{MPa}$
are plotted in figures 7 and 8, respectively.
(a)–(c) The 3-D shapes of biconcave RBC at three different time points: the initial state at
$t=0\,\mathrm{s}$
; the transient deformation at
$t=0.1\,\mathrm{s}$
; the steady deformation at
$t=1\,\mathrm{s}$
. (d)–(f) The corresponding second-order mean flow pattern and streamlines, with the black solid line outlining the RBC. (g)–(i) The corresponding r.m.s. first-order acoustic velocity distribution. The input acoustic pressure amplitude is
$p_{0}=0.58\,\mathrm{MPa}$
, and the time evolution of the transverse diameter and polar diameter during the deformation is shown in figures 5(a) and 5(b).
(a)–(c) The 3-D shapes of biconcave RBC at three different time points: the initial state at
$t=0\,\mathrm{s}$
; the transient deformation at
$t=0.01\,\mathrm{s}$
; the steady deformation at
$t=1\,\mathrm{s}$
. (d)–(f) The corresponding second-order mean flow pattern and streamlines, with the black solid line outlining the RBC. (g)–(i) The corresponding r.m.s. first-order acoustic velocity distribution. The input acoustic pressure amplitude is
$p_{0}=1.83\,\mathrm{MPa}$
, and the time evolution of the transverse diameter and polar diameter during the deformation is shown in figures 5(a) and 5(b).
First, the deformation shapes of biconcave RBCs under different acoustic pressure inputs are investigated. Comparing figure 7(a–c) with figure 8(a–c), it can be seen that at lower acoustic pressure input in figure 7(a–c), the concave region at the membrane poles expands outward slowly during the transient deformation of the cell, while at higher acoustic pressure inputs in figure 8(a–c), the concave region of the membrane poles expands outward rapidly during the transient deformation of the cell. In particular, comparing figures 7(b) and 8(b), the outward expansion of the RBC at 0.01 s under the high acoustic pressure input
$p_{0}=1.83\,\mathrm{MPa}$
is even higher than that at 0.1 s under the low acoustic pressure input
$p_{0}=0.58\,\mathrm{MPa}$
. As shown in figure 7(c), at lower acoustic pressure input, the stable deformation shape is still biconcave, with only a certain degree of outward expansion of the bipolar region. As shown in figure 8(c), at higher acoustic pressure input, the concave region on the membrane completely disappears when the steady state is reached, and the whole cell expands and stretches into a flattened disc shape.
Next, the second-order mean flow distribution pattern around the biconcave RBC during deformation is investigated. For this purpose, the streamlines and velocity contours of the second-order mean flow velocity
$\langle \boldsymbol{v}_{2}\rangle$
under different acoustic pressure inputs are shown in figures 7(d–f) and 8(d–f), where
$\langle \boldsymbol{v}_{2}\rangle$
includes both acoustic streaming and the Stokes flow driven by membrane deformation. As shown in figures 7(d) and 8(d), at the initial time
$t=0\,\mathrm{s}$
, the streamlines point outward at the pole, indicating that the Stokes flow driven by membrane deformation dominates the mean flow. As shown in figures 7(e) and 8(e), vortices are formed in the pole region and inward-pointing streamlines appear at the poles, because the acoustic streaming also affects the mean flow field in addition to the Stokes flow driven by membrane deformation. As shown in figures 7(f) and 8(f), when the RBC is gradually deformed to a steady state, the Stokes flow driven by membrane deformation disappears and the streamlines point inward, and then the mean flow
$\langle \boldsymbol{v}_{2}\rangle$
in the equilibrium state is only acoustic streaming. Figure 8(f) shows that under relatively large acoustic pressure input, the acoustic streaming structure around a flat disc-shaped RBC consists of two vortices inside and outside the cell.
Further, the first-order acoustic velocity distribution during the deformation of the biconcave RBC is investigated. As mentioned before, the time-averaged Piola–Kirchhoff stress is related to the first-order acoustic velocity
$\boldsymbol{v}_{1}$
, and therefore subsequent membrane stress analyses needs to be based on the first-order acoustic velocity. The first-order acoustic velocity
$\boldsymbol{v}_{1}$
is a time-dependent quantity that is not easy to visualize, whereas the root-mean-square (r.m.s.) acoustic velocity is used to visualize the first-order acoustic velocity. Figures 7(g–i) and 8(g–i) show the r.m.s. first-order acoustic velocity distribution during the deformation of the biconcave RBC under different acoustic pressure inputs. Comparing figures 7(g–i) and 8(g–i), it can be seen that there is a similarity in the distribution of the first-order acoustic velocities on the inner and outer sides of the cell under different acoustic pressure inputs. Near the equator, on the outer side of the cell, the red region indicates a larger velocity, and on the inner side of the cell, the blue region indicates a smaller velocity. This inner and outer velocity difference contributes to a larger time-averaged Piola–Kirchhoff stress difference. In addition, under different acoustic pressure inputs, there are also differences in the first-order acoustic velocity distributions on the inner and outer sides of the cell due to the different steady-state morphologies of the RBC after deformation. Taking the comparison between figures 7(i) and 8(i) as an example, the RBC in figure 7(i) is still biconcave after deformation, with lighter colours at the poles, which corresponds to larger first-order acoustic velocities, whereas the RBC in figure 8(i) is stretched into an ellipsoid shape, with relatively uniform blue colour inside the cell, which corresponds to relatively smaller and uniform first-order acoustic velocities. These differences in first-order velocities result in different time-averaged Piola–Kirchhoff stresses, which in turn correspond to different cell deformation morphologies, as will be further discussed in figure 9.
(a–b) Mean stress distribution (blue arrows) acting on the cell membrane for different input acoustic pressure amplitudes of
$p_{0}=0.58\,\mathrm{MPa}$
and
$p_{0}=1.83\,\mathrm{MPa}$
. The black line is the initial RBC shape and the red line is the deformed cell shape. (c–d) Normal stress (red solid line) and tangential stress (blue dotted line) of mean stress acting on the membrane versus the polar angle
$\theta =\arctan (r/z)$
at different input acoustic pressure amplitudes of
$p_{0}=0.58\,\mathrm{MPa}$
and
$p_{0}=1.83\,\mathrm{MPa}$
.
To further understand the cell deformation characteristics, figures 9(a) and 9(b) show the distributions of the time-averaged second-order mean stress on the cell membrane at the input acoustic pressure amplitude
$p_{0}=0.58\,\mathrm{MPa}$
and
$p_{0}=1.83\,\mathrm{MPa}$
, respectively. As shown in figure 9(a), when the input acoustic pressure amplitude is
$p_{0}=0.58\,\mathrm{MPa}$
, at both equators, the cell membrane is subjected to outward tensile mean stress, which is relatively large and symmetrically distributed. At both poles, the cell membrane is also subjected to outward tensile mean stress, which is relatively small and symmetrically distributed. As shown in figure 9(b), when the input acoustic pressure amplitude is
$p_{0}=1.83\,\mathrm{MPa}$
, the mean stress is mainly distributed at the equators, where the cell membrane is subjected to outward tensile mean stress, and the stress is large. At the poles, there is no outward mean stress acting on the cell membrane. In figures 9(a) and 9(b), the black line is the initial RBC shape and the red line is the deformed cell shape. Since the mean stress driving the cell deformation is mainly concentrated in the equatorial plane, the cell deformation mainly exhibits transverse stretching. When the input acoustic pressure amplitude is
$p_{0}=0.58\,\mathrm{MPa}$
, the cell shape remains biconcave and there is an outward mean stress difference at the pole, indicating that the mean stress does not exceed the elastic stress required to maintain the biconcave curvature of the membrane. When the input acoustic pressure amplitude is
$p_{0}=1.83\,\mathrm{MPa}$
, the mean stress increases by an order of magnitude, causing the cells to be stretched into an oblate disc shape. While there is no outward deformation driving force at the pole. Figures 9(a) and 9(b) show that different acoustic pressures inputs, resulting in different time-averaged second-order mean stresses distributed on the cell, can deform RBCs into distinct shapes.
Figures 9(c) and 9(d) further plot the normal and tangential stress components of the mean stress acting on the membrane as a function of the polar angle
$\theta =\arctan (r/z)$
. The positive normal component
$f_{n}=\boldsymbol{n}\boldsymbol{\cdot }\Delta \boldsymbol{f}$
points outward while the positive tangential component
$f_{t}=\boldsymbol{t}\boldsymbol{\cdot }\Delta \boldsymbol{f}$
points clockwise. As shown in figures 8(c) and 8(d), the normal stress is much larger than the tangential stress, indicating that the mean stress is dominated by the normal component. As mentioned above, the time-averaged Piola–Kirchhoff stress
$\langle \boldsymbol{P}_{2}\rangle$
driving the cell membrane deformation is closely related to the first-order acoustic velocity
$\boldsymbol{v}_{1}$
and the second-order mean flow velocity
$\langle \boldsymbol{v}_{2}\rangle$
. For the normal stress component of the Piola–Kirchhoff stress, in figures 9(a) and 9(b), the normal stress is maximum at the cell equator (
$\theta =\pm \pi /2$
), while in figures 7(i) and 8(i), the first-order acoustic velocity difference between the inside and outside of the cell membrane is also maximum at the cell equator (
$\theta =\pm \pi /2$
). Here, the normal stress in figure 9 is consistent with the first-order acoustic velocity distribution in figures 7 and 8, indicating that the first-order acoustic velocity plays a dominant role in the normal stress distribution. For the tangential stress component of the Piola–Kirchhoff stress, in figures 9(c) and 9(d), the maximum tangential stress also occurs near the equator of the RBC (
$\theta =\pm \pi /2$
), and correspondingly, the maximum second-order mean flow velocity occurs near the equator of the RBC (
$\theta =\pm \pi /2$
) in figures 7(f) and 8(f). The tangential stress distribution in figures 9(c) and 9(d) is consistent with the second-order mean flow velocity distribution in figures 7(f) and 8(f), indicating that the second-order mean flow velocity plays a dominant role in the tangential stress distribution. Actually, the tangential stress comes from the viscous shear stress generated by the second-order mean flow.
The effect of cytoplasm viscosity on the acoustic-induced cell deformation is investigated. Figure 10 shows the time evolution of the transverse and polar diameters of deformed RBCs at different cytoplasmic viscosities, both of which increase with increasing acoustic loading time and both of which eventually reach steady state deformation. As shown figures 10(a) and 10(b), as the cytoplasmic viscosity increases, the evolution of the transverse and polar diameters to their steady-state values slows down, i.e. the cytoplasmic viscosity affects the relaxation time. This is because during the transient deformation of the cell, significant acoustic streaming exists inside the cell (figures 7 and 8), accompanied by large internal velocity gradients and strong internal fluid shear forces. The internal cytoplasmic viscosity affects the viscous shear forces
$\mu (\langle \boldsymbol{\nabla }\boldsymbol{v}_{2}\rangle +\langle \boldsymbol{\nabla }\boldsymbol{v}_{2}^{T}\rangle )$
, thereby affecting the deformation relaxation time and cell deformation rate. As can be seen in figure 10(a), as the cytoplasmic viscosity increases, the steady-state value of the transverse diameter decreases. Due to the conservation of membrane enclosure volume, in figure 10(b), the steady-state value of the polar diameter increases with the increase of cytoplasmic viscosity to accommodate the change of the steady-state value of the transverse diameter. This is because during the steady-state deformation of the cell, internal acoustic streaming is significantly weakened, and internal viscosity can only affect cell deformation by generating small fluid shear forces through minor internal velocity gradients. Therefore, in the steady state, the internal cytoplasm viscosity has an effect on the steady-state deformation of biconcave RBCs, hindering transverse elongation and inducing polar expansion.
Time evolution of the acoustic-induced RBC deformation at different cytoplasm viscosities: (a) transverse diameter; (b) polar diameter. The input acoustic pressure amplitude is
$p_{0}=0.80\,\mathrm{MPa}$
.
Figure 11 shows the time evolution of the transverse and polar diameters of the biconcave RBC under different surrounding medium viscosities. Considering the actual liquid environment of RBCs (Bera et al. Reference Bera2022), a range of surrounding medium viscosity
$\mu _{{out}}$
from
$2\,\mathrm{mPa}\,\mathrm{s}$
to
$14\,\mathrm{mPa}\,\mathrm{s}$
is simulated. As shown in figures 11(a) and 11(b), the evolution of the transverse and polar diameters to their steady-state values slows down with increasing surrounding medium viscosity, i.e. the surrounding medium viscosity affects the relaxation time. This is because during the transient deformation of the cell, significant acoustic streaming exists outside the cell (figures 7 and 8), with large velocity gradients in the external acoustic streaming, leading to relatively strong external fluid shear forces. The viscosity of the surrounding medium influences these viscous shear forces, thereby affecting the deformation relaxation time and the cell deformation rate. As can be seen in figure 11(a), as the surrounding medium viscosity increases, the steady state value of transverse stretching decreases. As can be seen in figure 11(b), the larger the surrounding medium viscosity, the larger the steady-state value of the polar diameter. The effect of surrounding medium viscosity on cell deformation is similar to that of the cytoplasm. Due to the conservation of the membrane enclosure volume, the steady-state polar diameter increases with the increase of surrounding medium viscosity to accommodate the change in the steady-state value of the transverse diameter. This is because after the cell reaches steady-state deformation, the external acoustic streaming around the cell is significantly weakened, and the viscosity of the surrounding medium can only affect cell deformation by generating small fluid shear forces through minor external velocity gradients. Therefore, in the steady state, the surrounding medium viscosity has an effect on the steady-state deformation of biconcave RBCs, hindering transverse elongation and inducing polar expansion. Comparing figures 10 and 11, it can be seen that the surrounding medium viscosity has a greater effect on cell deformation than the cytoplasmic viscosity. From the mean flow field distributions in figures 7 and 8, it can be observed that at the end of deformation, i.e. when the cell deformation reaches a steady state, both internal and external acoustic streaming are significantly weakened, but the external acoustic streaming is stronger than the internal one, with corresponding external acoustic streaming velocity gradients also larger than the internal ones. Therefore, under the same conditions, the influence of external medium viscosity on cell deformation is greater than that of cytoplasmic viscosity.
Time evolution of the acoustic-induced RBC deformation at different buffer solution viscosities: (a) transverse diameter; (b) polar diameter. The input acoustic pressure amplitude is
$p_{0}=0.80\,\mathrm{MPa}$
.
Steady-state deformation of biconcave RBC with varying acoustic pressure amplitudes at different shear moduli
$G_{{s}}$
: (a) transverse diameter; (b) polar diameter. The bending stiffness is
$E_{{b}}=10k_{b}T$
.
The steady-state deformation of the biconcave RBC with varying acoustic pressure amplitudes at different shear moduli
$G_{{s}}$
is investigated, where the transverse diameter and polar diameter are shown in figures 12(a) and 12(b), respectively. Figure 12(a) shows that an increase in shear modulus leads to a decrease in the transverse diameter. Figure 12(b) shows that an increase in shear modulus also leads to a decrease in the polar diameter, which follows the same trend as the transverse diameter. As the acoustic input gradually increases, both the cell transverse diameter and polar diameter tend to steady state values, which is due to the limitation of cell surface area and cell volume. As shown in figure 12(b), for a given shear modulus (
$G_{\!s}=1.5, 3.5$
or
$6\,\unicode{x03BC} \mathrm{N}\,\mathrm{m}^{-1}$
), the polar diameter first increases with the increase of the acoustic input, while it suddenly jumps to a maximum value when the acoustic input continues to increase, and then decreases rapidly to converge to a steady state value. This complex deformation trend can be explained by the snap-through instability of RBCs (Das et al. Reference Das, Deshmukh and Thaokar2019). Similar to the snap-through instability phenomenon of RBCs in extensional flow, RBCs in an ultrasonic standing wave field also undergo snap-through instability as the acoustic input increases, i.e. the shapes of the two polar regions suddenly transition from biconcave to biconvex. The shear modulus can significantly affect the snap-through instability of RBCs. As shown in figure 12(b), as the shear modulus
$G_{\!s}$
increases, the critical acoustic pressure for instability increases and the onset of instability is delayed. This is because the higher the shear modulus, the greater the elastic energy barrier, and the greater the resistance that the acoustic-induced stress needs to overcome to trigger the shape transition from biconcave to biconvex.
To better understand the snap-through instability phenomenon of RBCs in an ultrasonic standing wave field, taking the shear moduli
$G_{{s}}=3.5\,\unicode{x03BC} \mathrm{N}\,\mathrm{m}^{-1}$
as an example, the steady equilibrium shape and mean stress distribution on the cell membrane at different acoustic pressure inputs are shown in figure 13. As shown in figure 13(a,b), when the acoustic pressure input is small, the RBC maintains a biconcave shape, and the variation of the polar diameter is dominated by the outward pointing stress at the cell poles. In figure 13(c), the snap-through instability of the RBC occurs at the critical value
$p_{0}^{2}=46\times 10^{10}\,\mathrm{Pa}^{2}$
, that is, the cell transforms into an oblate spheroidal shape, which is accompanied by a jump in the value of polar diameter. Fundamentally, the rapid enlargement of the concave areas, and even the outward protrusion of the poles, originates from the bistable mechanism of the phospholipid membrane in RBC membrane. Specifically, the initial oblate shape leads to a biconcave energy minimum shape, while the initial prolate shape leads to a peanut shape (Liu & Xin Reference Liu and Xin2022). In this study, the deformed RBC can be understood as a prolate shape. In figure 13(d–f), as the acoustic pressure input continues to increase, the polar diameter decreases due to the membrane-enclosed volume conservation.
Equilibrium shape and mean stress distribution on the cell membrane with the shear moduli
$G_{{s}}=3.5\,\unicode{x03BC} \mathrm{N}\,\mathrm{m}^{-1}$
at different acoustic inputs: (a)
$p_{0}^{2}=40\times 10^{10}\,\mathrm{Pa}^{2}$
, (b)
$p_{0}^{2}=45\times 10^{10}\,\mathrm{Pa}^{2}$
, (c)
$p_{0}^{2}=46\times 10^{10}\,\mathrm{Pa}^{2}$
, (d)
$p_{0}^{2}=50\times 10^{10}\,\mathrm{Pa}^{2}$
, (e)
$p_{0}^{2}=70\times 10^{10}\,\mathrm{Pa}^{2}$
, (f)
$p_{0}^{2}=80\times 10^{10}\,\mathrm{Pa}^{2}$
. The bending stiffness is
$E_{{b}}=10k_{b}T$
.
The effect of membrane bending stiffness is examined in figure 14, which depicts the steady-state deformation of the biconcave RBC with varying acoustic pressure amplitudes at different bending stiffness
$E_{{b}}$
, where the transverse diameter and polar diameter are shown in figures 14(a) and 14(b), respectively. Figure 14(a) shows that an increase in bending stiffness leads to a decrease in the transverse diameter of the steady state cell deformation, indicating that the bending stiffness greatly suppresses the transverse elongation. Figure 14(b) shows that, overall, an increase in bending stiffness leads to a decrease in the polar diameter of the steady-state cell deformation, reflecting a similar trend to the transverse diameter, suggesting that bending stiffness also suppresses the poleward elongation. When the acoustic input is gradually increased, both the cell transverse diameter and polar diameter tend to steady state values, which is due to the constraint of the cell surface area and cell volume.
Steady-state deformation of biconcave RBC with varying acoustic pressure amplitudes at different bending stiffness
$E_{{b}}$
: (a) transverse diameter; (b) polar diameter. The shear modulus is
$G_{{s}}=3.5\,\unicode{x03BC} \mathrm{N}\,\mathrm{m}^{-1}$
.
In fact, when the acoustic input is gradually increased, the trend of the pole diameter is a little different from that of the transverse diameter. As shown in figure 14(b), for a given bending stiffness (
$E_{{b}}=10k_{b}T$
or
$20k_{b}T$
), the polar diameter first increases with the increase of the acoustic input, and as the acoustic input continues to increase, the polar diameter suddenly jumps to a maximum value and then decreases rapidly to converge to a steady state value. As mentioned earlier, this phenomenon is due to the snap-through instability of RBCs. However, when the bending stiffness is large (
$E_{{b}}=30k_{b}T$
), the polar diameter increases monotonically with increasing acoustic input. Figure 15 shows steady equilibrium shape and mean stress distribution on the cell membrane with the bending stiffness
$E_{{b}}=30k_{b}T$
at different acoustic inputs. At the same acoustic input and shear stiffness, the cell morphologies of the RBCs in figure 15 are different from those in figure 13 due to the different bending stiffness. With increasing acoustic input, the biconcave RBC directly transforms into an oblate spheroid shape with no jump in polar diameter, and the outward pointing stress at the cell poles disappears as the membrane surface flattens. This can be explained by the fact that the snap-through instability of RBCs is essentially the result of a competition between the bending force
$\sum _{s}^{b}$
that maintains their shape and the acoustic-induced deformation force
$\Delta \boldsymbol{f}$
. The increase in membrane bending stiffness enhances the energy barrier that prevents shape evolution beyond the shape shown in figure 13(c). It can be concluded that the increase in bending stiffness can prevent snap-through instability of RBC due to its resistance to deformation.
Equilibrium shape and mean stress distribution on the cell membrane with the bending stiffness
$E_{{b}}=30k_{b}T$
at different acoustic inputs: (a)
$p_{0}^{2}=40\times 10^{10}\,\mathrm{Pa}^{2}$
, (b)
$p_{0}^{2}=45\times 10^{10}\,\mathrm{Pa}^{2}$
, (c)
$p_{0}^{2}=50\times 10^{10}\,\mathrm{Pa}^{2}$
, (d)
$p_{0}^{2}=60\times 10^{10}\,\mathrm{Pa}^{2}$
, (e)
$p_{0}^{2}=70\times 10^{10}\,\mathrm{Pa}^{2}$
, (f)
$p_{0}^{2}=80\times 10^{10}\,\mathrm{Pa}^{2}$
. The shear modulus is
$G_{{s}}=3.5\,\unicode{x03BC} \mathrm{N}\,\mathrm{m}^{-1}$
.
5. Concluding remarks
In this paper, the mean deformation dynamics of biconcave RBCs in a viscous fluid driven by an ultrasonic standing wave is studied by theoretical and numerical models. Considering the true physiological shape of RBCs with a biconcave geometric structure, this paper overcomes the challenge of modelling the nonlinear acoustic–solid coupling of complex biconcave curved shells and establishes the first computational model of the fully dynamic deformation of biconcave RBCs in a single standing wave acoustic field. In the ultrasonic standing wave field, the particles of the cells and fluids not only oscillate on a fast time scale, but also exhibits time-averaged response on a slow time scale due to an acoustic nonlinear effect. The acoustic perturbation method based on the generalized Lagrangian formulation is used to decompose the particle motion into first-order ultrasonic propagation and second-order time-averaged mean response. The equations governing the time-averaged mean dynamics of the fluid–cell system are derived from the Navier–Stokes equations, which are then reformulated into weak forms and solved using the finite element method. The developed computational model is able to handle well the coupling of three distinct physical fields, including the acoustic propagation field, cell deformation field and mean fluid flow field.
The study on acoustic deformation of biconcave RBCs reveals that when the acoustic pressure input is relatively small, the biconcave RBCs are stretched but still maintain the biconcave shape. When the acoustic pressure input is large enough, the biconcave RBCs are stretched into a flat disc shape. In addition, due to the bistable mechanism of biconcave RBCs, there exists a critical input acoustic pressure. When the acoustic pressure amplitude reaches the critical threshold, the RBCs undergo a significant shape transition, where the biconcave cavities at the poles open up and bulge outwards. This computational result marks the first prediction and discovery of the snap-through instability phenomenon in biconcave RBCs driven by ultrasound. By analysing the stress distribution on the biconcave RBC membrane, it is found that the stress driving the deformation of the cell membrane is mainly the normal component of the time-averaged Piola–Kirchhoff stress, which leads to the transverse elongation and polar expansion of the RBCs. Furthermore, the normal stress is mainly related to the first-order acoustic velocity, and the tangential stress is mainly related to the second-order mean flow velocity. Detailed analyses of the deformed RBC shape show that cytoplasmic viscosity and surrounding medium viscosity can both increase the viscous drag, thereby hindering transverse elongation in relaxation time and steady-state value. In addition to the properties of the cytoplasm and the surrounding medium, parametric studies are performed to investigate the mechanical properties of the membranes. The results show that both the shear modulus and the bending stiffness hinder the acoustic deformation of the cell by inhibiting transverse elongation and polar expansion to maintain its biconcave shape. The increase in bending stiffness can even prevent snap-through instability of RBC. This work develops a theoretical computational method for ultrasound-driven mechanical deformation of complex morphological cells in viscous fluids.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2025.10717.
Acknowledgements
This work was supported by the National Natural Science Foundation of China (52375122), the Natural Science Basic Research Program of Shaanxi Province (2023-JC-JQ-06 and 2025JC-QYXQ-001) and the Fundamental Research Funds for the Central Universities (LX6J013).
Declaration of interest
The author reports no conflict of interest.
Appendix. Spatial mesh density and time step convergence verification
To ensure calculation accuracy and efficiency, convergence verification of both spatial mesh density and time step is performed. The initial biconcave RBC is positioned at the acoustic pressure node. Prior to applying the ultrasonic standing wave field, the fluid is stationary with zero particle velocity (
$\boldsymbol{v}$
= 0). The input acoustic pressure amplitude is set to
$p_{0}=0.8\,\mathrm{MPa}$
.
First, the spatial mesh convergence verification of the model is performed. To ensure that the computational results accurately capture cell deformation phenomena, the mesh resolution around the cell membrane must be sufficiently fine, especially in the acoustic viscous boundary layer. The effects of the cell membrane line element mesh density N and fluid element mesh density M inside and outside the cell membrane on the transient and steady-state deformation of RBCs are shown in figure 16, where three distinct mesh configurations are tested. In figure 16, the mesh number N refers to the number of cell membrane line elements, while the mesh number M refers to the number of moving mesh elements in the fluid domain
$\mathcal{B}_{2}+\mathcal{B}_{3}$
inside and outside the cell in figure 3(a). In the process of mesh generation, we artificially set the size and density of the cell membrane line elements and the internal and external fluid mesh elements to be strictly synchronized and strongly correlated, i.e. the maximum rate of change in the size of adjacent mesh elements is controlled to be unified as 1.1, thereby achieving the synchronous growth of the mesh density of the fluid inside and outside the membrane with the line mesh density of membrane. Specifically, as shown in figure 17, when the number of membrane mesh line elements is
$N=100 , N=140$
and
$N=180$
, the corresponding number of fluid mesh elements is
$M=12\,157 , M=17\,478$
and
$M=22\,945$
, respectively. Therefore, the convergence verification of the cell membrane mesh for cell deformation is conducted simultaneously with the convergence verification of the fluid mesh densities inside and outside the cell membrane. Under different cell membrane line element densities and fluid element densities, figures 16(a) and 16(b) show the evolution of transverse diameter and polar diameter over time, while figure 16(c) shows the steady-state deformation profile of the RBC. The time evolution of the transverse diameter, polar diameter and steady-state deformation profile of the cell almost overlap at the three different mesh densities. That is, the transient evolution curves and steady-state shape profile at (
$N=100 , M=12\,157$
) are the same as those at (
$N=140 , M=17\,478$
) and (
$N=180 , M=22\,945$
), indicating that both the cell membrane mesh line density
$N=100$
and the fluid mesh density inside and outside the cell membrane
$M=12\,157$
meet the computational convergence requirements.
Effect of membrane line element count (N) and fluid element count (M) in the domain
$\mathcal{B}_{2}+\mathcal{B}_{3}$
(figure 3
a) on RBC deformation: (a) time evolution of transverse diameter; (b) time evolution of polar diameter; (c) steady-state deformation profile. The time step is
${\unicode[Arial]{x0394}} t=5\times 10^{-4}\,\mathrm{s}$
, and the input acoustic pressure amplitude is
$0.8\,\mathrm{MPa}$
.
Schematic diagram of different mesh densities in spatial mesh density convergence verification. This diagram shows an enlarged view of the mesh on both sides of the biconcave RBC membrane, corresponding to different mesh densities in figure 3(c).
Next, the time step convergence of the model is verified. As the choice of time step also affects the accuracy and efficiency of the calculation, it is necessary to verify the time step convergence of the model before performing numerical calculations. Figure 18 shows the effect of time step on the transient and steady-state deformation of RBCs, where the time steps
${\unicode[Arial]{x0394}} t=5\times 10^{-4}, 1\times 10^{-4}$
and
$0.5\times 10^{-4}\,\mathrm{s}$
are tested. At different time steps, figures 18(a) and 18(b) show the time evolution of the transverse diameter and polar diameter, while figure 18(c) shows the steady-state deformation profile of the cell. It can also be seen that the transient evolution curves and steady-state shape profile at
${\unicode[Arial]{x0394}} t=5\times 10^{-4}\,\mathrm{s}$
are the same as those at
${\unicode[Arial]{x0394}} t=1\times 10^{-4}\,\mathrm{s}$
and
${\unicode[Arial]{x0394}} t=0.5\times 10^{-4}\,\mathrm{s}$
, indicating that the time step
${\unicode[Arial]{x0394}} t=5\times 10^{-4}\,\mathrm{s}$
meets the computational convergence and accuracy requirements. Therefore, under the premise of satisfying computational convergence and accuracy, to pursue shorter computation time, the number of membrane line elements is selected as
$N=100$
, the number of fluid elements inside and outside the cell membrane is selected as
$M=12\,157$
and the time step is selected as
$t=5\times 10^{-4}\,\mathrm{s}$
in the calculations.
Effect of time step on RBC deformation: (a) time evolution of transverse diameter; (b) time evolution of polar diameter; (c) steady-state deformation profile. The membrane line element count is
$N=100$
, the fluid element count around cell membrane is
$M=12\,157$
, and the input acoustic pressure amplitude is
$0.8\,\mathrm{MPa}$
.
