Abstract
An experiment-based computational study that helps analyze blood flow behavior and wall shear stress (WSS) distribution is reported in this work. Large scale numerical analysis of hemodynamics in swine-specific stenosed carotid artery based on in vivo surgery is presented. A pressure stabilized domain decomposition method is used to symmetrize the linear systems of Navier-Stokes equations and the convection-diffusion equation. A numerical expression of swine blood flow and a detailed swine carotid vessel model with stenosis are newly proposed, and the empirical function of WSS was validated for the swine model. Two wall models, a rigid and another elastic, are compared in precisely modelling for pathological analysis of vascular disease like carotid atherosclerosis and hemangioma. The flexible wall performs better in representing experimental conditions while the stern wall is much more efficient. Numerical results show that the stenosis has a great influence on the behavior and characters of blood and its subsequent affect the WSS of the vessel; further details show how stenosis affect the distribution and magnitude of wall shear stress in an artery which lay a foundation for further medical study.
Keywords
Introduction
Shreds of evidence show that changes in arterial hemodynamic characteristics, especially in wall shear stress (WSS), play a significant role in the initiation and development of vascular diseases; it may cause impedance increases, 1 pattern changes,2,3 and intracranial atherosclerotic disease.4,5 Abnormal blood flow gives rise to various vascular disorders, such as atherosclerosis,6–8 plaque deposit,9,10 and accumulation of lipids. 11 Stenosis is a common biological phenomenon that is related to many disadvantages to the cerebral 12 and cardiac health. 13 Many previous studies based on observed data or human body test data have been conducted and achieved valuable findings in the effect of stenosis.2,5 However, the further understanding of its influence in blood vessel coupling is restricted due to the great danger and moral concerns of human experiments.
Consequently, most recent studies tend to use in vivo testing for non-human animals in experiments that seek to control the variables that affect the actions or biological system under study.14–21 The use of the mini-swine as a nonrodent species in research has continued to expand for over a decade, and they are becoming routinely used both in experimental pharmacology and as a therapeutic model for human diseases. 21 Chinese mini-swines, which has the most similar blood and vessel characteristics comparing to human beings, has been widely used in experimental hemodynamics study.8,19,22 However, in vivo animal experiments suffer from many factors in the pathological analysis, such as the long experimental period, high cost, and the difficulty in replicating due to the non-reproducibility of operation conditions.8,18 Computational fluid dynamic (CFD) approach uses the geometry and mechanical properties of the vasculature and the principles of conservation of mass and momentum to obtain the blood flow-pressure relation; with the development computational medical technology, CFD simulation becomes a prevalent method in studying hemodynamics for its high efficiency and low cost. 23 For the primary large arteries of the human vasculature, Randles et al. proposed a parallel simulation and with scaling efficiency or 96% on 1,572,864 cores.
Meanwhile, the difficulties in applying CFD simulation toward clinical applications show up. 24 As in this study, one challenge is the exact setting of parameters, which is utilized as material properties or initial boundary conditions for the simulation and can be further employed in the validation of numerical modeling. Accurate analysis of the hemodynamics of Chinese mini-swines requires high precision reconstruction of the vessel and corresponding blood flow characteristics as input parameters for the simulation process. The entire experimental process includes many stages, like animal culture, surgical operation, Doppler ultrasonography, and three-dimensional digital subtraction angiography. Due to the complexity of the vascular architecture and precise control of animal experiment mentioned above, much previous research employed empirical equations such as Womersley model25,26 and Hughes model 27 to simulate pulsatile blood flow behaviour in an artery; Leng et al. proposed a Eulerian-Lagrangian model for fluid-rotating structure interaction problems based on the arbitrary Lagrangian-Eulerian approach. Although this approximation may be reasonable under certain assumptions, the nonlinear effects must be considered when the Reynolds and Womersley numbers are higher 28 ; moreover, blood flow behaviour varies from species to species; to the best of our knowledge, seldom about the blood characteristics was documented for Chinese mini-swines. Hence, for accurate prediction of transient pressure and flow at any position for the swine model, a higher-precision description of the blood behavior of the Chinese mini-swines should be employed instead of a generalized empirical function.
On the other hand, another major problem in blood vessel modeling is that often, in reality, blood vessel wall materials are typically elastic and hence the vessel wall structure and shape tend to change due to the blood flow rushes in a rough or turbulent manner,29,30 especially in abnormality situations like narrowing or partially blocked. 30 The modeling such a coupled Multiphysics phenomena is typically difficult to achieve because of both modeling complexity and the increased computational time requirement. Nevertheless, although studies that assume blood vessel walls to be a rigid structure are valuable in terms of saving computational time, the validity of such simulations concerning the real blood flow behavior is questionable due to the lack of vessel deformation. 31 The development of a coupled flexible blood vessel to the blood flow dynamic model is becoming more feasible and should be introduced to improve the accuracy of the related state of the art models.
Besides, other problems need to be improved in a patient-specific simulation when involving the swine model developed in this work. Existing blood flow functions mostly represent the human blood flow profile and no precise numerical models being put forward for swine blood flow. Although swine models share similar blood flow characters with the human, there still exits differences in amplitude and heartbeat cycle and that a more precise relationship needs to be established in swine-model-based research. Moreover, since it is notably challenging to obtain the velocity gradient at locations very close to the vascular wall, researchers generally use an empirical function to calculate WSS in vivo experiments. However, this empirical WSS function has not been verified with experimental results, especially when applying to the mini-swine model. Its reliability needs further support from the numerical simulation.
A more realistic and high-fidelity hemodynamic numerical simulation of vascular stenosis is expected in this work by overcoming the difficulties mentioned above. By using the in vivo experiment and three-dimension digital subtraction angiography (DSA), a patient-specific computational model that describes the detailed anatomical configuration, for example, the asymmetric vessel wall enlargement after the stenosis will be precisely reconstructed. Based on the data obtained from the Doppler ultrasonography, a high-order numerical expression for the amplitude and period of Yucatan swines blood flow is expected. Two wall description models, one rigid and the other elastic, are compared in this study. Materials with various parameters are applied to find out the suitable simulation method of carotid artery stenosis in the swine model. Thickness is considered in the elastic wall model and thus, the Fluid-Solid Interaction (FSI) model is applicable for vascular simulation. Aiming to validate the feasibility of a patient-specific numerical simulation method and to determine the fundamental relationship between WSS and stenosis, which may lay a foundation for a further medical study, the data measured from carotid arteries of swine are employed to simulate the impulse of blood. The flow pattern, including the pressure, WSS, and velocity distribution in the artery, is analyzed; and the correlation with a vascular disease like angiogenesis and further stenosis will be examined.
The following of this paper is organized as follows: the experimental and numerical methods are discussed and elaborated in Section 2. Section 3 shows the experimental and simulation results, as well as the comparison for the validity of two numerical models. In Section 4, the conclusions of the current research are drawn, and its applications are discussed.
Methods
Swine model of carotid atherosclerosis
The experimental protocol is approved by the institutional animal research committee of Sun Yat-sen University, under policies set by the National Institutes of Health guidelines. Chinese mini-swines, 10–12 weeks of age and 15–20 kg body weight, are used to develop carotid atherosclerosis using the combination of partial surgical ligation and hyperlipidemia as described previously.18,31 Doppler ultrasonography is used to measure the diameter and blood flow velocity. Bilateral carotid arteries in each swine are partially ligated to create approximately 70% stenosis. Plasma cholesterol levels and blood viscosity are measured. WSS in the carotid arteries is calculated as described previously. 14 Three-dimensional rotational digital subtraction angiography (DSA) is performed to measure the diameter of the carotid artery and the degree of stenosis at the time of post ligation.
3D reconstruction
A DSA device from SIEMENS with a resolution of 0.1 mm is employed in this research to understand more details on the vertical axis. The target of the in vivo experiment is Chinese mini-swine, and the carotid artery of swine is stenosed by surgery. After injecting the contrast-medium into the vessel of Chinese mini-swines in vivo, a body scan by computed tomography (CT) is performed. A post-image processing is utilized to capture the 3D structure of vessels after the scanning. Syngo Neuro DSA CT, a dedicated post-processing application running on the SIEMENS workstation, which allows the removal of bone structures from CT-Angiography (CTA) data sets for improved visualization of the cerebral vasculature, gives images of the reconstruction as is shown in the reconstruction part of Figure 1.

In vivo DSA image and reconstruction for 3D modeling.
The algorithm of syngo Neuro DSA CT works fully automatic and makes this application easy to use; however, the workstation does not support the lossless output of the reconstructed model, and only black-white patterns of CT scanning are available; a customized 3D reconstruction according to the dimensions of vessels measured by CT image is implemented before the meshing, see the DSA part of Figure 1. In this work, the black-white patterns are imported into a post-processing software AMIRA and reconstruct a 3D vessel model. Due to the influence of the flaws caused by the scanning process, a growing regional method, thresholding method, and manual segmentation method based on grayscale value are applied to refine the quality of the model surface in AMIRA.
Numerical modeling
The fluid is assumed to be Newtonian, and parameters of material property are set to the values measured by the in vivo experiment. Since the subtle difference of the property parameters of the vessel is insignificant, average values are taken for the vessel properties. In the transient structural part, the vascular wall density is set to 1062
where
with the Kronecker delta
Hydrodynamic boundary conditions of constant pressure for the outlet and value-specified velocity for inlet were employed in this study. In order to simulate the realistic behavior of blood, we measured the periodic velocity change of blood in the target artery by Doppler ultrasonography. The inlet velocity was fitted according to the experimental data. The expression is a function that simulates a real period impulse of vascular flow.
where v is the velocity of blood and t represents time, see Figure 1. The coefficients
Coefficient settings.
The inlet velocity in a cardiac cycle is shown at the bottom of Figure 1; its range starts from 7.23 to 23.46 cm/s with a period of 0.8 s. A user-defined function (UDF) is applied in Fluent to simulate this velocity alternation periodically. Compared to existing pulsatile blood flow functions, the image-based velocity profile restores the maximum level of the trueness of real blood flow characteristics. This intrinsic correspondence between our patient-specific vessel model and velocity profile, rather than empirically drawing vessel wall prototype and velocity expressions from existing models, assures the best reality of the simulation and reliability of the results.
The boundary condition between the blood and the vascular wall is no slip. Meanwhile, the vessel was defined as a rigid and elastic wall to compare. For the rigid wall, the vessel is a stationary wall with no-slip conditions. While for the elastic wall, since physiological tissue has constrained to the vessel movement, apparent elastic modulus should be applied in the study. Given the tissue constraint and vessel’s elastic modulus, it is reasonable to define as follows, Young’s modulus of the wall that may deform equals to 4.66 MPa, Poisson’s ratio is 0.45, and thickness is 0.4 mm.
Meshing and domain decomposition
In order to achieve a better investigation of the hemodynamics of a carotid artery with stenosis inside a Chinese mini-swine, two numerical models for vessel wall function are built and compared in this work, namely the rigid wall model and the elastic wall model; an elastic wall with thickness is considered for the simulation of vessel wall, and FSI is employed to simulate the interaction between blood and vascular wall.
Aiming at achieving a reliable patient-specific numerical simulation method and discovering the fundamental relationship between WSS and stenosis, we employed the data measured from carotid arteries of swine to simulate the impulse of blood. An FSI model is introduced to the simulation to consider the thickness of the vascular wall in the elastic wall model.
In the swine model of carotid atherosclerosis, the average diameter of the healthy artery is 3.15 mm, and the smallest diameter of the stenosed part is 1.4 mm, which is approximately 44.4% of the normal one. The vascular wall function of the artery was established based on the information collected from the experiment. Although the thickness of the vascular wall varies slightly from the inlet to the outlet in the real situation, an average value of 0.4 mm is set for it in this work as the influence due to thickness variation can be neglected and the vascular is regarded as isotropy in this research.
The fluid and the solid domain must be processed separately for the simulation of the elastic wall model as required in the FSI method. For each case of meshing, a computational mesh built of tetrahedral elements, condensed in the region of stenosis. Additionally, inflation elements are introduced in the vicinity of the wall to reflect flow manners at vascular walls more precisely. A sample mesh is shown in the mesh part of Figure 1. The mesh independence tests have not yielded any significant differences that could affect the computation results. Therefore, due to time-consuming transient simulations, a middle-size density of 0.2 mm has been employed. The fluid domain has 800,000 elements, and the solid domain has 500,000 elements.
A domain decomposition method, 32 effective in large scale computation, is employed to show more details of dynamic behavior of the model, especially for the fluid structure coupling in the elastic wall model. The mesh before and after domain decomposition is given by Figure 2. Unstructured grids, which is derived from smoothing, layering and remeshing method, have excellent performance when confronting mesh deformation and are utilized in this work. Dynamic mesh enables the mesh to be reconstructed easily with less occurrence of errors and negative volumes. Text Command List (TUI) commands are utilized to open a smoothing-parameters setting to keep the dynamic mesh working smoothly.

Meshing and domain composition.
Results and discussions
The solutions are computed for five cardiac cycles at the beginning and the results show there is no significant discrepancy between the first cycle and the remaining. For the efficiency of the computation, two cycles are good enough to approach steady results. The average computation cost for the rigid wall model is 1.5 h by parallel computing with 16 nodes involved. However, for the elastic wall model, the average computation cost is 10.5 h.
To make sure the approach is reliable the results between the experimental data and the simulation results for both systole and diastole in a cardiac cycle are compared. Velocity vector, WSS and vorticity magnitude distribution are three main characters that we concern about. We also monitor the average velocity of the cross-section along the artery. By analyzing these characters, the relation between WSS and the corresponding disease can be studied.
Validation
The WSS obtained from an empirical function is compared with that extracted from numerical results. The difference between the two models at both vasoconstriction and vasodilation is also involved, as shown in Figure 3. The experimental value of WSS is calculated by
where

Comparison of WSS distribution (upper: rigid wall; lower: elastic wall).
According to Figure 3, WSS from experiment shares similar values and trends with numerical results, especially for elastic wall model: WSS increases gradually from location L1 to L7 and reaches the maximum at location L8; at the post-stenotic area, WSS decreases stepwise from location L8 to L15. The plot of experiment values always stays between vasodilation and vasoconstriction curves in the elastic wall case. It is also true for the rigid wall model except at location L4 to L7 during vasodilation, where WSS is slightly more significant than the experiment values. When comparing two different models, WSS values of the elastic wall model are slightly closer to the experiment values. During the blood vessel vasodilation, the maximum WSS occurs at location L8 in the elastic wall case. However, the maximum WSS location shifts a little bit toward the inlet direction in the rigid wall case.
Meanwhile, WSS in both elastic wall model and the experiment increases gently before the stenosis. It exists a minor variation before the stenosis in the rigid wall case. During the blood vessel vasoconstriction, both models work equally well; maximum WSS always occurs at location L8. Since WSS is a measure of force per unit area exerted by the blood flow in the direction that tangent to the vessel wall, it is most likely to have the greatest WSS occurring at or near the stenosis, where the flow velocity is expected to be the maximum, and remaining constant or displaying little variation when it is distant from stenosis. Along the several sections that immediately follows the stenosis, the velocity drops gradually as shown in Figure 4, which in theory indicates a stepwise decrease in WSS. Our observation from Figure 3 agrees with this theoretical results in all cases.

Comparison of the velocity profile of the two models.
Velocity distribution
Fifteen monitoring points are defined along the carotid artery during a Doppler ultrasound imaging process to obtain blood flow velocity. Fifteen cross-sectional planes at corresponding locations are created in the numerical model (Figure 1). The averaged velocity at each cross-section is calculated and compared between the two models (Figure 4). In general, the velocity increases gradually in the first half of the artery. It then increases much more rapidly and reaches the maximum at the stenosis. In the second half of the artery, the velocity decreases slowly and becomes constant at the end. Compared to two models, blood flow velocity profiles are similar in the first half of the artery, while the velocity decreases much more slowly in the elastic wall mode in the post-stenotic region. The velocities are very close to each other at the stenosis at most time points between two groups. However, the velocity in the elastic wall model is slightly lower than in the rigid wall model at the stepwise decreasing period (Figure 5).

Comparison of velocity at stenosis of two models.
The rate of stenosis at the beginning of each cycle (blood vessel vasoconstriction) quickly reaches the maximum and then presents a ladder-like descending trend and reaches a minimum at the end of the cycle (blood vessel vasodilatation). When the blood flow velocity rises, the results of the elastic wall model and the rigid wall model are similar; when the blood flow velocity begins to decrease, the velocity of the elastic wall model is slightly lower than the velocity of the rigid wall model. Due to the elasticity of the blood vessel wall model, the shape of the blood vessel changes as the blood flow velocity changes. That is, when the blood vessel vasodilatation happens, due to the elasticity of the vessel wall, the inner diameter of the wall of the elastic wall model becomes more extensive than that of the solid wall model; therefore, according to the law of conservation of mass, the velocity will decrease.
WSS and streamlines
As it can be seen in Figure 6, WSS changes along with the changes in blood flow velocity periodically; WSS presents different distribution at different time points within one heartbeat cycle. It has the most extensive distribution at systole while less distribution in the following time points when only looking at the WSS higher than 2.000e + 01 Pa. Moreover, from systole toward diastole, WSS decreases gradually in the region before stenosis, this confirms our results shown in Figure 5. The same situation also happens in the post-stenotic region, especially at the right-hand side vessel wall in both cases. At each time point, maximum WSS always occurs at the stenosis and pre-stenotic regions. However, it is noteworthy that although WSS reaches its local maximum at systole in the post-stenotic area, the streamline is quite uniform comparing to the other time points. In the following time points, the streamline shows vortices. One possible explanation could be: as the velocity component that is parallel to the vessel wall decreases, the blood flow can travel more distance along the direction that is perpendicular to the vessel wall, which further gives rise to the vortices we observed. The decrease in the velocity component tangent to the vessel wall matches the decrease in WSS in the post-stenotic area.

Comparison of WSS (upper) and streamline (lower) of two models at different stages in a cardiac cycle (R: rigid wall; E: elastic wall).
When comparing rigid wall and elastic wall cases, WSS distribution is quite similar in general. However, the maximum WSS area shifts slightly toward the inlet direction in the elastic wall case, and the area is also smaller comparing to the rigid wall. As for the streamline, both cases have uniform streamlines at the systole. The streamline of the elastic wall becomes disturbed faster than the rigid wall when it comes from systole toward diastole. Meanwhile, the streamline of the elastic wall also returns to a uniform pattern faster than the rigid wall case. It is noteworthy that, vessel wall enlargement occurs in those vortex regions. The extra room could be developed because of the stenosis in order to reduce the WSS as a physiological effect.
Comparison near the stenosis
Based on the above data and the results of the chart shown in Figure 7, both cases have observed significant vortices in the downstream region at the narrowing, but the locations are different in the case of rigid walls and elastic walls. This is consistent with what we observed in Figure 6: vortices shown by streamline in post-stenosis region. Maximum vorticity occurs at the stenosis in all cases, and it decreases along the outlet direction gradually. During vasoconstriction, rigid wall and elastic wall cases share similar vorticity distribution: the post-stenotic region has higher vorticity than the pre-stenotic region; maximum vorticity concentrates close to the vessel wall at stenosis and the post-stenotic area. During vasodilation, the post-stenotic region still has more significant vorticity than the pre-stenotic region. However, vorticity distributed much wider than vasoconstriction in both rigid and elastic wall cases comparing to vasoconstriction. What is more, a rigid wall always has mildly less vorticity in the middle of the pre-stenotic region comparing to the elastic wall, and this phenomenon is more distinct during the vasodilation.

Comparison of WSS (left) and vorticity (right) profiles (R: rigid wall; E: elastic wall).
Noteworthy, the subplot on the left shows a narrow band that indicates a sudden decrease in WSS distribution. This “trough” in WSS distribution is consistent with experiment results previously reported by Sadeghi et al. 33 The sudden decrease of WSS is due to the increase of wall circumferential stress (WCS). Flow immediate after the stenosis exerts a significant greater force perpendicular on the boundary, leading to a decrease in stress vector that tangent to the wall. However, this abrupt change in WSS disappears very quickly with the expansion of the vessel, resulting in a relatively high WSS and the following stepwise decrease.
Conclusion
In this study, blood flow and its dynamic behavior are simulated by large scale simulation. The in vivo experiment successfully replicates anatomical abnormality that could have occurred in a swine. With three-dimension DSA, a patient-specific computational model that describes the detailed anatomical configuration was reconstructed, and a high-order numerical expression for swine blood flow from ultrasonography was proposed, which presents precise amplitude and a period of swine blood flow. Measurements such as blood flow velocity and WSS obtained from in vivo tests provide precise boundary conditions and other information, which was utilized as the inputs and the benchmarks for model evaluation. The modeling was validated by comparing the results from simulation and experiment.
A comparison of a rigid wall and an elastic wall model was made in this work, and a domain decomposition method is employed to ensure the accuracy by large scale simulation. The dynamic behavior is more evident in the elastic wall than the rigid wall according to the distribution of WSS and streamline; however, not so much difference was observed in the distribution of magnitude. The rigid wall model is good enough to represent experiment results if time-efficiency and computational cost are considered. In future applications, one could use a rigid wall model to obtain a convincing result when involving a large sample size.
Physiologically, our model and simulation can help clinics to better obtain a full image of hemodynamics for a specific experiment subject, ranging from systole to diastole. A detailed spatio-temporal hemodynamic profile including velocity and WSS distribution, streamline profile as well as vorticity distribution can be obtained and studied prior to any further intervention. To be more specific, significant WSS not only existed even before the stenosis but also at part of the post-stenosis area and its magnitude changes across time within a complete cardiac cycle. Given the advantage of domain decomposition in large scale simulation, instructive simulation results for clinics can be approached in a short period of time.
Footnotes
Handling Editor: James Baldwin
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was funded by the National Key R&D Program for International Collaboration and for HPC, grant numbers 2018YFE9103900 and 2016YFB0200603. The Natural Science Foundation of China (NSFC), grant number 11972384 and Guangdong MEPP Fund, grant number GDOE[2019] A01, and a grant from Guangzhou Science and Technology Program, numbered 201704030089, also supported this work.
