Study of the mechanism of action of sand therapy on atherosclerosis based on the two-phase flow-Casson model

Abstract

BACKGROUND:

Sand therapy is a non-pharmacological physiotherapy method that uses the natural environment and resources of Xinjiang to treat through the heat transfer and magnetic effects of sand.

OBJECTIVE:

Employing the two-phase flow-Casson blood flow model, we investigate the mechanism of atherosclerosis prevention via sand therapy, offering a biomechanical theoretical rationale for the prevention of atherosclerosis through sand therapy via the prism of computational fluid dynamics (CFD).

METHODS:

Sand therapy experiments were conducted to obtain popliteal artery blood flow velocity, and blood was considered as a two-phase flow composed of plasma and red blood cells, and CFD method was applied to analyze the hemodynamic effects of Casson’s blood viscosity model before and after sand therapy.

RESULTS:

(1) The blood flow velocity increased by 0.24 m/s and 0.04 m/s at peak systolic and diastolic phases, respectively, after sand therapy; the axial velocity of blood vessels increased by 28.56% after sand therapy. (2) The average red blood cell viscosity decreased by 0.00014 Pa ⋅ s after sand therapy. (3) The low wall shear stress increased by 1.09 Pa and the high wall shear stress reached 41.47 Pa after sand therapy. (4) The time-averaged wall shear stress, shear oscillation index and relative retention time were reduced after sand therapy.

CONCLUSION:

The increase of blood flow velocity after sand therapy can reduce the excessive deposition of cholesterol and other substances, the decrease of erythrocyte viscosity is beneficial to the migration of erythrocytes to the vascular center, the increase of low wall shear stress has a positive effect on the prevention of atherosclerosis, and the decrease of time-averaged wall shear stress, shear oscillation index and relative retention time can reduce the occurrence of thrombosis.

Keywords

Sand therapy two-phase flow non-Newtonian fluid atherosclerosis erythrocytes

1. Introduction

Sand burial therapy (sand therapy for short) is a non-pharmacological physiotherapy method that uses the natural environment and resources of Xinjiang to treat through the heat transfer and magnetic effects of sand. Sand therapy can treat rheumatism and has analgesic, sedative, anti-inflammatory, decongestive, antihypertensive, hypolipidemic and antidiarrheal effects [1]. In the study of sand therapy on haemodynamics, Liu et al. [2] studied the haemodynamic effects of sand therapy on human femoral artery, and Qiao et al. [3] found that sand therapy has reduced the tendency of atherosclerotic plaque enlargement. Wu et al. [4] studied and statistically analysed the haemodynamic parameters of popliteal artery before and after sand therapy in 58 subjects, and found that sand therapy could increase the blood circulation rate and improve the blood circulation. However, current studies on the haemodynamic aspects of sand therapy consider blood to be a Newtonian fluid with constant viscosity; in fact, blood is a non-Newtonian fluid consisting of e.g. plasma, erythrocytes and platelets, and a single Newtonian fluid model may overestimate the wall shear stress.

Atherosclerosis is a lipid-driven, chronic inflammatory disease that forms plaques, usually in medium and large arteries, and is a major cause of ischaemic heart disease and stroke. An important and potential factor in the formation of atherosclerosis is wall shear stress, where low wall surface weakens the mass transfer between blood and the vessel wall, high wall shear stress may damage the vessel wall and endothelial cells ultimately resulting in thrombosis, and there is a strong association between local intimal thickening and wall shear stress. In addition, studies [5,6] have shown that atherosclerotic disease is a long-term developmental process, and that timely intervention at an early stage can be effective in controlling the disease or even improving the condition.

Although some achievements have been made in the study of the effect of sand therapy on haemodynamics, basically the blood is regarded as a single-phase Newtonian fluid and there is no specific analysis of the mechanism of sand therapy on atherosclerosis. In this paper, we use Phillips CX50 colour Doppler ultrasound to collect the popliteal artery inlet velocities before and after sand therapy, and the blood is regarded as a two-phase non-Newtonian model, which will be more able to reveal the mechanism of the effect of sand therapy on atherosclerosis. Majid Abbasian et al. [7] investigated the haemodynamics of Newtonian fluids and 15 non-Newtonian fluids, and found that the most accurate results were obtained by using the Casson blood viscosity model in computational fluid dynamics (CFD) calculations by comparing it with clinical experimental data. Therefore, this paper extends blood into a two-phase flow model composed of plasma (liquid phase) and red blood cells (solid phase), and other blood components such as platelets are ignored as their volume fraction is less than 1%, and analyses the changes in hemodynamic parameters of Casson blood viscosity model before and after sand therapy to explore the physiotherapy mechanism of sand therapy on atherosclerotic diseases, and to provide the theoretical basis for biomechanics in the prevention of atherosclerosis.

2. Materials and methods

2.1. Model reconstruction

Fifty-eight (age: 33 ± 14) healthy volunteers without popliteal artery stenosis and other known cardiovascular diseases were selected to participate in this study. Blood flow parameters such as popliteal artery blood flow velocity, intravascular diameter and other blood flow parameters before and after sand therapy were collected and statistically analyzed using Doppler ultrasonography with the subjects lying on their backs to calculate the Reynolds number of blood flow before and after sand therapy and to determine the state of blood flow. The DICOM data of the popliteal artery of a healthy volunteer without popliteal stenosis and other known cardiovascular diseases were obtained using CT scans, and the DICOM data were imported into Mimics 17.0, and the 3D geometric configuration of the vessel was reconstructed by selecting the popliteal artery (POA), anterior tibial artery (ATA), and posterior tibial artery (PTA) sites of the lower extremity. A solid model of the vessel wall with 1 mm thickness was constructed using Geomagic Studio, and the geometric model of the bifurcation vessel of the popliteal artery in the subject was obtained as shown in Fig. 1.

Fig. 1.

Establishment of geometric model. A. Construction of rouge artery model based on CT data. B. Rouge artery model after construction and smoothing.

2.2. Governing equations

In this paper, the multiphase model is extended to a two-phase flow consisting of plasma and red blood cells (RBC), and the multiphase CFD method proposed by Gidaspow [8] and Anderson and Jackson [9] is applied to describe the multiphase hemodynamics of blood. The multiphase model includes the volume fraction of each phase, and the mechanisms of momentum, mass and energy exchange between the phases.

The equation of continuity between the phases (k = plasma, RBC) is given by the following equation: $\begin{eqnarray}\frac{\partial ({\rho}_{k}{\varepsilon}_{k})}{\partial t}+{\nabla}\cdot ({\rho}_{k}{\varepsilon}_{k}\overset{\rightarrow }{v}_{k})=0\end{eqnarray}$ (1) where 𝜌 is the density, 𝜀 is the volume fraction, and t is time. Also the volume fractions of each phase must add up to 1. $\begin{eqnarray}\mathop{\mathop{\sum }^{\mathit{np}}}_{k=1}{\varepsilon}_{k}=1\end{eqnarray}$ (2) Where np is 2 and the volume fraction of one phase cannot be occupied by other phases.

The momentum equation for each phase (k = plasma, RBC) is given by the following equation [10]: $\begin{eqnarray}\frac{\partial ({\rho}_{k}{\varepsilon}_{k}\overset{\rightarrow }{v}_{k})}{\partial t}+{\nabla}\cdot ({\rho}_{k}{\varepsilon}_{k}\overset{\rightarrow }{v}_{k}\overset{\rightarrow }{v}_{k})=-{\varepsilon}_{k}{\nabla}p+{\nabla}\cdot {\varepsilon}_{k}\overline{\overline{{\tau}}}_{k}+{\varepsilon}_{k}{\rho}_{k}\overset{\rightarrow }{g}+\mathop{\sum }_{l\neq k}{\beta}_{kl}(\overset{\rightarrow }{v}_{l}-\overset{\rightarrow }{v}_{k})+\overset{\rightarrow }{F}_{k}\end{eqnarray}$ (3) Where p is the pressure, $\bar{\bar{{\tau}}}_{k}$ is the stress tensor, $\vec{g}$ is the gravitational force, and $\vec{F}_{k}$ is the extemal force (lift, virtual mass force), and 𝛽_kl is the interphase momentum exchange coefficient.

The mixed density of blood is derived from the sum of the volume fraction weighted densities of plasma and RBCs [11]: $\begin{eqnarray}{\beta}_{mix}=\mathop{\mathop{\sum }^{\mathit{np}}}_{k=1}{\varepsilon}_{k}{\beta}_{k}.\end{eqnarray}$ (4) In this study, a blood mixture density of 1080 kg∕m³ [12 ] and a plasma density of 1000 kg∕m³, with an inlet RBC blood volume ratio of 𝜑(RBC) = 0.45 was used [13].

The stress tensor for each stage is given in Newtonian form as follows [7]: $\begin{eqnarray}\overline{\overline{{\tau}}}=2{\mu}_{k}\overline{\overline{D}}_{k}-\frac{2}{3}{\mu}_{k}tr(\overline{\overline{D}}_{k}).\end{eqnarray}$ (5) The deformation tensor $\bar{\bar{D}}_{k}$ is given by the following equation [14]: $\begin{eqnarray}\overline{\overline{D}}_{k}=\frac{1}{2}[{\nabla}\overset{\rightarrow }{v}_{k}+({\nabla}\overset{\rightarrow }{v}_{k})^{T}].\end{eqnarray}$ (6) Viscosity is an important parameter of blood and the Carreau–Yasuda viscosity model is used to describe the behavior of non-Newtonian fluids. Thus, the blood mixture viscosity 𝜂 is [15]: $\begin{eqnarray}{\eta}=\frac{{\varepsilon}_{\mathit{RBC}}{\mu}_{\mathit{RBC}}+{\varepsilon}_{\mathit{Plasma}}{\mu}_{\mathit{Plasma}}}{{\mu}_{\mathit{plasma}}}=m[1+({\lambda}\dot{{\rm\gamma}})^{2}]^{\frac{n-1}{2}}.\end{eqnarray}$ (7)

Plasma is considered as a Newtonian fluid with constant viscosity 𝜇_plasma = 0.0035 Pa ⋅ s [16], $\dot{{\gamma}}$ is the shear rate (s ⁻¹), 𝜆 is a constant time parameter of 0.110 s [12], and m and n are parameters related to the volume of red blood cells. $\begin{eqnarray}\displaystyle n & = & \displaystyle 0.8092{\varepsilon}_{\mathit{RBC}}^{3}-0.8246{\varepsilon}_{\mathit{RBC}}^{2}-0.3503{\varepsilon}_{\mathit{RBC}}+1\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle m & = & \displaystyle 122.28{\varepsilon}_{\mathit{RBC}}^{3}-51.213{\varepsilon}_{\mathit{RBC}}^{2}+16.305{\varepsilon}_{\mathit{RBC}}+1.\end{eqnarray}$ (9)

Particles moving in a viscous fluid are subjected to drag forces and the Schiller–Naumann model was used to calculate the trapping forces between the continuous phase (cp) and the discrete phase (dp) [17]. $\begin{eqnarray}{\beta}=\frac{3}{4}C_{d}\frac{{\rho}_{\mathit{cp}}{\varepsilon}_{\mathit{cp}}{\varepsilon}_{dp}|\overset{\rightarrow }{v}_{\mathit{cp}}-\overset{\rightarrow }{v}_{dp}|}{d_{dp}{\phi}}\end{eqnarray}$ (10) where C _d is related to the Reynolds number Re: $\begin{eqnarray}\displaystyle C_{d} & = & \displaystyle \frac{24}{\mathit{Re}}[1+0.15\mathit{Re}^{0.687}]\quad \mathit{Re}<1000\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle C_{d} & = & \displaystyle 0.44\quad \mathit{Re}\geq 1000\end{eqnarray}$ (12) $\begin{eqnarray}\displaystyle \mathit{Re} & = & \displaystyle \frac{{\rho}_{\mathit{cp}}d_{dp}|\overset{\rightarrow }{v}_{\mathit{cp}}-\overset{\rightarrow }{v}_{dp}|{\phi}}{{\mu}_{\mathit{cp}}}\end{eqnarray}$ (13) where d is the diameter of the blood cell and is the dynamic shape factor. The erythrocyte shape factors used to calculate the resistance are assumed to be spherical, i.e., 𝛷 = 1 in the absence of damage.

The other forces acting on each phase include the lift force [18] and the virtual mass force [19], and the lift force expression is given by: $\begin{eqnarray}\displaystyle \displaystyle \overset{\rightarrow }{F}_{\mathit{lift},k} & = & \displaystyle -0.5{\varepsilon}_{m}{\rho}_{k}(\overset{\rightarrow }{v}_{k}-\overset{\rightarrow }{v}_{m})\times \left.\left({\nabla}\times \overset{\rightarrow }{v}_{k}\right)\right|_{k\neq m}\nonumber\\ \displaystyle \displaystyle \overset{\rightarrow }{F}_{\mathit{lift},k} & = & \displaystyle -\overset{\rightarrow }{F}_{\mathit{lift},m}.\end{eqnarray}$ (14)

The virtual mass force expression is given by: $\begin{eqnarray}\displaystyle \overset{\rightarrow }{F}_{vm,k} & = & \displaystyle 0.5{\varepsilon}_{m}{\rho}_{k}\left.\left(\frac{d_{k}\overline{v}_{k}}{\text{d}t}-\frac{d_{m}\overset{\rightarrow }{v}_{m}}{\text{d}t}\right)\right|_{k\neq m}\nonumber\\ \displaystyle \overset{\rightarrow }{F}_{vm,k} & = & \displaystyle -\overset{\rightarrow }{F}_{vm,m}\end{eqnarray}$ (15) $\begin{eqnarray}\displaystyle \frac{\text{d}\overset{\rightarrow }{v}}{\text{d}t} & = & \displaystyle \frac{\partial \overset{\rightarrow }{v}}{\partial t}+(\overset{\rightarrow }{v}\cdot {\nabla})\overset{\rightarrow }{v}.\end{eqnarray}$ (16)

Casson blood viscosity model [20]: $\begin{eqnarray}\sqrt{{\mu}}=\sqrt{\frac{{\tau}_{0}}{\dot{{\gamma}}}}+\sqrt{{\mu}_{0}}\end{eqnarray}$ (17) where 𝜇 is the kinetic viscosity, $\dot{{\gamma}}$ is the shear strain rate, is the yield stress 𝜏₀ = 0.004 Pa, and 𝜇₀ is the initial blood viscosity 𝜇₀ = 0.0035 Pa ⋅ s.

2.3. Fluid-structure coupling

The elastic vessel wall will deform periodically under the action of pulsating blood flow, while the pressure and stress on the vessel wall surface will change periodically, and in turn the deformation of the vessel wall will affect the change of the blood flow field, i.e., the flow-structure coupling between the pulsating blood flow and the elastic vessel wall.

The arbitrary Lagrangian–Euler method (ALE) is the most commonly used method for solving the fluid–solid coupling problem [21,22]. For the solid region, the control equations satisfy the law of conservation of momentum. $\begin{eqnarray}{\rho}_{s}\frac{\partial ^{2}\vec{X}_{s}}{\partial t^{2}}={\nabla}\cdot \overline{\overline{{\sigma}_{s}}}+\vec{f}_{s}^{B}\end{eqnarray}$ (18) where the s subscript denotes the solid portion of the vessel wall, 𝜌_s is the density, $\vec{X}_{s}$ is the displacement vector, $\bar{\bar{{\sigma}}}_{s}$ is the Cauchy stress tensor, and $\vec{f}_{s}^{B}$ is the volume force exerted on the solid.

For the fluid region, the continuity equation and the momentum conservation equation are still satisfied. Considering that this paper applies a two-fluid model to describe the blood flow, Eq. (3) can be rewritten in an arbitrary Lagrangian–Eulerian form as follows: $\begin{eqnarray}\frac{\partial ({\rho}_{k}{\varepsilon}_{k}\vec{v}_{k})}{\partial t}+{\rho}_{k}{\varepsilon}_{k}(\vec{v}_{k}-\vec{v}_{g}{\nabla}\cdot \vec{v}_{k})=-{\varepsilon}_{k}{\nabla}p_{k}+{\nabla}\cdot \overline{\overline{{\tau}}}_{k}+{\varepsilon}_{k}{\rho}_{k}\vec{g}+\sum _{l\neq k}^{}{\beta}_{kl}(\vec{v}_{l}-\vec{v}_{k})+\vec{F}_{k}+\vec{f}_{f}^{B}\end{eqnarray}$ (19) where $\vec{v}_{g}$ is the grid movement velocity in the ALE method and $\vec{f}_{f}^{B}$ denotes the volume force on the fluid.

In this paper, the finite element software ANSYS Transient Structural and CFX modules are used to calculate the solid and fluid domains, respectively, and the data are exchanged between the two under the ANSYS workbench platform through a predefined intersection interface, i.e., the inner wall surface of the vessel, thus realizing the fluid–solid coupling calculation. The solid and fluid regions satisfy the boundary conditions of no slip and no infiltration at the intersection: $\begin{eqnarray} \displaystyle \overset{\rightharpoonup}{d_{f}} & = & \displaystyle \overset{\rightharpoonup}{d_{s}}(x,t)\in S_{i} \end{eqnarray}$ (20) $\begin{eqnarray} \displaystyle \overset{\rightharpoonup}{{\sigma}_{f}} & = & \displaystyle \overset{\rightharpoonup}{{\sigma}_{s}}(x,t)\in S_{i} \end{eqnarray}$ (21) where the subscripts f, s denotes fluid domain and solid domain respectively, $\stackrel{\rightharpoonup}{d_{f}},\stackrel{\rightharpoonup}{d_{s}}$ denotes displacement, $\stackrel{\rightharpoonup}{{\sigma}_{f}},\stackrel{\rightharpoonup}{{\sigma}_{s}}$ denotes normal traction, and S _i denotes fluid–solid coupling intersection.

2.4. Numerical methods

The reconstructed popliteal artery model was loaded in the CFD software CFX 17.0, and the popliteal artery inlet velocity waveforms before and after sand therapy were fitted according to the experimental data, as shown in Fig. 2. The vessel wall is an isotropic linear elastic material with a wall density of 1150, an elastic modulus of 5 MPa, and a Poisson’s ratio of 0.45 [23]. The reference pressure value in the computational domain is 13339 Pa [24], and the relative pressure at the outlet is 0 Pa [25], the contact surface between blood and the vessel wall is a non-permeable, non-slip boundary, and the collision between red blood cells and the vessel wall is a completely elastic collision.

Fig. 2.

Popliteal artery inlet velocity before and after sand therapy.

Fig. 3.

Mesh model.

When generating the computational mesh, a tetrahedral mesh was used to mesh the model with a cell size of 0.4 mm, and the mesh was encrypted near the contact between the blood and the vessel wall with three layers of encryption, and the mesh effect is shown in Fig. 3. In order to eliminate the influence of the initialization error on the calculation results, three consecutive cardiac cycles were calculated, and the calculation results of the third cardiac cycle were taken for analysis.

Fig. 4.

Blood flow velocity distribution. A: Schematic diagram of popliteal artery cross-section; B: mean velocity distribution of blood flow before and after sand therapy, subscript is cardiac cycle before sand therapy, superscript is cardiac cycle after sand therapy; C: axial velocity distribution of blood flow at peak systole.

3. Results

3.1. Velocity of flow

Due to the increase in heart rate after sand therapy, the cardiac cycle decreased from 0.789 s to 0.6516 s [3], and to compare the changes in blood flow velocity after sand therapy, the mean velocity distribution of popliteal artery during the cardiac cycle before and after sand therapy was plotted as shown in Fig. 4B. The blood flow velocity increased after sand therapy, and at the peak of systole, the blood flow velocity was 0.78 m/s before sand therapy and 1.02 m/s after sand therapy, with an increase in blood flow velocity of 0.24 m/s. At the peak of diastole, the blood flow velocity was 0.25 m/s before sand therapy and 0.29 m/s after sand therapy, with an increase in blood flow velocity of 0.04 m/s. To facilitate comparison of the axial velocity distribution, the vessel diameter was normalized and the axial velocity distribution of the popliteal artery cross-section at peak systole was plotted as shown in Fig. 4C. The axial velocity of blood increased by 28.56% after sand therapy.

Figure 5 shows the flow diagrams of blood at characteristic moments before and after sand therapy, with flat and stable blood flow before and after sand therapy at peak systole and peak diastole. Figure 5B shows that unstable blood flow occurs at the vascular bifurcation and vascular outlet at the systolic trough, and Fig. 5D shows that blood flow velocity decreases with the prolongation of diastolic phase, and vortex flow occurs in the vascular bifurcation vessels before sand therapy and decreases in the bifurcation vessels after sand therapy.

Fig. 5.

Comparison of flow lines before and after sand therapy. A: Blood flow line during peak systole (T = 1.80 s before sand therapy; T = 1.58 s after sand therapy); B: blood flow line during systolic trough (T = 2.07 s before sand therapy; T = 1.82 s after sand therapy); C: blood flow line during peak diastole (T = 2.16 s before sand therapy; T = 1.86 s after sand therapy); D: blood flow line during mid-diastole (T = 2.30 s before sand therapy; T = 2.0 s after sand therapy).

3.2. Erythrocyte viscosity

Figure 6 shows the variation of erythrocyte viscosity with radial position at the popliteal artery cross-section, which shows that the erythrocyte viscosity near the vessel wall is relatively low, and the erythrocyte viscosity in the central part of the vessel is high, reaching a maximum viscosity of 0.0052 Pa ⋅ s at a dimensionless radial position of about 0.5. The average erythrocyte viscosity decreased by 0.00014 Pa ⋅ s after sand therapy in the interval of −0.5 to 0.5 dimensionless radial position.

Fig. 6.

RBC viscosity of popliteal artery cross section.

3.3. Wall shear stress

In order to study the distribution of wall shear stress (WSS) in the vessel wall, 56 data extraction points were set along the vessel as shown in Fig. 7, with data points 1∼19 for the lateral anterior tibial artery, data points 20∼36 for the lateral posterior tibial artery, and data points 37∼56 for the medial branch vessel. The WSS at the peak of systole before and after sand therapy were extracted separately and plotted as shown in Fig. 7. The trend of WSS distribution before and after sand therapy was more or less the same, and the WSS increased after sand therapy, with a low WSS at the bending part of the vessel (Point6∼Point14 and Point26∼Point28) with a minimum value of 0.13 Pa (Point 27) before sand therapy and a WSS of 1.22 Pa after sand therapy. The maximum value was 30.39 Pa (Point 47) before sand therapy and reached 41.47 Pa after sand therapy.

Fig. 7.

Distribution of data points.

3.4. Time-averaged wall shear stress, oscillatory shear index and relative retention time

WSS is an important parameter of hemodynamics, and researchers have further derived WSS-based parameters to describe different characteristics of WSS, such as Time-averaged Wall Shear Stress (TAWSS), Oscillatory Shear Index (OSI), and Relative Residence Time (RRT).

TAWSS refers to the average value of WSS over the entire cardiac cycle, while WSS is the tangential frictional force of blood flow on endothelial cells, which is decisive for the histophysiological changes of endothelial cells in the lesioned area and is generally considered to be a key factor contributing to atheromatous plaques. The distribution of TAWSS is more clinically significant than the transient wall shear. The mathematical expression is defined as: $\begin{eqnarray}\mathit{TAWSS}=\frac{1}{T}\int _{0}^{T}|{\tau}_{w}|\text{d}t.\end{eqnarray}$ (22)

OSI is the force at transient shear angles greater than 90° during the cardiac cycle and ranges from 0 to 0.5. Higher OSI values are thought to be determinants of the cyclic stresses that lead to functional deficits and altered cell structure in endocardial cells. The mathematical expression is defined as: $\begin{eqnarray}\mathit{OSI}=0.5\times \left(1-\frac{\left|{\displaystyle \frac{1}{T}}\displaystyle \int _{0}^{T}{\tau}_{w}\text{d}t\right|}{{\displaystyle \frac{1}{T}}\displaystyle \int _{0}^{T}|{\tau}_{w}|\text{d}t}\right).\end{eqnarray}$ (23)

The calculation of RRT combines the WSS parameters and the OSI parameters to describe the relative residence time of the blood flow. Relative residence time is closely related to thrombus formation, with higher values of this value indicating a longer residence time and a higher probability of thrombus formation. The mathematical expression is defined as: $\begin{eqnarray}RRT=\frac{1}{(1-2\times \mathit{OSI}\times \mathit{TAWSS})}\end{eqnarray}$ (24) where T is the pulsation period and 𝜏_w is the shear force vector.

4. Discussion

As a non-pharmacological therapy, sand therapy has been used for a long time to improve blood circulation, and for the study of sand therapy hemodynamics, Brgeby [26] analyzed the effect of sand therapy on femoral artery hemodynamics, which provided a theoretical basis for the subsequent sand therapy studies. Qiao et al. [3] analyzed the effect of sand therapy on hemodynamic parameters of different stenosis models of femoral bifurcation vessels and found that the trend of atherosclerotic plaque enlargement was attenuated after sand therapy. Based on this paper, the two-phase flow of blood and the non-Newtonian characteristics of Casson model were considered, and the physiotherapeutic mechanism of sand therapy on atherosclerotic diseases was further investigated to provide a theoretical basis for the prevention of atherosclerosis.

Figure 4B shows that the blood flow velocity increased after sand therapy, and the increase of blood flow velocity during systole was the largest, and the increase of blood flow velocity can promote blood circulation thus treating cardiovascular diseases caused by poor blood circulation, while the increase of blood flow velocity. Figure 4C shows that the blood flow velocity increases in the center of the blood vessel after sand therapy, and according to Bernoulli’s equation, the increase in blood flow velocity causes a decrease in pressure in the center of the blood vessel, and the substances in the blood that cause atherosclerosis (cholesterol, etc.) will move from the high-pressure area to the low-pressure area, i.e., from the vessel wall to the center of the vessel, reducing the atherosclerosis caused by excessive deposition of cholesterol and other substances on the vessel wall [27].

Figure 5 shows that during the diastolic phase, the flow velocity of the bifurcated vessels is low, the blood flow is disturbed, and the phenomenon of vortex and reflux is obvious, which tends to cause insufficient blood flow dynamics, and cholesterol and other substances in the blood will be stagnant and deposited in the vortex area, providing conditions for early atherosclerosis [28]. After sand therapy, the blood flow velocity increases and the vortex and reflux phenomenon decreases, and Fig. 5D shows that at the middle of the diastolic phase, the vortex disappears at the vascular bifurcation, which is conducive to reducing the accumulation of cholesterol and other substances, thus reducing the risk of atherosclerosis generation.

Figure 6 shows that the high viscosity of erythrocytes in the central region of the vessel is due to the fact that blood is a shear-thinning fluid with low viscosity at the wall surface velocity gradient and high viscosity at the center with a small velocity gradient [29]. The viscosity of erythrocytes in the central vascular region decreases after sand therapy, and according to Jung et al. the high viscosity zone tends to impede the flow of erythrocytes, causing them to migrate preferentially through the boundary layer, and the decrease in the viscosity of erythrocytes in the central vascular region after sand therapy facilitates the migration of erythrocytes from the vascular center, reducing the possibility of overcrowding of erythrocytes in the vessel wall leading to the formation of atherosclerosis [30].

Fig. 8.

Wall shear stress distribution of popliteal artery.

Fig. 9.

Comparison of TAWSS before and after sand therapy.

Fig. 10.

Comparison of OSI before and after sand therapy.

Fig. 11.

Comparison of RRT before and after sand therapy.

Figure 8 illustrates that the overall wall shear stress of the vessels experienced an increase subsequent to sand therapy. Atherosclerosis predominantly occurred in the low wall shear stress area as well as the high oscillatory wall shear stress area, where the vessels bend and bifurcate. This is because the low wall shear stress area can trigger intimal hyperplasia, causing harm to the vascular intima, and ultimately leading to the development of atherosclerosis [31]. In vitro studies [32,33] have also indicated that a shear stress of 1.0–1.5 Pa promotes the expression of arterioprotective endothelial genes, whereas shear stress below 0.4 Pa leads to an atherosclerotic phenotype. Moreover, following sand therapy, low wall shear stress levels increased from 0.13 Pa to 1.22 Pa, which may play a significant role in reducing vascular endothelial damage, promoting the expression of arterioprotective endothelial genes, and mitigating the occurrence of atherosclerosis.

Figure 9 shows that the maximum TAWSS value decreased after sand therapy and the area of high TAWSS decreased, which indicates that sand therapy reduces the occurrence of atheromatous plaques. Figure 10 shows that the OSI values decreased after sand therapy. As can be seen in Fig. 11, RRT values are reduced after sand therapy. Sand therapy decreases the probability of popliteal artery thrombosis by reducing the functional deficiency of endothelial cells and altering the cytoarchitectural cyclic stress.

5. Conclusion

This paper presents a computational study of flow-solid coupling using a two-phase flow-Casson model in a realistic vascular model of the human popliteal artery. Haemodynamic parameters were obtained before and after sand therapy, and conclusions were drawn that provide a theoretical foundation for the biomechanics of atherosclerosis prevention.

The increase of blood axial velocity after sand therapy can promote the movement of cholesterol and other substances toward the center of blood vessels and reduce the risk of the development of atherosclerosis caused by excessive deposition of cholesterol and other substances on the surface of blood vessel walls.

The reduction of vortex flow at the bifurcated blood vessels after sand therapy helps to reduce the accumulation of cholesterol and other substances.

The reduction of erythrocyte viscosity in the center of the vessel after sand therapy facilitates the migration of erythrocytes from the center of the vessel.

The increase of low wall shear stress after sand therapy can reduce the occurrence of atherosclerosis.

The decrease of time-averaged wall shear stress, shear oscillation index and relative retention time after sand therapy can reduce the occurrence of thrombosis.

Footnotes

Conflict of interest

None to report.

Funding

This research was financially supported by the National Natural Science Foundation of China (no. 31460245).

References

Bugebi

Adayi

, Hemodynamic effects of VMC sand therapy on femoral artery branches, Medical Biomechanics 31(5) (2016), 431–436.

L. Liu, Dilinaer-Mahmuti, R. Fu et al., Experimental study on the effect of indoor sand therapy on human femoral artery hemodynamics, Journal of Xinjiang Medical University 37(1) (2014), 29–31.

Qiao

Y.Q.

R.C.C.

and Zhang

L.H.

, Effect of sand therapy on the pulsatile dynamics of femoral artery bifurcation vessels with different degrees of stenosis, Chinese Tissue Engineering Research 24(8) (2020), 1218–1224.

and Yang

, Effects of sand therapy on fluid dynamics of blood in the popliteal artery based on the euler multiphase modle, Journal of Mechanics in Medicine and Biology 22(10) (2022), 2250044.

Cardiovascular Disease Prevention and Control Research Center, Ministry of Health. Summary of the Chinese cardiovascular disease report 2010, Journal of Nursing Management 1 (2012), 14–15.

Redberg

R.F.

Benjamin

E.J.

Bittner

, ACCF/AHA 2009 performance measures for primary prevention of cardiovascular disease in adults, Journal of the American College of Cardiology 54(14) (2009), 1364–1405.

Abbasian

Shams

Valizadeh

, Effects of different non-Newtonian models on unsteady blood flow hemodynamics in patient-specific arterial models with in-vivo validation, Computer Methods and Programs in Biomedicine 186 (2020), 105185.

Gidaspow

, Multiphase flow and fluidization: continuum and kinetic theory description, Journal of Non-newtonian Fluid Mechanics 55(2) (1994), 207–208.

Anderson

T.B.

and Jackson

, A fluid mechanical description of fluidized beds, I & EC Fundamentals 6(4) (1967), 527–539.

10.

Ling

Tang

and Liu

, Numerical investigation of two-phase non-Newtonian blood flow in bifurcate pulmonary arteries with a flow resistant using Eulerian multiphase model, Chemical Engineering Science 233 (2021), 116426.

11.

Jung

Lyczkowski

R.W.

Panchal

C.B.

, Multiphase hemodynamic simulation of pulsatile flow in a coronary artery, Journal of Biomechanics 39(11) (2006), 2064–2073.

12.

Gijsen

F.J.

Vosse

and Janssen

J.D.

, The influence of the non-Newtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model, Journal of Biomechanics 32(7) (1999), 705–713.

13.

Dill

D.B.

and Costill

D.L.

, Calculation of percentage changes in volumes of blood, plasma, and red cells in dehydration, Journal of Applied Physiology 37(2) (1974), 247–248.

14.

Yong

H.K.

Vande vord

P.J.

and Lee

J.S.

, Multiphase non-Newtonian effects on pulsatile hemodynamics in a coronary artery, International Journal for Numerical Methods in Fluids 58 (2010), 803–825.

15.

Bird

R.B.

Armstrong

R.C.

Hassager

, Dynamics of polymeric liquids, vols. 1 and 2, Physics Today 31(2) (1978), 54–57.

16.

Lightfoot

E.N.

, Transport Phenomena and Living Systems: Biomedical Aspects of Momentum and Mass Transport, University of Nottingham, 1985.

17.

ANSYS FLUENT User's Guide 14.0, ANSYS, Inc., 2011.

18.

Jung

and Hassanein

, Three-phase CFD analytical modeling of blood flow, Medical Engineering & Physics 30(1) (2008), 91–103.

19.

Kim

Y.H.

Vandevord

P.J.

and Lee

J.S.

, Multiphase non-Newtonian effects on pulsatile hemodynamics in a coronary artery, International Journal for Numerical Methods in Fluids 58(7) (2008), 803–825.

20.

Morales

H.G.

Larrabide

Geers

A.J.

, Newtonian and non-Newtonian blood flow in coiled cerebral aneurysms, Journal of Biomechanics 46(13) (2013), 2158–2164.

21.

Gao

Watanabe

and Matsuzawa

, Stress analysis in a layered aortic arch model under pulsatile blood flow, Biomedical Engineering Online 5(1) (2006), 1–11.

22.

Qiao

Yin

and Chu

, Numerical simulation of fluid-structure interaction in bypassed DeBakey III aortic dissection, Computer Methods in Biomechanics and Biomedical Engineering 18(11) (2015), 1173–1180.

23.

Cong

Zhao

, Hemodynamic analysis of the right coronary artery originating from the left coronary sinus, Medical Biomechanics 35(3) (2020), 284–288.

24.

Luo

, Multi-Component Multiphase Hemodynamic Study in Stenotic Vessels, Nanchang University, 2017.

25.

Chaichana

Sun

and Jewkes

, Computation of hemodynamics in the left coronary artery with variable angulations, Journal of Biomechanics 44(10) (2011), 1869–1878.

26.

Bugebi

Adayi

, Study on the pulsatile hemodynamics of bidirectional flow-solid coupling in the femoral artery under vitamin therapy, Journal of Interventional Radiology 26(1) (2017), 49–54.

27.

Chandran

K.B.

Rittgers

S.E.

and Yoganathan

A.P.

, Human circulatory system, in: Biofluid Mechanics, 2015.

28.

Huo

Wischgoll

and Kassab

G.S.

, Flow patterns in three-dimensional porcine epicardial coronary arterial tree, The American Journal of Physiology-Heart and Circulatory Physiology 293(5) (2007), H2959–H2970.

29.

Zhang

and Xie

H.W.

, Comparative analysis of blood flow simulation methods, Practice and Understanding of Mathematics 48(1) (2018), 162–167.

30.

Zhao

Ping

, Computational fluid dynamics simulation of time-resolved blood flow in Budd–Chiari syndrome with inferior vena cava stenosis and its implication for postoperative efficacy assessment, Clinical Biomechanics 82 (2021), 105256.

31.

Caro

C.G.

Pedley

T.J.

Schroter

R.C.

The Mechanics of Circulation, Oxford University Press, 1978.

32.

Malek

A.M.

Alper

S.L.

and Izumo

, Hemodynamic shear stress and its role in atherosclerosis, Journal of the American Medical Association 282(21) (1999), 2035–2042.

33.

Topper

J.N.

Cai

Falb

, Identification of vascular endothelial genes differentially responsive to fluid mechanical stimuli: Cyclooxygenase-2, manganese superoxide dismutase, and endothelial cell nitric oxide synthase are selectively up-regulated by steady laminar shear stress, Proceedings of the National Academy of Sciences of the United States of America 93(19) (1996), 10417–10422.