Topological Optimization of Auxetic Coronary Stents Considering Hemodynamics

This paper is to design a new type of auxetic metamaterial-inspired structural architectures to innovate coronary stents under hemodynamics via a topological optimization method. The new architectures will low the occurrence of stent thrombosis (ST) and in-stent restenosis (ISR) associated with the mechanical factors and the adverse hemodynamics. A multiscale level-set approach with the numerical homogenization method and computational fluid dynamics is applied to implement auxetic microarchitectures and stenting structure. A homogenized effective modified fluid permeability (MFP) is proposed to efficiently connect design variables with motions of blood flow around the stent, and a Darcy-Stokes system is used to describe the coupling behavior of the stent structure and fluid. The optimization is formulated to include three objectives from different scales: MFP and auxetic property in the microscale and stenting stiffness in the macroscale. The design is numerically validated in the commercial software MATLAB and ANSYS, respectively. The simulation results show that the new design can not only supply desired auxetic behavior to benefit the deliverability and reduce incidence of the mechanical failure but also improve wall shear stress distribution to low the induced adverse hemodynamic changes. Hence, the proposed stenting architectures can help improve safety in stent implantation, to facilitate design of new generation of stents.

Introduction

Percutaneous coronary intervention (PCI) has been popular as a common treatment for coronary artery disease, but the risks of stent thrombosis (ST) and in-stent restenosis (ISR) still threaten the safety of stent implantation and represent a serious clinical shortfall. Although various reasons accounting for these adverse biological responses that have not been fully understood, the stent essentially serves as a mechanical structure and is believed to have a major effect on the ST and ISR complications. In our previous work (Xue et al., 2020), the self-expanding (SE) auxetic stent has been demonstrated to have capability in supplying adaptive deformation to help overcome the mechanical failures and low the incidence of the ST and ISR complications. In this research, the stent-induced hemodynamic changes will also be combined into stents to further improve stenting performance from both mechanical and hemodynamic aspects, because the stenting structures are primally influenced not only by mechanical failure including inadequate expansion, fracture, malapposition, foreshortening and dogbone, but also by the hemodynamic changes induced in stent implantation.

The adverse hemodynamic changes caused by implanting stents are also related to high risks of ISR and ST (Tomaszewski et al., 2020; Gharleghi et al., 2021). In 1995, Rogers and Edelman (1995) noted that the stenting structures might affect restenosis and thrombosis in the arteries through actions other than the mechanical injures. Then, neointimal hyperplasia independent of arterial injury was first investigated in stented rabbit iliac arteries (Garasic et al., 2000). The results indicated that stenting structures might adversely affect fluid motions in the arterial lumen and induce intimal thicken. Nearly in the same period, Berry et al. (2000) shown that enlarging strut spacings can facilitate blood flow and reduce the cumulation of neointima in a stented coronary artery. Then, the inverse proportion relationship between neointimal hyperplasia and the change of wall shear stress (WSS) after stent implantation was found by Wentzel et al. (2001) and further studied by LaDisa et al. (2005).

In hemodynamics, the WSS on the endothelial cells of coronary arteries has been demonstrated as a significant factor connected with atherosclerotic disease development (Ku, 1997). The protrusions of stent struts inside the arterial lumen disturb blood flow, so as to change the local flow environment, and further affect the distributions of wall stress in the stented segment. The shear stress is generated when blood flows through the endothelium, which is proportional to the blood viscosity and the velocity gradient. The WSS is related to the vascular hemodynamics and correlated with radial responses, intimal thickening, and platelet thrombosis (Gijsen et al., 2019). The WSS under a steady laminar flow can promote endothelial cells to release factors to inhibit coagulation, migration of leukocytes, and smooth muscle proliferation, to benefit the healing of narrowed vessels (Traub and Berk, 1998). However, unusual WSS can cause neointimal hyperplasia and atherosclerotic plaque formation. Many studies have shown that the areas with WSS lower than 0.5 Pa are prone to intimal thickening and atherosclerosis, while the WSS higher than 2.5 Pa may increase the risk of plaque rupture and thrombosis (Beier et al., 2016; Gharleghi et al., 2021). The distributions of WSS will change after implanting stents due to the altered velocity profile of the blood flow in the stented segment (Gijsen et al., 2019). Therefore, WSS has been widely adopted as a metric to evaluate the impact of stent implantation on hemodynamic changes. Sometimes, for nonuniform flow and pulsatile flow, wall shear stress gradient (WSSG) and time-averaged wall shear stress (TAWSS) are also used to evaluate the stent, respectively (Gharleghi et al., 2019). Although the stent induced hemodynamic changes can be quantified via the WSS distribution, how the stenting structures affect the WSS in the stented segment is still unclear.

We can find that the stent protrusion is related to the WSS and affects the adverse hemodynamic. The influence of the stent protrusion can be reduced but it is inevitable in stent implantation. For instance, the flow disturbance can be eliminated by fully embedding stent architectures into artery walls (Gharleghi et al., 2019). However, by doing this, the over-sized stent can lead to high tension stress on the vessel wall and may cause serious injury (Nolan and Lally, 2018). As a result, the influences of strut geometric characteristics on flow patterns have been widely investigated (Putra et al., 2019; Wei et al., 2019; Gharleghi et al., 2021) to minimize the adverse hemodynamic effect on stents via designs. Among stenting geometric characteristics, strut thickness contributes the most concern because it directly determines the degree of stent protrusion and is strongly connected to hemodynamic changes. Thicker struts were found to significantly promote intimal thickening (Stiehm et al., 2019). Strut spacing is another factor, as a large spacing can restore disturbed flow (Bekal et al., 2018). Since the adverse hemodynamic changes are usually presented as flow disturbance, streamlined stent strut cross-section profiles were suggested to facilitate blood flow (Melzer et al., 2019). Also, varying strut angles (Beier et al., 2016) were also found to benefit hemodynamic designs, such as aligning the orientations of struts and connectors with the flow direction (Bekal et al., 2018).

Most current stent designs usually utilize empirical or surrogate model-based methods to study the influence of variables, such as struts thickness, width, angle, and spacings on WSS. The earlier works were mainly focused on a single objective, such as minimizing low WSS areas (Gundert et al., 2012). However, the stent usually works in a complex environment associated with structural mechanics and fluid dynamics (Bressloff et al., 2016). Thus, several studies have tried to perform design by accounting for both aspects, simultaneously. For example in the work (Putra et al., 2019) the design objectives were defined based on both mechanical and hemodynamic metrics, including recoil, flexibility, and WSS. Objectives related to drug release were also taken into account in stent design (Vijayaratnam et al., 2019). However, these optimization formulations mostly adopt a single objective to consider hemodynamic and a small number of structural design variables related to sizes, which in practice are hard to capture the impacts of stenting structures on the blood flow alterations in the arteries. Hence, two objectives, such as mean square WSS and mean swirl value (Doutel et al., 2018), or recirculation zone length and struts reattachment distance (Gharleghi et al., 2019) have been defined to implement the design, although the improvements were limited. After that, researchers attempted to include more parameters with strut geometric characteristics into the optimization (Gharleghi et al., 2021), while too many variables and multiple design objectives lead to very complex optimization formulations.

The above designs narrow the freedom to generate new stents by only allowing the change of stent geometrical parameters. Particularly, the conventional designs are not able to change structural topologies and shapes for uncovering novel architectures to improve stenting performance. Hence, new design methodologies are in demand. Recently, topology optimization is experiencing popularity as an efficient tool for systematic design of new structures. Since stents subject to hemodynamic changes are sensitive to both the shape and topology of stenting structures, this paper employs a parametric level set method (PLSM) (Luo et al., 2007; Luo et al., 2009) to implement topological design of stents. PLSM can generate structures with clear and smooth boundaries, more suitable for hemodynamic stent designs. PLSM has been successfully applied designs of mechanical metamaterials and auxetic structures (Wang et al., 2014; Wu et al., 2017; Li et al., 2018). It is noted that a stent usually consists of a number of uniform and repetitive structural architectures. In this paper the numerical homogenization method is used to evaluate effective properties (e.g., the elasticity tensor) of underline microstructures. This paper integrates the numerical homogenization method and Computational fluid dynamics (Gundert et al., 2012) with the PLSM approach (Luo et al., 2007) to fulfill topological design of stents using auxetic metamaterials. Auxetics are engineered mechanical metamaterials consisting of periodic microstructures with effective negative Poisson’s ratio (NPR) (Lakes, 1993), exhibiting counterintuitive property in deforming structures: contract or expand in a transverse direction when compressed or stretched in an axial direction. Such behaviour in elastic deformation will considerably benefit stents, e.g., energy absorption, adaptive stiffness, indentation resistance, and enhanced fracture toughness, which gives the stenting structures enhanced flexibility and adaptability to automatically match the artery shape during the deformation.

A modified fluid permeability (MFP) is proposed to quantify obstructions to fluid at different directions in stenting microstructures, to include the hemodynamic changes into stent design to describe the stent-induced fluid changes. A fluid permeability is often used to quantify the ability of a porous medium to allow fluids passing through it. It is mainly determined by the porosity, shapes of pores, and their distributions. The porous medium with higher permeability can enable the fluid to easily move through it. Implanting stent changes flow environments due to the obstructions from stent struts (Putra et al., 2019). The obstructions can cause flow disturbances, flow separation, recirculation zones, resulting in adverse clinical outcomes (Bukač et al., 2019). Then, reducing the adverse hemodynamic changes around stents can be obtained by minimizing the MFP in directions of blood flow.

Assumptions in Stent Designs

To perform a hemodynamic topology optimization for the stent, the simulation of blood flow in the stented artery is necessary. The complex endovascular flow environment and the interactions between blood, stent, and endothelial cells are challenging in practice. Although CFD provides an efficient way to describe fluid motion by solving Naiver-Stokes (NS) equations, a complex blood flow model still leads to a high computational cost in solving the full NS equations. Some assumptions necessary to simply the fluid computation and facilitate the design optimization are given as follows.

Blood Flow

Blood is an incompressible fluid, so the NS equations with incompressible conditions are utilized to describe the motion of blood flow. Eq. 1 is the momentum equation, Eq. 2 is the incompressibility condition, and Eq. 3 is the no-slip condition on the boundary.

$\rho \frac{\partial \mathbf{u}}{\partial t}+\rho \left(\mathbf{u}\cdot \nabla \right)\mathbf{u}-\mu \Delta \mathbf{u}=-\nabla p+\rho \mathbf{g}\text{in}{\text{Ω}}_{\text{f}}(1)$

$\nabla \cdot \mathbf{u}=0\text{in}{\text{Ω}}_{\text{f}}(2)$

$\mathbf{u}=0\text{on }{\text{Γ}}_{\text{f}}(3)$

where Ω_f is a fluid domain with boundary Γ_f. ρ and μ denote fluid density and viscosity, respectively. u is the flow velocity, t is time, p is the pressure, and g is body accelerations acting on the fluid, such as gravity and inertial. ρg denotes the external body force vector including inertia forces. ∇, Δ and ∇⋅ are the gradient, Laplacian, and divergence operators, respectively. The assumptions are:

1) Geometric assumptions: the stent is assumed as a straight, rigid body with constant thickness and uniformed horizontal cross-section profile, and no stent deformations are generated under the blood flow. These assumptions will reasonably simplify the CFD model (Balossino et al., 2008).

2) Blood flow is assumed as an incompressible Newtonian fluid with a steady flow state, so the unsteady item in Eq. 1 is neglected. Based on the study of (Johnston et al., 2006), only about 30% cardiac cycle was obviously affected by the non-Newtonian model. Blood flows in all arteries present pulsatile motions due to the cyclic behaviors of systole and diastole (Ku, 1997). Hence, the blood flow with a steady-state (Mattace-Raso et al., 2006) is usually adopted, and the inlet velocity is defined as a time-averaged value over a cardiac cycle.

3) The optimization intends to improve the blood flow in this worst region near the artery walls. Since the blood flow is relatively slow near the artery walls, the convective and inertia items can be ignored. Hence, the blood flow is consequently simplified to a Stokes flow in the computational domain, where the momentum equations are given in Eq. 4.

$\begin{array}{l}-\mu \Delta \mathbf{u}=-\nabla p+\rho \mathbf{g}\text{in}{\text{Ω}}_{\text{f}}\\ \nabla \cdot \mathbf{u}=0\text{in}{\text{Ω}}_{\text{f}}\\ \mathbf{u}=0\text{on}{\text{Γ}}_{\text{f}}\end{array}(4)$

Design Domain

In the stented segment of the artery, blood flows through the inner lumen and fills the gaps between stent struts. The fluid is separated by the stent and divided into different regions. A 3D computational domain based on one stenting microstructure is established for the CFD model, as shown in Figure 1. The domain has four regions: three blood regions with red colour and one stent region with blue colour. In the middle of the computational domain, two thin layers with the same thicknesses are defined as the stenting microstructures, with one microstructure defined in the stent region at the bottom layer. Blood flows through the whole top layer and fills the gaps between stent struts in the bottom layer. Hence, the stent region consists of one stent microstructure filled with the blood in gaps. Two blood regions with the same sizes are defined in the two ends of the microstructure to avoid the impacts of the flow boundaries. Since the fluid convection is ignored, various lengths of inlet and outlet will affect the fluid. The CFD simulation is then performed under the following conditions:

1) A constant inlet velocity of blood flow specified on the left side of the computational domain.

2) The outlet boundary with nearly zero pressure located on the right side.

3) No-slip conditions applied on the outside surfaces of the simulation region.

4) No additional conditions to the solid-fluid interfaces due to the restriction to the flow in solid areas.

FIGURE 1

FIGURE 1. The 3D computational domain.

As the assumption of constant thickness, a 2D horizontal cross-section profile is extracted from a stent microstructure and defined as the design domain. The outer contours of the design domain are presented with white lines, and the bottom of the stent region is defined as the artery wall, as given in Figure 2. The velocities at the design domain are extracted to evaluate the index MFP. After that, the design variables are updated by the optimization method. Then, the material layout in the stent region can be correspondingly updated by the new microstructural profile. In this way, the associations between the computation and design domains are established.

FIGURE 2

FIGURE 2. The 2D design domain.

Parametric Level Set Method

LSM implicitly defines the structural boundary as the zero-level set of a higher-dimensional scalar function Φ (x). For example, a two-dimensional (2D) design, given in Figure 3, is defined by

$\{\begin{array}{l}\mathit{\Phi }\left(x\right)>0x\in \mathit{\Omega }\backslash \partial \mathit{\Omega }\left(\text{Solid}\right)\\ \mathit{\Phi }\left(x\right)=0x\in \partial \mathit{\Omega }\left(\text{Boundary}\right)\\ \mathit{\Phi }\left(x\right)<0x\in D\backslash \left(\mathit{\Omega }\cup \partial \mathit{\Omega }\right)\left(\text{Void}\right)\end{array}(5)$

where x is a point in space D. Ω and ∂Ω are the design domain and its boundary, respectively. The dynamic motion of the level set function then drives the topological shape changes of the structure. Differentiating zero level set function with respect to a pseudo time t leads to motion of interface as

$\frac{\partial \mathit{\Phi }\left(x,t\right)}{\partial t}-v_n\left|\nabla \mathit{\Phi }\left(x,t\right)\right|=0(6)$

FIGURE 3

FIGURE 3. (A) 3D level set surface; 2D level set boundary.

Thus, the level set function can be solved by finding a velocity field v_n to maintain the first-order Hamilton-Jacobi partial differential equation (H-J PDE) in Eq. 6. In the PLSM (Luo et al., 2007), the level set function is interpolated via centrally positioned CSRBFs at a set of given fixed knots N over the whole design domain by

$\mathit{\Phi }\left(x,t\right)=\phi {\left(x\right)}^T\alpha \left(t\right)={\displaystyle \sum _{i=1}^N\phi \left(x\right){\alpha }_i\left(t\right)}(7)$

The vector with the CSRBFs functions is defined as:

$\phi \left(x\right)={\left[{\varphi }_1\left(x\right),{\varphi }_2\left(x\right),\mathrm{...},{\varphi }_N\left(x\right)\right]}^T(8)$

The expansion coefficient vector is given by:

$\alpha \left(t\right)={\left[{\alpha }_1\left(t\right),{\alpha }_2\left(t\right),\mathrm{...},{\alpha }_N\left(t\right)\right]}^T(9)$

In this paper, the CSRBFs with C2 continuity is used (Luo et al., 2007; Luo et al., 2009). Because of CSRBF knots fixed in space, the shape functions φ (x) are spatial only and the expansion coefficients α (t) are time-dependent only, leading to a separation of the time and space in the level set function. Thus, the H-J PDE is transformed as the following ordinary differential equation system:

$\phi {\left(X\right)}^T\dot{\alpha }\left(t\right)-v_n\left|{\left(\nabla \phi \right)}^T\alpha \left(t\right)\right|=0(10)$

In this way, the normal velocity field is only determined by the parameters α (t):

$v_n=\frac{\phi {\left(X\right)}^T}{\left|{\left(\nabla \phi \right)}^T\alpha \left(t\right)\right|}\dot{\alpha }\left(t\right),\text{ where }\dot{\alpha }\left(t\right)=\frac{d\alpha \left(t\right)}{dt}(11)$

Darcy-Stokes Coupling

The application of the Darcy-Stokes coupling approach in saturated porous mediums (Guest and Prévost, 2006) provides an efficient way to deal with the moving-boundary no-slip conditions for solid-fluid interfaces. In the system, the solid region is treated as a porous medium with flow governed by Darcy’s law, and the fluid area is treated as Stokes flow. A unified analysis model that combines the Darcy and Stokes equations can then be established. In the model, a penalization to flow in the solid region is taken by assigning a low permeability to drive the velocity to close zero. Thus, the no-slip conditions along the solid-fluid interface are maintained (Guest and Prévost, 2006; Wang et al., 2016). By using the Darcy-Stokes coupling method, no additional boundary conditions are needed for solid-fluid interfaces, which saves the computational time and benefits the convergence of the optimization.

Stokes flow is usually used to describe a steady viscous fluid with slow velocity by ignoring convective and inertia items in the momentum equations. A typical Stokes flow is formulated in Eq. 4. Darcy flow, derived from the NS equations, usually describes a fluid through a porous medium by the homogenization method (Whitaker, 1986). In this work, the solid regions in the design domain can be treated as porous mediums govern by Darcy flow, where the nodal velocities are close to zero. The typical equations of Darcy flow are given as:

$\begin{array}{l}\mathbf{u}=-\frac{k}{\mu }\left(\nabla p-\mathbf{f}\right)\text{in}{\text{Ω}}_{\text{f}}\\ \nabla \cdot \mathbf{u}=0\text{in}{\text{Ω}}_{\text{f}}\\ \mathbf{u}=0\text{on}{\text{Γ}}_{\text{f}}\end{array}(12)$

where k is the fluid permeability, and f is the vector of external force. When using the Darcy-Stokes method, the velocity and pressure are dependent variables but with different orders in both flows. The incompressible condition has a stability requirement for the combination of the interpolation functions. Here, the stabilized mixed finite element methods (Hughes et al., 1986) is used for Stokes flow, and the method proposed by Masud and Hughes (2002) is applied for Darcy flow. Both stabilized mixed methods can avoid restrictions of the Babuška-Brezzi conditions (Babuška, 1971).

The stabilized matrix form for Darcy-Stokes coupling is given by

$\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{ds}}& -{\mathbf{G}}_{\mathrm{d}\mathrm{s}}\\ {\mathbf{G}}_{\text{ds}}^{\text{T}}& {\mathbf{M}}_{\mathrm{d}\mathrm{s}}\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_{\mathrm{d}\mathrm{s}}\\ {\mathbf{H}}_{\mathrm{d}\mathrm{s}}\end{array}\right](13)$

Combined within the PLSM, the matrices in the Darcy-Stokes system are defined as:

$\begin{array}{l}{\mathbf{K}}_{\text{ds}}\left(x\right)=A\left\{H\left(\text{Φ}\left(x\right)\right){\mathbf{K}}_{\text{d}}+\left[1-H\left(\text{Φ}\left(x\right)\right)\right]{\mathbf{K}}_{\text{s}}\right\}\\ {\mathbf{G}}_{\text{ds}}\left(x\right)=A\left\{H\left(\text{Φ}\left(x\right)\right){\mathbf{G}}_{\text{d}}+\left[1-H\left(\text{Φ}\left(x\right)\right)\right]{\mathbf{G}}_{\text{s}}\right\}\\ {\mathbf{G}}_{\text{ds}}^{\text{T}}\left(x\right)=A\left\{H\left(\text{Φ}\left(x\right)\right){\mathbf{G}}_{\text{d}}^{\text{T}}+\left[1-H\left(\text{Φ}\left(x\right)\right)\right]\left({\mathbf{L}}_{\text{s}}+{\mathbf{G}}_{\text{s}}^{\text{T}}\right)\right\}\\ {\mathbf{M}}_{\text{ds}}\left(x\right)=A\left\{H\left(\text{Φ}\left(x\right)\right){\mathbf{M}}_{\text{d}}+\left[1-H\left(\text{Φ}\left(x\right)\right)\right]{\mathbf{M}}_{\text{s}}\right\}\end{array}(14)$

where A is the standard finite element routine. The subscript “ds” denotes the matrices defined in the Darcy-Stokes system. The subscripts “d” and “s” emphasize Darcy and Stokes flow, respectively. u is the velocity vector, and p is the pressure vector. K is the viscosity stiffness matrix, G is the gradient matrix, G^T is the divergence matrix, L is the consistency matrix, and M is the stabilization matrix. F and H are the nodal forces. H [Φ(x)] is the Heaviside function (Wang et al., 2003; Luo et al., 2007) of Φ(x) at point x. The interpretations indicate solid elements have Darcy stiffness, and void elements have Stokes stiffness.

Homogenization of the Stent Properties

The stenting structure consists of periodic array of identical microstructures. A near wall horizontal cross-section profile of the microstructure is specified as the 2D design domain. The numerical homogenization method is adopted to compute the effective elasticity tensor and the MFP tensor.

Homogenization of Elasticity

The effective elasticity tensor D^H_ijkl of the microstructure is assessed by the homogenization method as:

$D_{ijkl}^{\text{H}}=\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\ast \left(ij\right)}\right)D_{pqrs}\left({\epsilon }_{pq}^{0\left(kl\right)}-{\epsilon }_{pq}^{\ast \left(kl\right)}\right)dY}(15)$

where Y is the micro design domain and denotes one unit cell of the stent, and |Y| is the area of the unit cell. i, j, k, l = 1, 2. D_pqrs is the elasticity tensor of the solid material in the design domain. ${\mathit{\epsilon }}_{\mathit{pq}}^{\mathit{0}\mathit{(ij)}}$ is the test unit strain field, where (1,0,0)^T, (0,1,0)^T and (0,0,1)^T are used in 2D models. ${\mathit{\epsilon }}_{\mathit{pq}}^{\mathit{*}\mathit{(ij)}}$ is the locally varying strain fields and defined by:

${\epsilon }_{pq}^{\ast \left(ij\right)}=\frac{1}{2}\left(\frac{\partial u_p^{ij}}{\partial y_q}+\frac{\partial u_q^{ij}}{\partial y_p}\right)(16)$

The displacement field u^(ij) is calculated by applying the periodical boundary conditions as

${\displaystyle {\int }_Y\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\ast \left(ij\right)}\left(u^{\left(ij\right)}\right)\right)D_{pqrs}{\epsilon }_{rs}^{\ast \left(kl\right)}\left(v^{\left(kl\right)}\right)dY}=0,\forall v^{\left(kl\right)}\in U_p\left(Y\right)(17)$

where ν^(kl) is a virtual displacement field in U_p(Y). Note that U is the kinematically admissible displacement space comprised of periodic Y.

The effective stiffness matrix of the material that has four independent components is defined by

$D^{\text{H}}=\left[\begin{array}{ccc}{D}_{11}^{\text{H}}& D_{12}^{\text{H}}& 0\\ D_{12}^{\text{H}}& D_{22}^{\text{H}}& 0\\ 0& 0& D_{33}^{\text{H}}\end{array}\right](18)$

The auxetic properties can be obtained through the design having effective negative Poisson’s ratios (NPR). Here, two Poisson’s ratios Mu1 and Mu2 in two directions are used to define the related design objectives, and the calculations of the two Poisson’s ratios are given by:

$\{\begin{array}{c}Mu1=D_{12}^{\text{H}}/D_{11}^{\text{H}}\\ Mu2=D_{12}^{\text{H}}/D_{22}^{\text{H}}\end{array}(19)$

Homogenization of Modified Fluid Permeability

A permeability in fluid mechanics is usually used to quantify the ability of a porous material for fluids to go through it. The permeability can also be applied to measure the obstructions caused by stents. Considering the minimization of the stenting hemodynamic alterations, the optimization aims to reduce the obstruction effect of the stent on the blood flow. Since the blood flow is above the stent layer, rather than across through the stent region, the MFP is different from the typical fluid permeability. By applying the homogenization method in Stokes equations, Darcy’s law can be obtained:

$\mathbf{U}=-\frac{1}{\mu }{\mathbf{K}}^{\text{H}}\left(\nabla p-\mathbf{f}\right)(20)$

where U is the vector of average velocities. μ and f are the viscosity and the external body force, respectively. ∇p is the pressure gradient. K^H is the homogenized effective fluid permeability tensor. The calculation of a typical permeability can be formulated as:

${\mathbf{K}}^{\text{H}}=\left[K_{\text{ij}}^{\text{H}}\right]=\frac{1}{\left|Y\right|}{\left({\mathrm{w}}^{\left(i\right)}\right)}^T{\mathbf{K}}_{\text{ds}}{\mathrm{w}}^{\left(j\right)}(21)$

where the effective fluid permeability tensor K^H is assembled by the components K^H_ij in the principal direction. w is the velocity vector solved from the above Darcy-Stokes coupling system via Eq. 13. K_ds is the Darcy-Stokes viscosity matrix defined in Eq. 14. After that, the MFP tensor is proposed for the 2D design domain. The formulation of the MFP derived from Eq. 21 can be defined as:

${\mathbf{K}}_{2\text{D}}^{\text{H}}=\left[K_{2\text{D}\left(ij\right)}^{\text{H}}\right]=\frac{1}{\left|Y_{2D}\right|}{\left({\mathbf{w}}_{2D}^{\left(i\right)}\right)}^T{\mathbf{K}}_{\text{ds}}^{2\text{D}}{\mathbf{w}}_{2D}^{\left(j\right)}(22)$

where the subscript “2D” denotes the extracted 2D design domain. K^H_2D is the MFP tensor. w_2D is the velocity vector in the 2D design domain, which is calculated in the 3D computational domain. K^2D_ds is the Darcy-Stokes viscosity matrix defined for the 2D design domain and derived from the mixed viscosity matrix K_ds. The velocity vector w_2D contains two directional components: one is in the blood flow direction, and another one is normal to the blood flow. Thus, the purpose of reducing the stenting obstructions can be obtained by minimizing the vertical component of the MFP.

Optimization Formulation and Sensitivity Analysis

This section will formulate the multiscale multi-objective topology optimization problem, and the sensitivity analyse of the optimization objectives will also be derived. Since different objectives will be associated with different design variables. The mechanical objectives are auxetic properties and stent stiffness, while the MFP is related to CFD simulation. After that, the connection between the MFP and design variables can be built, and the sensitivity analysis can be performed.

Optimization Model

The numerical optimization starts from an initial guess in the 2D design domain, where the expansion coefficients of the CS-RBF interpolation in the PLSM are defined as design variables. Then, the three objectives are computed simultaneously. With respect to the mechanical objectives, the compliance of the macroscopic stenting structure can be computed by the obtained effective elasticity tensor. Regarding the fluidic objective, the 3D computational domain is first established based on the initial 2D design domain. The material distribution of the 3D microstructure is determined by the 2D profile in the design domain. The surrounding blood regions are specified the same as the initial definition in Figure 1. The coupled Darcy-Stokes system can be solved, and the velocity vector of the 3D fluid model is then achieved. After that, the velocity vector of the 2D design domain is extracted from the 3D fluid model to calculate the MFP. Thus, all three objectives and related sensitivities are obtained. They are then normalized and weighted as an equivalent objective function. Based on the sensitivity information, the PLSM is utilized to update the structural topological shapes.

The proposed multiscale topology optimization model is formulated as:

$\{\begin{array}{l}Find{\alpha }_{2D,n}^{MI}\left(n=1,2,\mathrm{...},N\right)\\ MinJ\left({\alpha }_{2D}^{MI}\right)=W_1{J}_1\left({\alpha }_{2D}^{MI}\right)+W_2{J}_2\left({\alpha }_{2D}^{MI}\right)+W_3{J}_3\left({\alpha }_{2D}^{MI}\right)\\ =W_1{J}_A^{MI}\left({\alpha }_{2D}^{MI}\right)+W_2{J}_P^{MI}\left({\alpha }_{2D}^{MI}\right)+W_3{J}^{MA}\left({\alpha }_{2D}^{MI}\right)\\ S.t.V\left({\alpha }_{2D}^{MI}\right)={\displaystyle {\int }_{{\Omega }^{MI}}H\left({\mathit{\Phi }}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)d{\mathit{\Omega }}_{2D}^{MI}-V_{2D}^{max}\le 0}\\ F_{2D}^{MI}\left(u_{2D}^{MI},w_{2D}^{MI},{\alpha }_{2D}^{MI}\right)=L_{2D}^{MI}\left(w_{2D}^{MI},{\alpha }_{2D}^{MI}\right),\forall w_{2D}^{MI}\in \overline{U}\left({\mathit{\Omega }}_{2D}^{MI}\right)\\ F_{3D}^{MI}\left(u_{3D}^{MI},p_{3D}^{MI},v_{3D}^{MI},q_{3D}^{MI},{\alpha }_{2D}^{MI}\right)=L_{3D}^{MI}\left(v_{3D}^{MI},q_{3D}^{MI}\right),\forall v_{3D}^{MI},q_{3D}^{MI}\in \overline{U}\left({\mathit{\Omega }}_{3D}^{MI}\right)\\ F_{2D}^{MA}\left(u_{2D}^{MA},w_{2D}^{MA},D_{ijkl}^H\right)=L_{2D}^{MA}\left(w_{2D}^{MA}\right),\forall w_{2D}^{MA}\in \overline{U}\left({\mathit{\Omega }}_{2D}^{MA}\right)\\ {\alpha }_{2D,min}^{MI}\le {\alpha }_{2D,n}^{MI}\le {\alpha }_{2D,max}^{MI}\end{array}(23)$

where $J_A^{MI}={\left(Mu1\left({\alpha }_{2D}^{MI}\right)+1\right)}^2+{\left(Mu2\left({\alpha }_{2D}^{MI}\right)+1\right)}^2$, $J_P^{MI}={\mathbf{K}}_{2\text{D}}^{\text{H}\left(\mathrm{2,2}\right)}\left({\alpha }_{2D}^{MI}\right)$ and $J^{MA}=\frac{1}{2}{\displaystyle {\int }_{{\Omega }_{2D}^{MA}}{\epsilon }_{ij}^T\left(u_{2D}^{MA}\right)D_{ijkl}^H\left({\alpha }_{2D}^{MI}\right){\epsilon }_{kl}\left(u_{2D}^{MA}\right)d{\mathit{\Omega }}_{2D}^{MA}}$.

J is the objective function, consisting of the auxetic property $J_A^{MI}$, the fluid objective $J_p^{MI}$, and macro compliance J^MA, where W₁, W₂, W₃ are corresponding weight factors. The superscript “MA” and “MI” denote parameters in macroscale and microscale, respectively; the subscript “2D” refers to parameters in the 2D design domain, and “3D” means parameters in 3D domain. N is the total number of fixed knots in the 2D design domain. The coefficients of the interpolation ${\alpha }_{2D,n}^{MI}$ are the design variables, ranging from ${\alpha }_{2D,\mathit{min}}^{MI}$ to ${\alpha }_{2D,\mathit{max}}^{MI}$. The subscript “A” and “P” denote auxetic and permeability, respectively. K^H_2D is the MFP, and K^H(2,2)_2D is the vertical component. V is the volume constraint with the upper limit $V_{\mathrm{2D}}^{Max}$, applied to the 2D microstructure but can also restrict the 3D microstructure.

The MFP related to the design variables can be calculated by:

${\mathbf{K}}_{2\text{D}}^{\text{H}}\left({\alpha }_{2D}^{MI}\right)=\left[{\mathbf{K}}_{2D\left(ij\right)}^{\text{H}}\left({\alpha }_{2D}^{MI}\right)\right]=\frac{1}{\left|{\mathit{\Omega }}_{2D}^{MI}\right|}{\left({\mathbf{w}}_{2D}^{\left(i\right)}\right)}^T{\mathbf{K}}_{\text{ds}}^{2\text{D}}\left(H^f\left({\mathit{\Phi }}_{2D}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)\right){\mathbf{w}}_{2D}^{\left(j\right)}(24)$

where the velocity vector w_2D is extracted from the velocity vector ${\mathbf{u}}_{\mathrm{3D}}^{MI}{\mathbf{u}}_{\mathrm{3D}}^{MI}$, which is calculated by the following momentum equation performed in the 3D computational domain.

$F_{3D}^{MI}\left(u_{3D}^{MI},p_{3D}^{MI},v_{3D}^{MI},q_{3D}^{MI},{\alpha }_{2D}^{MI}\right)=L_{3D}^{MI}\left(v_{3D}^{MI},q_{3D}^{MI}\right)(25)$

where $v_{3D}^{MI}$ and $P_{3D}^{MI}$ are the boundary conditions defined in the CFD model, given by

$\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{d}\mathrm{s}}\left(H^f\right)& -{\mathbf{G}}_{\mathrm{d}\mathrm{s}}\left(H^f\right)\\ {\mathbf{G}}_{\mathrm{d}\mathrm{s}}^{\mathrm{T}}\left(H^f\right)& {\mathbf{M}}_{\mathrm{d}\mathrm{s}}\left(H^f\right)\end{array}\right]\left[\begin{array}{c}\mathbf{u}\left(u_{3D}^{MI}\right)\\ \mathbf{p}\left(p_{3D}^{MI}\right)\end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_{\mathrm{d}\mathrm{s}}\left(w_{3D}^{MI}\right)\\ {\mathbf{H}}_{\mathrm{d}\mathrm{s}}\left(q_{3D}^{MI}\right)\end{array}\right](26)$

where H^f denotes the Heaviside function H applied to the Darcy-Stokes coupling in the CFD model. Θ in H^f is equal to zero to identify fluid material in the design domain. The $D_{ijkl}^H$ used for the evaluation of the mechanical properties in the optimization model can be calculated by:

$D_{ijkl}^{\text{H}}\left({\alpha }_{2D}^{MI}\right)=\frac{1}{\left|{\mathit{\Omega }}_{2D}^{MI}\right|}{\displaystyle {\int }_{{\mathit{\Omega }}_{2D}^{MI}}\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\ast }\left(u_{2D}^{\left(ij\right)}\right)\right)D_{pqrs}\left({\epsilon }_{rs}^{0\left(kl\right)}-{\epsilon }_{rs}^{\ast }\left(u_{2D}^{\left(kl\right)}\right)\right)H\left({\mathit{\Phi }}_{2D}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)d{\mathit{\Omega }}_{2D}^{MI}}(27)$

where D_pqrs is the elasticity tensor of the solid. ${\mathit{\epsilon }}_{\mathit{pq}}^{\mathit{0}}$ is the test unit strain field, where (1,0,0)^T, (0,1,0)^T and (0,0,1)^T are used. ${\mathit{\epsilon }}_{\mathit{pq}}^{\mathit{*}}$ is the strain field related to the displacement u_2D that can be calculated by

${\displaystyle {\int }_{{\Omega }_{2D}^{MI}}\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\ast }\left(u_{2D}^{MI\left(ij\right)}\right)\right)D_{pqrs}{\epsilon }_{rs}^{\ast }\left(w_{2D}^{MI\left(kl\right)}\right)H\left({\Phi }_{2D}^{MI}\right)d{\Omega }_{2D}^{MI}}=0,\forall w\in \overline{U}\left({\Omega }_{2D}^{MI}\right)(28)$

where w is the virtual displacement field.

The bilinear energy and the linear load forms in the 2D microscale domain are given by

$F_{2D}^{MI}\left(u_{2D}^{MI},w_{2D}^{MI},{\alpha }_{2D}^{MI}\right)={\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\epsilon }_{ij}^{\ast }\left(u_{2D}\right)D_{pqrs}{\epsilon }_{kl}^{\ast }\left(w_{2D}^{MI}\right)H\left({\Phi }_{2D}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)d{\Omega }_{2D}^{MI}}(29)$

$L_{2D}^{MI}\left(w_{2D}^{MI},{\alpha }_{2D}^{MI}\right)={\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\epsilon }_{ij}^{0\left(ij\right)}\left(u_{2D}\right)D_{pqrs}{\epsilon }_{kl}^{\ast }\left(w_{2D}^{MI}\right)H\left({\Phi }_{2D}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)d{\Omega }_{2D}^{MI}}(30)$

The bilinear energy and the linear load forms in the 2D macroscale domain can be computed by:

$F_{2D}^{MA}\left(u_{2D}^{MA},w_{2D}^{MA},D_{ijkl}^H\right)={\displaystyle {\int }_{{\Omega }_{2D}^{MA}}{\epsilon }_{ij}\left(u_{2D}^{MA}\right)D_{ijkl}^H{\epsilon }_{kl}\left(w_{2D}^{MA}\right)d{\Omega }_{2D}^{MA}}(31)$

$L_{2D}^{MA}\left(w_{2D}^{MA}\right)={\displaystyle {\int }_{{\Omega }_{2D}^{MA}}p{w}_{2D}^{MA}d{\Omega }_{2D}^{MA}}+{\displaystyle {\int }_{{\Omega }_{2D}^{MA}}\tau w_{2D}^{MA}d{\Gamma }_{2D}^{MA}}(32)$

where p is the body force, and τ is the traction of the boundary ${\Gamma }_{\mathrm{2D}}^{MA}$ in the 2D macroscale.

Sensitivity Analysis

Based on the topology optimization formulation, the sensitivities can be obtained via the first-order derivatives of the objective functions with respect to the design variables ${\alpha }_{2D}^{MI}$.

The sensitivity of the macroscopic compliance is computed by:

$\frac{\partial J^{MA}}{\partial {\alpha }_{2D}^{MI}}=\frac{1}{2}{\displaystyle {\int }_{{\Omega }_{2D}^{MA}}{\epsilon }_{ij}^T\left(u_{2D}^{MA}\right)\frac{\partial D_{ijkl}^H\left({\alpha }_{2D}^{MI}\right)}{\partial {\alpha }_{\xi }^{MI}}{\epsilon }_{kl}\left(u_{2D}^{MA}\right)d{\Omega }_{2D}^{MA}}(33)$

where the first-order derivatives of the effective elasticity tensor $D_{ijkl}^H$ with respect to t is:

$\frac{\partial D_{ijkl}^H}{\partial t}=\frac{1}{\left|{\Omega }_{2D}^{MI}\right|}{\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\beta }_{2D}^{MI}{v}^n\left|{\left(\nabla {\Phi }_{2D}^{MI}\right)}^T\right|\delta \left({\Phi }_{2D}^{MI}\right)d{\Omega }_{2D}^{MI}}(34)$

${\beta }_{2D}^{MI}=\left({\epsilon }_{pq}^{0\left(ij\right)}-{\epsilon }_{pq}^{\ast }\left(u_{2D}^{MI\left(ij\right)}\right)\right)D_{pqrs}\left({\epsilon }_{rs}^{0\left(kl\right)}-{\epsilon }_{rs}^{\ast }\left(u_{2D}^{MI\left(kl\right)}\right)\right)(35)$

where vⁿ is determined in Eq. 11 and is then substituted into Eq. 34:

$\frac{\partial D_{ijkl}^H}{\partial t}=\left(\frac{1}{\left|{\Omega }_{2D}^{MI}\right|}{\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\beta }_{2D}^{MI}{\phi }_{2D}^{MI}{\left(x\right)}^T\delta \left({\Phi }_{2D}^{MI}\right)d{\Omega }_{2D}^{MI}}\right){\dot{\alpha }}_{2D,n}^{MI}\left(t\right)(36)$

In terms of the chain rule, the first-order derivatives of DH ijkl with respect to t can also be given as:

$\frac{\partial D_{ijkl}^H}{\partial t}=\frac{\partial D_{ijkl}^H}{\partial {\alpha }_{2D}^{MI}}{\dot{\alpha }}_{2D,n}^{MI}\left(t\right)(37)$

Comparing Eq. 36 and Eq. 37, the first-order derivatives of $D_{ijkl}^H$ with respect to the expansion coefficients ${\alpha }_{2D}^{MI}$ can then be obtained as:

$\frac{\partial D_{ijkl}^H}{\partial {\alpha }_{2D}^{MI}}=\frac{1}{\left|{\Omega }_{2D}^{MI}\right|}{\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\beta }_{2D}^{MI}{\phi }_{2D}^{MI}{\left(x\right)}^T\delta \left({\Phi }_{2D}^{MI}\right)d{\Omega }_{2D}^{MI}}(38)$

Therefore, the sensitivity of the objective J^MA can be calculated by substituting Eq. 38 into Eq. 33.

Since the sensitivity of the objective JMI A is also based on Eq. 38, it can be calculated by:

$\frac{\partial J_A^{MI}}{\partial {\alpha }_{2D}^{MI}}=\frac{\partial {\left(Mu1+1\right)}^2}{\partial {\alpha }_{2D}^{MI}}+\frac{\partial {\left(Mu2+1\right)}^2}{\partial {\alpha }_{2D}^{MI}}=\frac{\partial {\left(D_{12}^H/D_{11}^H+1\right)}^2}{\partial {\alpha }_{2D}^{MI}}+\frac{\partial {\left(D_{12}^H/D_{22}^H+1\right)}^2}{\partial {\alpha }_{2D}^{MI}}(39)$

The sensitivity of the MFP can be derived from the Darcy-Stokes coupling system, so the first-order derivatives of $J_P^{MI}$ with respect to the expansion coefficients ${\alpha }_{2D}^{MI}$ can be given by:

$\frac{\partial J_P^{MI}}{\partial {\alpha }_{2D}^{MI}}=\frac{\partial {\mathbf{K}}_{2\text{D}}^{\text{H}\left(\mathrm{2,2}\right)}\left({\alpha }_{2D}^{MI}\right)}{\partial {\alpha }_{2D}^{MI}}(40)$

where the sensitivity of $\mathbf{K}$ can be further obtained according to the first order of the Heaviside function H^f (Dirac function) (Luo, et al., 2007) in the fluid-solid coupling system by

$\frac{\partial {\mathbf{K}}_{2\text{D}}^{\text{H}}\left({\alpha }_{2D}^{MI}\right)}{\partial {\alpha }_{2D}^{MI}}=\frac{1}{\left|{\Omega }_{2D}^{MI}\right|}{\displaystyle {\int }_{{\Omega }_{2D}^{MI}}{\left({\mathbf{w}}_{2D}^{\left(i\right)}\right)}^T{\phi }_{2D}^{MI}{\left(x\right)}^T{\delta }^f\left({\Phi }_{2D}^{MI}\right)\left({\mathbf{K}}_{\text{d}}^{2\text{D}}-{\mathbf{K}}_{\text{s}}^{2\text{D}}\right){\mathbf{w}}_{2D}^{\left(j\right)}d{\Omega }_{2D}^{MI}}(41)$

As discussed previously, the volume constraint is only defined in the 2D design domain but can control the volume of the 3D model. Thus, the sensitivity of the volume constraint can be calculated by:

$\frac{\partial V}{\partial {\alpha }_{2D}^{MI}}={\displaystyle {\int }_{{\Omega }_{\xi }^{MI}}{\phi }_{2D}^{MI}{\left(x\right)}^T\delta \left({\mathit{\Phi }}_{2D}^{MI}\left({\alpha }_{2D}^{MI}\right)\right)d{\mathit{\Omega }}_{2D}^{MI}}(42)$

After that, the non-dimensional sensitivities for the three objective functions can be given by Eq. 43, which benefits the multi-objective optimization (Marler and Arora, 2004).

$\frac{\partial J_i^{Non}}{\partial \alpha }=\frac{\partial J_i/\partial \alpha }{\left|\partial J_i^{\max }/\partial \alpha \right|}i\in \left(\mathrm{1,2,3}\right)(43)$

where i denotes the three individual objective functions.

Numerical Results

The numerical optimization is implemented via MATLAB 2018b to obtain the stenting microstructures. The multi-domain involved in the optimization model is displayed in Figure 4. From the figure, we can see three domains are defined in the model. The stent is deployed as a 2D structure, and it is the macro design domain that consists of 20 × 20 uniform microstructures. The macroscopic boundary and loads conditions: the vertical degree of freedoms of the top and bottom edges are fixed, and horizontal unit displacements are applied on the left and right edges. The stenting microstructures is then defined as the micro design domain. The 2D micro design domain can be treated as the constant horizontal cross-section profile of a 3D microstructure, which is utilized to establish a 3D micro computational domain for the related CFD analysis. In the 3D field, a scaled velocity of blood flow is specified as the inlet condition, while the zero-pressure boundary is defined as the outlet condition. By considering the computational efficiency, a square element with four nodes is used for the 2D micro and macro design domains; an 8-node brick element is used in the 3D domain.

FIGURE 4

FIGURE 4. Multi-domain of the numerical model.

The multi-objective optimization is formulated by weighting three single objective functions. Different combinations of weight factors can lead to different results. All the three optimization objectives are critical to the design. Compared with the macroscopic compliance, the other two objectives determine the material layout to have auxetic and MFP properties, so larger weighting factors are specified for these two objective functions. Thus, three weight factors, 0.35, 0.35 and 0.30, are used.

The volume fraction of the microstructure in the optimization is another factor in determining the optimization results. It should not be too big to block blood flow, or too small to support the artery. Thus, the volume fractions from 25 to 60% are used to investigate the optimization results. The optimized results are summarized in Table 1, which includes the microstructure of each case, the related 3 × 3 microstructural array, and the blood flow streamlines in the design domain. In the table, the black colour denotes stenting microstructures, and the red colour shows blood.

TABLE 1

TABLE 1. The results of various volume fractions.

All results in Table 1 show auxetic characteristics: re-entrant connections between struts. The solid material in all designs concentrate more along the direction of the blood flow, which is beneficial for reducing the stenting obstructions. To distinguish different struts in the stenting microstructures, the struts along the blood flow are named horizontal struts, while the others are called vertical struts. When the designs have relatively large volume fractions, such as the cases 1–4, more solid material will be assigned to the horizontal struts to facilitate the blood flow, leading to less or thin vertical struts. When volume fractions are reduced, such as the cases 5–7, the material distributions become more uniform. In these cases, the two horizontal struts are connected by only one vertical strut location on the right edge, exhibiting an opening at the upstream side. This can lead to a long space for blood through and facilitate recovering the disturbed flow. Besides that, the other two vertical struts are also near the right edge to facilitate blood flow. As the volume fraction was further reduced in case 8, the two vertical struts near the horizontal struts become thinner. They are used to connect the stenting unit cells in the same circumference, and two weak connections may lead to stenting fractures. Within the results, case 7 with a 30% volume fraction has more material distribution and few obstructions on the blood flow.

The homogenized effective elasticity tensor and the MFP are given by

$D^H=\left[\begin{array}{ccc}\text{0}\text{.}1082& -\text{0}\text{.}0339& 0\\ -\text{0}\text{.}0339& \text{0}\text{.}0324& 0\\ 0& 0& \text{0}\text{.}0022\end{array}\right]\text{where},\{\begin{array}{c}Mu1=-0.3133\\ Mu2=-1.0463\end{array}(44)$

and

${\mathbf{K}}_{2\text{D}}^{\text{H}}=\left[\begin{array}{cc}1\text{.}0593\times {10}^{-3}& 0\\ 0& 2\text{.}5233\times {10}^{-4}\end{array}\right](45)$

It can be found that the values of NPRs are −0.3133 along the flow direction and −1.0463 in the vertical direction, which indicates the microstructure has a relatively smaller vertical stiffness. It means the stenting radial stiffness is smaller than the axial stiffness by adopting the optimized microstructure. The less radial stiffness of the stent can enhance the adaptive deformation and lead to better flexibility. On the other hand, both directions present the auxetic properties that can obtain a smaller volume after compression, so as to facilitate deliverability. The effective MFP can be found in Eq. 45. With respect to the MFP, the design reduces the vertical blood flow which can also leads to directional differences. Hence, the directional differences in both NPR and MFP improve the stenting performances.

To further investigate numerical case 7, the dynamic evolutions of the microstructure are presented in Figure 5. We can see the solid material is gradually assigned to the horizontal struts in the optimization process, leading to fewer obstructions to the blood flow. Re-entrant characteristics exhibiting auxetic behaviour can also be created, simultaneously. As shown in Figure 5F, the blood flow is mainly disturbed by the vertical struts, but the normal flow patterns can be recovered behind the struts due to large fluid spaces. Thus, local recirculation regions may occur around the vertical struts but will be small. The convergences of the related objective functions in the optimization are shown in Figure 6. All effective properties, objective functions, and volume constraints are gradually converged to their solutions. We can see that the two MFPs increase from zero to the highest value at the beginning of the optimization, because the solid material fills the design domain at the initial stage. In a word, the optimization process of the case 7 also shows the effectiveness of the optimization.

FIGURE 5

FIGURE 5. The results of case 7: (A) Initial design; (B–E) Four intermediate results; (F) Final design.

FIGURE 6

FIGURE 6. The convergent history of case 7: (A) Poisson’s ratios; (B) Modified permeabilities; (C) Objective and volume fraction.

Numerical Validations

This section is to numerically validate the performances of the optimized stent. For instance, the minimization of the stent induced adverse hemodynamic can be obtained by reducing the obstructions of the stent to the blood flow. WSS and related derivatives are usually used as the widely accepted metric to measure the stented hemodynamics. Thus, the optimized microstructure is simulated in ANSYS 2019R3 to validate the macroscopic performances. For the optimized stent, the auxetic properties, disturbances to the blood flow, and WSS distributions are numerical validated, respectively.

Auxetics

The auxetic property is one of the design objectives for the stent in the proposed design. Due to the directional differences of the stenting properties, two compression tests are performed: one is under a radial compression load, and the other is under an axial compression load. To create the stenting geometry for the tests, the optimized microstructures are first obtained from MATALB via the STL format, and then ANSYS SpaceClaim module is utilized to finish the creation of the geometry (Figure 7). The stent thickness is 0.1 mm, and the length and the diameter are 12 and 4 mm. In the model, 20 cells are configured along the peripheral direction and 20 cells are arranged along the axial direction.

FIGURE 7

FIGURE 7. The geometry of the optimized stent.

The deformation of the radial compression test is shown in Figures 8A,B, where both axonometric and side views are utilized to show its deformation behaviours. A pressure load is applied on the stenting surface to compress the diameter about 50% in the test. From the figure, we can see the stent contracts slightly in the axial direction due to larger axial stiffness. When the stent is compressed in the radial direction to adapt to the various artery shapes, the stenting axial length is relatively stable. It can keep enough length of the stent to cover the lesion even though the stenting radial dimension is changed. The axial compression test result is presented in Figures 8C,D. The test result shows that the stent contracts significantly in the radial direction when applying compression loads in the axial direction from the right side, indicating the NPR behaviour that can improve deliverability.

FIGURE 8

FIGURE 8. The deformation results of the compression test: (A) Axonometric view of the radial compression test result; (B) Side view of the radial compression test result; (C) Axonometric view of the axial compression test result; (D) Side view of the axial compression test result.

Fluid Validations

The disturbance of stents on blood flow is another essential issue in the design. Thus, the simulation for blood flow in the stented segment is performed via the software CFX in ANSYS v2019R3. This research intends to reduce the adverse hemodynamic changes by minimizing the stenting obstructions on the blood flow. In the validations, the streamlines generated by the blood flow are used to investigate the stenting obstructions. WSS distribution is also computed in the simulation.

To save computational cost, the CFD model is established based on 1 cell along the peripheral direction while 20 cells in the axial direction, as shown in Figure 9, where the black one is the stent, and the red is blood. The stent thickness is 0.1 mm, and the length and the diameter are 12 and 4 mm, respectively. By considering the impact of the boundary conditions on flow, the distance between the inlet and the stent struts is specified as 3 mm and another 3 mm between the outlet and the stent. The blood in the model is defined as an incompressible Newtonian viscous fluid. The density is 1,050 kg/m3, and the dynamic viscosity is 3.5 mPa⋅s. The flow is assumed to be steady state with 0.35 m/s inflow velocity and 0 Pa outlet boundary condition. No-slip conditions are applied on the fluid-solid interfaces and the fluid walls. The cylindrical symmetric condition is defined in the peripheral direction.

FIGURE 9

FIGURE 9. The CFD model.

The results of blood flow streamlines are shown in Figure 10, where the lower velocities are denoted with darker colours. In the figure, the maximum velocity of the flow is greater than the inflow condition, because the flow can be accelerated around the stent struts. Due to the applied no-slid boundary conditions, the minimum velocity is zero. The result shows that the velocities of the streamlines are affected by the stenting struts. Since relatively narrow spaces are formed between the horizontal struts in each unit cell, the blood flows are obstructed around these spaces and exhibits lower velocities. And the blood flows between each unit cell can maintain higher speeds due to fewer obstructions and larger fluid areas. Overall, blood flow can keep a high velocity, which can be demonstrated that many streamlines present light colours through the whole fluid region.

FIGURE 10

FIGURE 10. The streamlines of blood around the stent F-7.

The overall influences of stenting obstructions can be evaluated by blood flow changes in the upstream side and downstream side of the stent. The streamlines in the proximal and distal struts are presented in Figure 11. In the figure, the fluid region around the stent can be divided into two different parts: one is the inner space of each individual cell, and the other is the connecting space between cells. When blood through these spaces, the flow can be accelerated around the struts and may lead to recirculation zones. In the connecting spaces of the proximal struts, the disturbances of the flow focus on the four corners and cause small recirculation zones. Also, a few vertical flows occur along the left struts in the connecting spacings. However, small recirculation zones are formed around the vertical struts in the inner spaces of proximal struts due to the narrow fluid spaces. Compared with the proximal struts, similar blood flow behaviours can be found around the distal struts but exhibit relatively smaller flow disturbances and recirculation zones. Although the recirculation zones may induce adverse effects, they only cover very few regions around the struts. Hence, the result shows that the optimization for the MFP can effectively reduce the vertical flow and the stenting obstructions to benefit blood flow.

FIGURE 11

FIGURE 11. The streamlines of blood around the proximal and distal struts: (A) Proximal streamlines (Axonometric view); (B) Proximal streamlines (Side view); (C) Distal streamlines (Axonometric view); (D) Distal streamlines (Side view).

The WSS is an important metric to quantify the hemodynamic impacts in stent implantation. The unusual WSS is often associated with adverse effects including ST (extremely high WSS) and ISR (extremely low WSS). When the magnitude of WSS is higher than 2.5 Pa increases the risk of ST, while lower than 0.5 Pa increases ISR incidence. Based on the flow direction, the WSS is mainly determined by the axial component, so the axial WSS is utilized to evaluate the stent. The axial WSS of the stent is illustrated in Figure 12. It can be found that the regions covered by the stent have higher WSS. In the results, the maximum axial WSS is 2.195 Pa, which is less than the high adverse stress 2.5 Pa and indicates a lower ST incidence. However, the regions surrounded by the stent struts exhibit lower and even negative axial WSS, which are recirculation zones. The same recirculation zones are also found in the streamline results. Among them, the connecting spaces between unit cells have higher WSS compared with their inner spaces due to larger fluid regions. The average WSS in most regions of the stent is around 0.9 Pa and greater than the lower bound of the threshold 0.5 Pa, but the areas close to the vertical struts exhibit smaller, even negative WSS. It is because the struts protrude into the artery lumen and separate the laminar flow, resulting in lower shear stress around the struts. This kind of effect is determined by the shape and size of the strut profiles and cannot be eliminated unless the stent is fully embedded into the artery walls. The regions with lower WSS cannot be completely avoided, but the impacted areas are small. Moreover, the WSS distribution of the stent is relatively uniform, which can avoid structural fatigue breakage to extend stent lifetime in use and is able to low occurrence of ST and ISR.

FIGURE 12

FIGURE 12. The WSS distributions of the stent: (A) Covered by the stent; (B) Without the stent.

Hence, from the above numerical examples and discussions we only found limited disturbance in blood flow around the area of the stent, which indicates less obstructions of the newly designed stent to the blood flow. This shows the efficiency of the optimally designed structural architectures in stents.

Conclusion

This paper has developed a family of auxetic architectures for stenting structures using a topological optimization approach and its associated numerical methods, which is especially beneficial to self-expanding stents. The new design combines the hemodynamic effect with auxetic structures to improve the stent performances particularly from both mechanical and hemodynamic aspects. Several numerical examples and their simulation results have shown the efficiency of the topologically optimized stenting architectures. The homogenized MFP can quantify the disturbances of the stenting microstructures to the blood flow, through which the reduction of the adverse blood flows can be transformed to the reduction of the obstructions of the stent to the blood flow. The obstructions have been found to associate with the material layouts of the microstructures. The fluent simulation results present fewer obstructions of the optimized stent to the blood flow and fewer adverse WSS distributions associated with ST and ISR risk factors. Hence, the newly developed auxetic architectures can benefit stenting performance by reducing occurrence of mechanical failures as well as the influence of adverse hemodynamics, which will help low the incidence of ST and ISR complications due to mechanical structural issues in stent implantation and treat heart disease in clinic practice.

In this study, it is noted that the pulsatile blood flow state is not considered in the design. Furthermore, this paper is focused on generic design of a new type of stenting architectures sharing common geometric profiles, not patient-specific and customised designs, through a topological optimisation method using idealized models. The success rate of PCI therapy has been significantly improved over the past, but the ST and ISR problems have not been completely resolved. Various factors may account for the incidence of ST and ISR in stent implantation. However, this work emphases the importance of mechanical structural failures and hemodynamics, rather than biomaterial, patient, clinical and operational aspects. Moreover, technical details for developing end-user products such as compressing stents into sheaths, injection and release of stents in implantation are outside of the scope of this paper.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

HX undertook the research and completed most detailed contents of this manuscript. SS was provided useful advice and assistance particularly CFD simulation. SB given useful discussions in wear shearing stress and hemodynamic. NJ contributed professional advice for stenting design and implantation. ZL was the principal supervisor who has supplied the major research ideas and financial support and helped HX complete this manuscript. All authors have made solid contributions to the completion of this manuscript.

Funding

This research was partially supported by the Australian Research Council (ARC)—Discovery Projects (DP210101353).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Babuška, I. (1971). Error-Bounds for Finite Element Method. Numer. Math. 16 (4), 322–333. doi:10.1007/BF02165003

CrossRef Full Text | Google Scholar

Balossino, R., Gervaso, F., Migliavacca, F., and Dubini, G. (2008). Effects of Different Stent Designs on Local Hemodynamics in Stented Arteries. J. Biomech. 41 (5), 1053–1061. doi:10.1016/j.jbiomech.2007.12.005

CrossRef Full Text | Google Scholar

Beier, S., Ormiston, J., Webster, M., Cater, J., Norris, S., Medrano-Gracia, P., et al. (2016). Hemodynamics in Idealized Stented Coronary Arteries: Important Stent Design Considerations. Ann. Biomed. Eng. 44 (2), 315–329. doi:10.1007/s10439-015-1387-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bekal, C., Shetty, R., and Shenoy, S. (2018). Numerical Investigation of Influence of Number of Stent Cells and Type of Link on Expansion and Haemodynamic Behaviour of Balloon-Expandable Coronary Stent. Sādhanā. 43 (6), 1–12. doi:10.1007/s12046-018-0883-1

CrossRef Full Text | Google Scholar

Berry, J. L., Santamarina, A., Moore, Jr., J. E., Roychowdhury, S., and Routh, W. D. (2000). Experimental and Computational Flow Evaluation of Coronary Stents. Ann. Biomed. Eng. 28 (4), 386–398. doi:10.1114/1.276

PubMed Abstract | CrossRef Full Text | Google Scholar

Bressloff, N. W., Ragkousis, G., and Curzen, N. (2016). Design Optimisation of Coronary Artery Stent Systems. Ann. Biomed. Eng. 44 (2), 357–367. doi:10.1007/s10439-015-1373-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Bukač, M., Čanić, S., Tambača, J., and Wang, Y. (2019). Fluid-Structure Interaction between Pulsatile Blood Flow and a Curved Stented Coronary Artery on a Beating Heart: A Four Stent Computational Study. Computer Methods Appl. Mech. Eng. 350, 679–700. doi:10.1016/j.cma.2019.03.034

CrossRef Full Text | Google Scholar

Doutel, E., Carneiro, J., Campos, J. B. L. M., and Miranda, J. M. (2018). Experimental and Numerical Methodology to Analyze Flows in a Coronary Bifurcation. Eur. J. Mech. - B/Fluids. 67, 341–356. doi:10.1016/j.euromechflu.2017.09.009

CrossRef Full Text | Google Scholar

Garasic, J. M., Edelman, E. R., Squire, J. C., Seifert, P., Williams, M. S., and Rogers, C. (2000). Stent and Artery Geometry Determine Intimal Thickening Independent of Arterial Injury. Circulation. 101 (7), 812–818. doi:10.1161/01.cir.101.7.812

PubMed Abstract | CrossRef Full Text | Google Scholar

Gharleghi, R., Wright, H., Khullar, S., Liu, J., Ray, T., and Beier, S. (2019). “Advanced Multi-Objective Design Analysis to Identify Ideal Stent Design,” in Machine Learning and Medical Engineering for Cardiovascular Health and Intravascular Imaging and Computer Assisted Stenting (Springer), 193–200. doi:10.1007/978-3-030-33327-0_23

CrossRef Full Text | Google Scholar

Gharleghi, R., Wright, H., Luvio, V., Jepson, N., Luo, Z., Senthurnathan, A., et al. (2021). A Multi-Objective Optimization of Stent Geometries. J. Biomech. 125, 110575. doi:10.1016/j.jbiomech.2021.110575

CrossRef Full Text | Google Scholar

Gijsen, F., Katagiri, Y., Barlis, P., Bourantas, C., Collet, C., Coskun, U., et al. (2019). Expert Recommendations on the Assessment of Wall Shear Stress in Human Coronary Arteries: Existing Methodologies, Technical Considerations, and Clinical Applications. Eur. Heart J. 40 (41), 3421–3433. doi:10.1093/eurheartj/ehz551

PubMed Abstract | CrossRef Full Text | Google Scholar

Guest, J. K., and Prévost, J. H. (2006). Optimizing Multifunctional Materials: Design of Microstructures for Maximized Stiffness and Fluid Permeability. Int. J. Sol. Structures. 43 (22-23), 7028–7047. doi:10.1016/j.ijsolstr.2006.03.001

CrossRef Full Text | Google Scholar

Gundert, T. J., Marsden, A. L., Yang, W., and LaDisa, J. F. (2012). Optimization of Cardiovascular Stent Design Using Computational Fluid Dynamics. J. Biomech. Eng. 134 (1), 011002. doi:10.1115/1.4005542

CrossRef Full Text | Google Scholar

Hughes, T. J. R., Franca, L. P., and Balestra, M. (1986). A New Finite Element Formulation for Computational Fluid Dynamics: V. Circumventing the Babuška-Brezzi Condition: a Stable Petrov-Galerkin Formulation of the Stokes Problem Accommodating Equal-Order Interpolations. Computer Methods Appl. Mech. Eng. 59 (1), 85–99. doi:10.1016/0045-7825(86)90025-3

CrossRef Full Text | Google Scholar

Johnston, B. M., Johnston, P. R., Corney, S., and Kilpatrick, D. (2006). Non-Newtonian Blood Flow in Human Right Coronary Arteries: Transient Simulations. J. Biomech. 39 (6), 1116–1128. doi:10.1016/j.jbiomech.2005.01.034

CrossRef Full Text | Google Scholar

Ku, D. N. (1997). Blood Flow in Arteries. Annu. Rev. Fluid Mech. 29 (1), 399–434. doi:10.1146/annurev.fluid.29.1.399

CrossRef Full Text | Google Scholar

LaDisa, J. F., Olson, L. E., Molthen, R. C., Hettrick, D. A., Pratt, P. F., Hardel, M. D., et al. (2005). Alterations in Wall Shear Stress Predict Sites of Neointimal Hyperplasia after Stent Implantation in Rabbit Iliac Arteries. Am. J. Physiology-Heart Circulatory Physiol. 288 (5), H2465–H2475. doi:10.1152/ajpheart.01107.2004

CrossRef Full Text | Google Scholar

Lakes, R. (1993). Advances in Negative Poisson's Ratio Materials. Adv. Mater. 5 (4), 293–296. doi:10.1002/adma.19930050416

CrossRef Full Text | Google Scholar

Li, H., Luo, Z., Gao, L., and Walker, P. (2018). Topology Optimization for Functionally Graded Cellular Composites With Metamaterials by Level Sets. Computer Methods Appl. Mech. Eng. 328, 340–364. doi:10.1016/j.cma.2017.09.008

CrossRef Full Text | Google Scholar

Luo, Z., Tong, L., and Kang, Z. (2009). A Level Set Method for Structural Shape and Topology Optimization Using Radial Basis Functions. Comput. Structures. 87 (7-8), 425–434. doi:10.1016/j.compstruc.2009.01.008

CrossRef Full Text | Google Scholar

Luo, Z., Tong, L., Wang, M. Y., and Wang, S. (2007). Shape and Topology Optimization of Compliant Mechanisms Using a Parameterization Level Set Method. J. Comput. Phys. 227 (1), 680–705. doi:10.1016/j.jcp.2007.08.011

CrossRef Full Text | Google Scholar

Marler, R. T., and Arora, J. S. (2004). Survey of Multi-Objective Optimization Methods for Engineering. Struct. Multidisciplinary Optimization. 26 (6), 369–395. doi:10.1007/s00158-003-0368-6

CrossRef Full Text | Google Scholar

Masud, A., and Hughes, T. J. (2002). A Stabilized Mixed Finite Element Method for Darcy Flow. Comput. Methods Appl. Mech. Eng. 191 (39-40), 4341–4370. doi:10.1016/S0045-7825(02)00371-7

CrossRef Full Text | Google Scholar

Mattace-Raso, F. U. S., van der Cammen, T. J. M., Hofman, A., van Popele, N. M., Bos, M. L., Schalekamp, M. A. D. H., et al. (2006). Arterial Stiffness and Risk of Coronary Heart Disease and Stroke. Circulation. 113 (5), 657–663. doi:10.1161/CIRCULATIONAHA.105.555235

PubMed Abstract | CrossRef Full Text | Google Scholar

Melzer, H.-S., Ahrens, R., Guber, A. E., and Dohse, J. (2019). Numerical Simulation and In Vitro Examination of the Flow Behaviour Within Coronary Stents. Curr. Dir. Biomed. Eng. 5 (1), 541–544. doi:10.1515/cdbme-2019-0136

CrossRef Full Text | Google Scholar

Nolan, D. R., and Lally, C. (2018). An Investigation of Damage Mechanisms in Mechanobiological Models of In-Stent Restenosis. J. Comput. Sci. 24, 132–142. doi:10.1016/j.jocs.2017.04.009

CrossRef Full Text | Google Scholar

Putra, N. K., Palar, P. S., Anzai, H., Shimoyama, K., and Ohta, M. (2019). Multiobjective Design Optimization of Stent Geometry With Wall Deformation for Triangular and Rectangular Struts. Med. Biol. Eng. Comput. 57 (1), 15–26. doi:10.1007/s11517-018-1864-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Rogers, C., and Edelman, E. R. (1995). Endovascular Stent Design Dictates Experimental Restenosis and Thrombosis. Circulation. 91 (12), 2995–3001. doi:10.1161/01.cir.91.12.2995

PubMed Abstract | CrossRef Full Text | Google Scholar

Stiehm, M., Wüstenhagen, C., Siewert, S., Ince, H., Grabow, N., and Schmitz, K.-P. (2019). Impact of Strut Dimensions and Vessel Caliber on Thrombosis Risk of Bioresorbable Scaffolds Using Hemodynamic Metrics. Biomed. Tech. (Berl). 64 (3), 251–262. doi:10.1515/bmt-2017-0101

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomaszewski, M., Sybilski, K., Małachowski, J., Wolański, W., and Buszman, P. P. (2020). Numerical and Experimental Analysis of Balloon Angioplasty Impact on Flow Hemodynamics Improvement. Acta Bioeng. Biomech. 22 (3), 169. doi:10.37190/ABB-01660-2020-03

PubMed Abstract | CrossRef Full Text | Google Scholar

Traub, O., and Berk, B. C. (1998). Laminar Shear Stress. Atvb. 18 (5), 677–685. doi:10.1161/01.ATV.18.5.677

CrossRef Full Text | Google Scholar

Vijayaratnam, P. R. S., Reizes, J. A., and Barber, T. J. (2019). Flow-Mediated Drug Transport From Drug-Eluting Stents Is Negligible: Numerical and In-Vitro Investigations. Ann. Biomed. Eng. 47 (3), 878–890. doi:10.1007/s10439-018-02176-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M. Y., Wang, X., and Guo, D. (2003). A Level Set Method for Structural Topology Optimization. Comput. Methods Appl. Mech. Eng. 192 (1-2), 227–246. doi:10.1016/S0045-7825(02)00559-5

CrossRef Full Text | Google Scholar

Wang, Y., Luo, Z., Zhang, N., and Kang, Z. (2014). Topological Shape Optimization of Microstructural Metamaterials Using a Level Set Method. Comput. Mater. Sci. 87, 178–186. doi:10.1016/j.commatsci.2014.02.006

CrossRef Full Text | Google Scholar

Wang, Y., Luo, Z., Zhang, N., and Qin, Q. (2016). Topological Shape Optimization of Multifunctional Tissue Engineering Scaffolds With Level Set Method. Struct. Multidisc Optim 54 (2), 333–347. doi:10.1007/s00158-016-1409-2

CrossRef Full Text | Google Scholar

Wei, L., Leo, H. L., Chen, Q., and Li, Z. (2019). Structural and Hemodynamic Analyses of Different Stent Structures in Curved and Stenotic Coronary Artery. Front. Bioeng. Biotechnol. 7, 366. doi:10.3389/fbioe.2019.00366

PubMed Abstract | CrossRef Full Text | Google Scholar

Wentzel, J. J., Krams, R., Schuurbiers, J. C. H., Oomen, J. A., Kloet, J., van der Giessen, W. J., et al. (2001). Relationship Between Neointimal Thickness and Shear Stress after Wallstent Implantation in Human Coronary Arteries. Circulation. 103 (13), 1740–1745. doi:10.1161/01.CIR.103.13.1740

PubMed Abstract | CrossRef Full Text | Google Scholar

Whitaker, S. (1986). Flow in Porous Media I: A Theoretical Derivation of Darcy's Law. Transp Porous Med. 1 (1), 3–25. doi:10.1007/BF01036523

CrossRef Full Text | Google Scholar

Wu, J., Luo, Z., Li, H., and Zhang, N. (2017). Level-set Topology Optimization for Mechanical Metamaterials Under Hybrid Uncertainties. Computer Methods Appl. Mech. Eng. 319, 414–441. doi:10.1016/j.cma.2017.03.002

CrossRef Full Text | Google Scholar

Xue, H., Luo, Z., Brown, T., and Beier, S. (2020). Design of Self-Expanding Auxetic Stents Using Topology Optimization. Front. Bioeng. Biotechnol. 8, 736. doi:10.3389/fbioe.2020.00736

PubMed Abstract | CrossRef Full Text | Google Scholar