Impact of Modelling Surface Roughness in an Arterial Stenosis

Arterial stenosis is a problem of immediate significance, as cardiovascular disease is the number one leading cause of death worldwide. Generally, the study of stenotic flow assumes a smooth, curved stenosis and artery. However, the real situation is unlikely to present an infinitely smooth-surfaced arterial stenosis. Here, the impact of surface roughness on the flow in an arterial stenosis was studied via a computational fluid dynamics analysis. A patient-specific geometry with a smooth surface was reconstructed, and a partially rough model was built by artificially adding random roughness only on the stenotic region of the smooth model. It was found that the flow was oscillatory downstream of the stenosis in the models. A slightly lower velocity near the wall and more oscillatory flows were observed due to the presence of the roughness in the stenotic region. However, the pressure distributions did not vary significantly between the smooth and rough models. The differences in the wall shear metrics were slight in the stenotic region and became larger in the downstream region of the models.


Introduction
Cardiovascular disease and associated arterial stenosis is the number one leading cause of death worldwide, contributing to 30% of all annual fatalities [1]. Stenosis development is strongly linked to spatial and temporal variations in wall shear stress (WSS) [2], which can be caused by flow separation, reversed flow and unsteadiness [3]. Turbulent flow results in significant temporal and spatial variations in velocity, which cause large WSS fluctuations; this is complicated by the pulsatile flow presented in arterial systems, where turbulent and transitional flows may appear for some parts of the cycle. Pulsatile flow in pipes features different flow characteristics from those of the steady state, as the velocity across the diameter of the pipe can vary considerably, and therefore turbulence can occur near the walls while the core of the flow remains laminar [4]. Typically, based on a nonconstricted region diameter, the peak Reynolds number is around 1000 for these types of flows; however, actual velocities through the stenosis cause the local Reynolds numbers to be significantly larger in this region. The critical Reynolds numbers for the transition to turbulence decrease with smaller roughness-wavelength-to-tube-diameter ratios, and it is well known that the inner layer of a turbulent boundary layer is significantly altered due to the surface roughness [5].
There exists a significant body of work regarding stenotic flow [6][7][8][9][10]; however, the majority of studies assume a smooth, curved (usually axisymmetric) stenosis and artery. However, reality is unlikely to present an infinitely smooth-surfaced arterial stenosis. In the early stage of atherosclerosis, the luminal surface of an arterial wall becomes rough as a result of endothelial damage [11], with a roughness of the order of the cell height, i.e., from 10-20 µm up to a few hundred µm [12]. A stenosed blood vessel may contain rough irregularities that resemble small valleys and ridges [13], and irregular plaque deposits are more likely to be encountered than smooth stenoses in pathological conditions [14].
Stenosis surface irregularity can lead to complex flow patterns [13] and therefore a variation in the local WSS, although the pressure drop is unaffected as the geometric effect of the stenosis is much more significant [13,15]. A smooth surface can also result in a reduction of as much as 80% in the rate of hemolysis [16].
Mustapha et al. [17] modelled an irregular stenosis couple and noted that the multiple stenoses, as well as their surface roughness, prompted early separation of flow. Bhaganagar [5] conducted a direct numerical simulation study with four types of roughness in a channel flow, representative of four types of arterial roughness. The study revealed that the turbulence generated due to stenosis is strongly dependent on the shape of this stenosis, and the turbulent structures exhibit characteristic structures and mechanisms very specific to the different roughness types. In particular, the ridge/valley representation resulted in strong near-wall velocity gradients with roll-up of the shear layer near the peak locations and recirculating regions near the valley locations. The peak turbulence production was also higher compared to the other three types. WSS effects were not reported.
Bark Jr. and Ku [18] used a random roughness with an amplitude of 200 µm to simulate surface roughness. In their computational model (using a laminar Navier-Stokes solution), it was found that shear-rate peaks developed along the hills with small recirculation eddies in the valleys of the rough surface, causing platelets to experience varying shear stresses near the apex of a stenosis. The irregular thrombus surface resulted in focal areas of increased shear rate relative to a smooth surface. However, Xiong and Taylor [19] found that surface roughness had only a minor effect on the WSS, although the local flow was affected.
A complex pattern of vessel roughness is difficult to describe in an essentially qualitative angiographic report, prompting Cho et al. [20] and Back et al. [21] to take a silicone casting from a cadaveric coronary artery (with significant but non-obstructive diffuse atherosclerotic disease) to obtain the detailed geometric variation up to an order of accuracy of 0.1 mm, and these data have been used in subsequent numerical studies [22]. Stroud et al. [14] used high-resolution imaging of surgically removed plaque specimens from carotid and coronary arteries to define a range of characteristics for arterial stenosis models. These were then modelled with a laminar Navier-Stokes solver, to evaluate the isolated effects of curvature and surface irregularity. These features were observed to have substantial effects on hemodynamics descriptions such as WSS. Fully three-dimensional (3D) irregularities significantly increased the shear stress at the wall, and an oscillating pattern of WSS was observed, even far from the stenosis.
Although the impact of stenotic surface irregularity has been studied, a hemodynamic analysis for a patient-specific model is still lacking. In this paper, a patient-specific model with a smooth surface was clinically obtained and computed, and a rough model was created by manually adding roughness only in the stenotic region of the smooth model. We aimed to investigate the impact of the surface roughness of the stenosis on the blood fluid dynamics. Comparisons are presented between the smooth and rough models via computational fluid dynamics (CFD) studies.

Materials and Methods
A patient-specific model with a smooth surface of a left anterior descending coronary artery was reconstructed from intravascular ultrasound images, which were obtained in a catheterization laboratory at 15 frames per second (Axiom Artis, Siemens, Germany) ( Figure 1A). The diameter at the inlet was 2.97 mm, while the diameter at the outlet was 0.65 mm. The area at the narrowing region was about 0.57 mm 2 , and the area reduction in this case was 79% (percentage diameter reduction of about 54%). The lesion length was 13.14 mm, and the total length was 63.24 mm. Random roughness was constructed manually, only in the stenotic region, using SolidWorks 2019 ( Figure 1B). Small faces in the stenotic region were recognized and segmented automatically after importing the geometry into the software. Ridges and valleys were created along the direction of the vector of each sketch face. The lengths and widths of the ridge and valley were varied and determined based on the curvature of the selected face, but the roughness height remained the same at 50 µm. The surfaces of the upstream and downstream segments in the rough model were unchanged and smooth, as shown in Figure 1B. Finite elements with five prism layers were created using Ansys ICEM, and the minimum mesh quality was 0.37. To capture the flow variations in the rough model, a finer mesh was created in the stenotic segment, of the order of the roughness height, in contrast to the upstream and downstream regions (zoomed-in view in Figure 1). manually, only in the stenotic region, using SolidWorks 2019 ( Figure 1B). Small faces in the stenotic region were recognized and segmented automatically after importing the geometry into the software. Ridges and valleys were created along the direction of the vector of each sketch face. The lengths and widths of the ridge and valley were varied and determined based on the curvature of the selected face, but the roughness height remained the same at 50 µ m. The surfaces of the upstream and downstream segments in the rough model were unchanged and smooth, as shown in Figure 1B. Finite elements with five prism layers were created using Ansys ICEM, and the minimum mesh quality was 0.37.
To capture the flow variations in the rough model, a finer mesh was created in the stenotic segment, of the order of the roughness height, in contrast to the upstream and downstream regions (zoomed-in view in Figure 1). Simulations were performed in Ansys Fluent 2020 using a high-performance computing cluster (Gadi, National Computational Infrastructure, Australia). The fluid was assumed to be Newtonian and incompressible, with a density of 1050 kg/m 3 and a viscosity of 0.0035 Pa·s. These values are typical for arterial simulations [23,24]. The patient-specific pulsatile waveform of the velocity was prescribed at the inlet ( Figure 1C), while zero pressure was applied at the outlet. A fully-developed parabolic flow can be represented at the inlet by applying a user-defined function in the software [25]. The peak Reynolds number was around 330, and the flow was assumed to be laminar in this study. The artery was considered as a rigid conduit with a non-slip wall. The Reynolds-averaged Navier-Stokes equations were solved for the modelling of the blood flow, and the PISO (pressure-implicit with splitting of operators) scheme was used for the pressure-velocity coupling, with the accuracy set to be second order. To reach a Courant number lower than 1 and ensure numerical stability [26], the time-step size was set to be 2 × 10 −6 . The maximum number of iterations of each step was 20, and the convergence level was set at 1 × 10 −3 . Results were obtained and analyzed for numerical stability and accuracy after three cardiac cycles. The time-averaged wall shear stress (TAWSS), oscillatory shear index (OSI), relative residence time (RRT) and transverse wall shear stress (transWSS) were calculated by executing user-defined functions that are hooked into the iterative process. The definitions and descriptions of these shear metrics are given in Table 1. Table 1. Definitions and brief descriptions of four shear metrics in this study.

Equation Description
Time-averaged wall shear stress (TAWSS) Oscillatory variation in the WSS [27] Relative residence time (RRT) RRT ∼ 1 (1−2 * OSI) * TAWSS Near-wall flow stagnation [28] Transverse wall shear stress (transWSS) dt Multi-directionality of the flow field [29] A grid independency study was performed in the rough model to obtain a convergent and reliable result in the simulations. The meshing size in the model was around 6,700,000. The parameters for the grid generation were then decreased to half, and the total element size was increased to 9,800,000. The values of the area-weighted average WSS at the peak velocity were selected for the mesh independency study and compared after the computation. These values were 74.51 and 74.64 Pa in the models with a normal and a fine grid, respectively. The difference was observed to be 0.17% for the rough model. Therefore, the normal element grids were deemed to be appropriate for use in the study. Figure 2 presents the velocity streamlines at the peak of the cardiac cycle. The maximal velocity, which is larger than 8 m/s, occurs at the center of the stenosis. Oscillatory flow is observed in the region downstream of the stenosis in both models. To further visualize the flow contour near the ridges and valleys, zoomed-in views of the slices across the stenotic region at the peak in the smooth and rough model are shown in Figure 3. A velocity cut-off value of 1 m/s is shown in the contours. The near-wall velocity in the rough model is lower than in the smooth model, especially near the ridges. Meanwhile, the velocity near the ridge is lower than that near the valley. As can be seen, the near-wall velocity gradient in the smooth model is slightly larger than in the rough model. visualize the flow contour near the ridges and valleys, zoomed-in views of the slices across the stenotic region at the peak in the smooth and rough model are shown in Figure 3. A velocity cut-off value of 1 m/s is shown in the contours. The near-wall velocity in the rough model is lower than in the smooth model, especially near the ridges. Meanwhile, the velocity near the ridge is lower than that near the valley. As can be seen, the near-wall velocity gradient in the smooth model is slightly larger than in the rough model.  visualize the flow contour near the ridges and valleys, zoomed-in views of the slices across the stenotic region at the peak in the smooth and rough model are shown in Figure 3. A velocity cut-off value of 1 m/s is shown in the contours. The near-wall velocity in the rough model is lower than in the smooth model, especially near the ridges. Meanwhile, the velocity near the ridge is lower than that near the valley. As can be seen, the near-wall velocity gradient in the smooth model is slightly larger than in the rough model.

Pressure
The pressure contours of the smooth and rough model at the peak inlet velocity are presented in Figure 4. The pressure drop from the inlet to the outlet was about 8770 Pa in the smooth model and slightly smaller (8650 Pa) in the rough model. Over the cardiac cycle, the mean pressure dropped to 1500 Pa and 1480 Pa in the smooth and rough models, respectively. The overestimation rate of the smooth model was 1.33% compared to the rough model, calculated as the ratio (1500 − 1480)/1500. The difference between the mean pressure drops in the two models was not significant compared with other factors such as the degree of the stenosis [30]; this finding is consistent with the results published by Yakhot, Grinberg and Nikitin [13], as well as those of Johnston and Kilpatrick [15].

Pressure
The pressure contours of the smooth and rough model at the peak inlet velocity are presented in Figure 4. The pressure drop from the inlet to the outlet was about 8770 Pa in the smooth model and slightly smaller (8650 Pa) in the rough model. Over the cardiac cycle, the mean pressure dropped to 1500 Pa and 1480 Pa in the smooth and rough models, respectively. The overestimation rate of the smooth model was 1.33% compared to the rough model, calculated as the ratio (1500 − 1480)/1500. The difference between the mean pressure drops in the two models was not significant compared with other factors such as the degree of the stenosis [30]; this finding is consistent with the results published by Yakhot, Grinberg and Nikitin [13], as well as those of Johnston and Kilpatrick [15].

WSS
The TAWSS was calculated over the cardiac cycle and is presented in Figure 5. There is a minor difference in the upstream segments of the two models, and the high WSS occurs primarily in the middle of the stenotic region in both cases. The distributions of the TAWSS are generally similar, yet differences may be observed in the stenotic region due to the presence of the surface roughness.

WSS
The TAWSS was calculated over the cardiac cycle and is presented in Figure 5. There is a minor difference in the upstream segments of the two models, and the high WSS occurs primarily in the middle of the stenotic region in both cases. The distributions of the TAWSS are generally similar, yet differences may be observed in the stenotic region due to the presence of the surface roughness.  Figure 6 shows the oscillation of the WSS, and a cut-off value of 0.2 is present in the contours. A locally elevated OSI occurs primarily downstream of the stenosis, and its distributions are different in the models, as can be observed in the zoomed-in views. A high OSI can also be found at the ridges, highlighted by the arrows in the figure.   Figure 6 shows the oscillation of the WSS, and a cut-off value of 0.2 is present in the contours. A locally elevated OSI occurs primarily downstream of the stenosis, and its distributions are different in the models, as can be observed in the zoomed-in views. A high OSI can also be found at the ridges, highlighted by the arrows in the figure. The RRT incorporates the effects of TAWSS and OSI, i.e., low and oscillatory flow, which is an indicator of the vascular endothelial permeability (Figure 7). A high RRT is observed near the inlet in the models, mainly due to the low TAWSS at that position and the mathematical calculations based on the equation for RRT. The values of the RRT in the downstream segments in the models are larger than those in the stenotic segments, where the high-RRT region corresponds to the oscillatory flow in Figure 2. The RRT incorporates the effects of TAWSS and OSI, i.e., low and oscillatory flow, which is an indicator of the vascular endothelial permeability (Figure 7). A high RRT is observed near the inlet in the models, mainly due to the low TAWSS at that position and the mathematical calculations based on the equation for RRT. The values of the RRT in the downstream segments in the models are larger than those in the stenotic segments, where the high-RRT region corresponds to the oscillatory flow in Figure 2.
The TransWSS represents the multi-directionality of the WSS during the cardiac cycle. A high transWSS mainly occurs in the stenotic and downstream segments of the models (Figure 8). The zoomed-in views show that the distributions of transWSS are different, and a high transWSS can be found in the ridges in the rough model. The TransWSS represents the multi-directionality of the WSS during the cardiac cycle. A high transWSS mainly occurs in the stenotic and downstream segments of the models ( Figure 8). The zoomed-in views show that the distributions of transWSS are different, and a high transWSS can be found in the ridges in the rough model. The total surface area of the stenotic region was changed after adding the roughness, with values of 46.77 mm 2 and 48.11 mm 2 in the smooth and rough models, respectively. To compare the difference in WSS metrics more specifically, the ratio of the surface area under different thresholds to the total surface area in the stenotic region is presented in Table 2. The proportion of low TAWSS in the smooth model was less than that in the rough model (4.51% vs. 5.16%, threshold < 5 Pa). A threshold of > 0.2 was set to calculate the proportions of high OSI in the stenotic region, which were 1.97% and 2.47% in the smooth and rough models, respectively. Moreover, the proportions of high RRT (> 0.5 Pa −1 ) in the stenotic region of the smooth and rough models were 1.22% and 1.62%, yet the proportion of high transWSS in the stenotic region was slightly smaller in the rough model compared with the smooth model (4.95% vs. 4.44%, threshold > 10 Pa). From the statistics in the table, the maximum difference in WSS metrics induced by the roughness was 4.69%. The total surface area of the stenotic region was changed after adding the roughness, with values of 46.77 mm 2 and 48.11 mm 2 in the smooth and rough models, respectively. To compare the difference in WSS metrics more specifically, the ratio of the surface area under different thresholds to the total surface area in the stenotic region is presented in Table 2. The proportion of low TAWSS in the smooth model was less than that in the rough model (4.51% vs. 5.16%, threshold < 5 Pa). A threshold of > 0.2 was set to calculate the proportions of high OSI in the stenotic region, which were 1.97% and 2.47% in the smooth and rough models, respectively. Moreover, the proportions of high RRT (> 0.5 Pa −1 ) in the stenotic region of the smooth and rough models were 1.22% and 1.62%, yet the proportion of high transWSS in the stenotic region was slightly smaller in the rough model compared with the smooth model (4.95% vs. 4.44%, threshold > 10 Pa). From the statistics in the table, the maximum difference in WSS metrics induced by the roughness was 4.69%.   In summary, the near-wall velocity in the rough model was found to be lower, mainly due to the trapping effects of the flow inside the ridges and valleys. This also correlates with a larger RRT [13]. However, the pressure distributions of the two models did not vary significantly. Over the cardiac cycle, although the difference in the proportion under various WSS metrics in the stenotic region was slight, based on the calculations, the distributions of the WSS contours were different. A lower WSS was found in the rough model from the results for TAWSS and RRT. Moreover, a slightly higher OSI but lower transWSS may demonstrate a more oscillatory flow from the axial direction in the rough model. It is known that a low and oscillatory flow is associated with plaque formation and progression [31,32]. The results above indicate a slightly higher potential for atherosclerotic change in the rough model.
Bhaganagar [5] concluded that the region in the vicinity of rough plaque is directly affected by the plaque morphology, and the regions away from the rough wall are indistinguishable. In our study, the stenotic region was directly affected by the roughness of the plaque, and the differences in the upstream segment were minor and could be neglected. However, the velocity and WSS distributions in the downstream segment of the artery changed due to the rough plaque morphology, especially regarding the contours of OSI, RRT and transWSS. To specify the variations in WSS in the downstream region, the areas of high shear stress in the smooth and rough models were calculated under the same threshold, as mentioned previously (OSI > 0.2, 4.41 vs. 4.15 mm 2 ; RRT > 0.5 Pa −1 , 1.93 vs. 2.54 mm 2 ; transWSS > 10 Pa, 6.75 vs. 4.19 mm 2 ), and the differences were 5.90%, 31.60% and 37.93%, respectively. It can therefore be concluded that the stenotic roughness not only influences the flow in the stenotic segment but also changes the hemodynamics downstream of the stenosis.
Our study demonstrated that stenotic model simulations with a smooth surface may contribute to an inaccurate estimation of the WSS metrics. Apart from the smooth surface assumption, other assumptions or simplifications are often introduced, which may affect the outcomes of a hemodynamic study. The WSS variabilities resulting from the surface roughness in this study were as significant as the variations found from other factors such as the periodicity of the flow boundary waveform [33], blood viscosity [34] and inlet flow condition [35].
Currently, a number of papers focus on the study and evaluation of the functional severity of the coronary stenosis in terms of fractional flow reserve [36][37][38][39], but few analyze the plaque itself. The characteristics of the plaque can reveal its formation, development and even rupture. Herein, we analyzed the impact of the surface roughness on the hemodynamics, and the flow and WSS in the stenotic and downstream regions were changed if the effect was included.
One major limitation of this study is the assumption of a rigid wall. It remains unclear how the deformation of the arterial wall would influence the comparisons between smooth and rough surface simulations, and this may be addressed in future studies by considering both wall deformation and roughness. In addition, there is currently no standard for the selection of the threshold for the WSS metrics. The values of the threshold in our study were based on our simulation results and selected to present a better comparison. The roughness in the stenotic region was randomly generated, based on the stenosis morphology in this study. Realistic patient-specific roughness may be included in the future. Moreover, this study is only a preliminary analysis of the patient-specific roughness. More patient-specific models and experimental validations should be performed in the future to further elucidate the characteristics of the flow.

Conclusions
A patient-specific geometry with a smooth surface was constructed, and a random roughness was added only in the stenotic region to create a partially rough model. Oscillatory flow was demonstrated downstream in both models, and a lower near-wall velocity was observed in the rough model. The pressure distributions did not vary significantly depending on the presence of the surface roughness. In the stenotic region, the variations in the proportions of WSS metrics under different thresholds were slight, contributing to a slightly lower and oscillatory flow in the axial direction in the rough model. However, the differences in the WSS metrics in the downstream region became larger. The findings suggest that taking the stenotic surface roughness into account could be important for