Numerical study on the impulsive growth of a gaseous plug inside a cylindrical vein with compliant coating

Geometry of the problem

The targeted vein of a human body within which the dynamic behavior of the explosive growth of a gaseous plug has to be investigated was assumed to have a cylindrical shape. At the first step the cylindrical vein was considered as a rigid wall, and it has a rigid wall with an internal compliant coating. Fig. 1 illustrates the schematic presentation of a rigid cylindrical vein, while Fig. 2 illustrates a rigid cylindrical vein with an internal compliant coating. The axis of symmetry is indicated by z and is along the length of the cylinder, while r represents the radial axis.

Fig. 1

Schematic presentation of a rigid cylindrical vein.

×

Fig. 1.

Fig. 2

Schematic presentation of a rigid cylindrical vein with an internal compliant coating.

×

Fig. 2.

Hydrodynamic equations of the problem

The liquid domain inside the cylindrical vein has been assumed to be inviscid, incompressible, and irrotational. The liquid around the bubble is actually blood if the bubble is gaseous plug inside a vein of the human body. But in this research for simplicity, the liquid around the bubble was assumed as water and then its viscosity was very low and could be neglected. Then by considering the liquid flow around the gaseous plug as an incompressible irrotational flow, we can assume the liquid around the gaseous plug as potential flow and employ the Green’s integral formula and boundary element method. The assumption of potential flow around and inside the veins of the human body in complex problems was an acceptable assumption.¹⁵ Therefore, the liquid flow around the gaseous plug was assumed to be a potential flow. Thus, the Green’s integral formula was given as the governing hydrodynamic equation of the problem, where S_gp the boundary of the gaseous was plug and S_cw was the internal wall of the cylindrical vein:

\[c(p)\varphi (p) + \int\limits_{{S_{gp}} + {S_{cw}}}^{} {\varphi (p)\frac{\partial }{{\partial n}}} \left( {\frac{1}{{\left| {p - q} \right|}}} \right)dS = \int\limits_{{S_{gp}} + {S_{cw}}}^{} {\frac{\partial }{{\partial n}}(\varphi (q))\left( {\frac{1}{{\left| {p - q} \right|}}} \right)} dS\] (Eq. 1)

The internal wall of the cylindrical vein has been continued up to the physical infinity, where the explosive growth of the bubble has negligible effect on the hydrodynamic behavior of the liquid domain around the bubble. In equation (1), φ is the velocity potential, p is any point on the boundary or inside of the liquid domain, q is any point on the boundary, and n is the normal direction to the boundary and

\[c(p) = \left\{ \begin{gathered} 2\pi ,if\;p\;is\;on\;the\;boundary\;including\;the \hfill \\ wall\;and\;the\;bubble\;surfrace \hfill \\ 4\pi ,if\;p\;is\;inside\;the\;liquid\;domain \hfill \\ \end{gathered} \right.\]

The unsteady Bernoulli equation provides the possibility for evaluation of the velocity potentials on the boundaries of the liquid domain at any instant as time advances and is given in the Eulerian form as:

\[\frac{{\partial \varphi }}{{\partial t}} = \frac{{{P_W} - {P_\infty }}}{\rho } - \frac{1}{2}{\left| {\nabla \varphi } \right|^2}\] (Eq. 2)

Since in this problem the bubble size was very small, the buoyancy forces were negligible. The Lagrangian form of the unsteady Bernoulli equation which must be employed on the moving boundaries is given as:

\[\frac{{D\varphi }}{{Dt}} = \frac{{{P_W} - {P_\infty }}}{\rho } - \frac{1}{2}{\left| {\nabla \varphi } \right|^2}\] (Eq. 3)

In equations (2) and (3), P_w is the pressure on the internal wall of the cylindrical vein, P_∞ is the pressure at the far field, ρ is the liquid density, g is the gravity acceleration, and r represents the distance from the axis of symmetry.

Mathematical modeling of the dynamic behavior of a compliant coating

For the purpose of more realistic mathematical modeling of a vein of the human body, a compliant surface was assumed as an internal coating of the rigid cylindrical vein. The compliant coating has been assumed to be a thin rigid plate-spring model and the mathematical equation of its dynamic behavior has been given as:

\[\frac{{{\partial ^2}\eta }}{{\partial {t^2}}} = \frac{{T{\partial ^2}\eta }}{{m\partial {x^2}}} - d\frac{{\partial \eta }}{{\partial t}} - \frac{B}{m}\frac{{{\partial ^4}\eta }}{{\partial {y^4}}} - \frac{K}{m}\eta + f\] (Eq. 4)

Equation (4) has been employed by Carpenter and Garrad¹⁶ and the other researchers^17-20 to simulate Kramer-type coatings. This relatively simple model contains characteristics of a broad range of surfaces. Equation (4) mathematically models the translational displacement of a Kramer type compliant coating and has been given in the Cartesian coordinates. This equation provides the geometrical data for that part of the half axisymmetric profile of the gaseous plug which is attached to the compliant coating. Once the geometrical shape of the gaseous plug’s axisymmetric half profile was obtained, and used in the axisymmetric form of the Green’s integral formula. T is the longitudinal tension and B is the flexural rigidity of the thin plate. Also, m is the mass per unit area, d is the damping coefficient, K is the spring constant and f is the external forcing term.

In this research, the compliant wall was assumed to be cylindrical, which was attached to the internal surface of a rigid cylinder. The cylindrical compliant coating was attached to the rigid boundary along its outer edges. By assuming a very thin membrane and by ignoring the effect of damping, the equation of the cylindrical compliant coating may be written in the form of:

\[\frac{{{\partial ^2}\eta }}{{\partial {t^2}}} = \frac{{T{\partial ^2}\eta }}{{m\partial {x^2}}} - \frac{k}{m}\eta + ({P_w}(z,t) - {P_\infty })\] (Eq. 5)

Where,η =η (r,t) is the translational displacement of the cylindrical compliant coating which is in the radial direction and towards or away from the axis of symmetry. P_w is the pressure on the compliant surface. In the absence of the bubble, for (t≤ 0) the springs are compressed uniformly due to the steady uniform pressure, P_∞ exerted by the internal liquid domain of the cylinder. The displacement of the thin cylindrical plate at this instant has been taken as zero. The cylindrical thin plate was attached to the rigid wall at z = ± z_cs, so the inward or outward radial displacements of the thin plate at z = ± z_cs were zero. During the explosive growth phase of the gaseous plug, the liquid domain inside the cylindrical vein and the surface of its internal compliant coating were coupled using the linearized equation for the pressure and the velocity in the both systems of the liquid domain inside the cylindrical vein and the surface of the compliant coating that attached to its internal wall. These equations which are satisfied at the undisturbed position of the internal surface of the cylindrical compliant coating,(r=r_cs), are:

\[\frac{{\partial \varphi }}{{\partial r}} = \frac{{\partial \eta }}{{\partial t}} = \psi \] (Eq. 6)

\[{P_W}(z,t) - {P_\infty } = - \rho \frac{{\partial \phi }}{{\partial t}}\] (Eq. 7)

It should be noted that equation (7) is the linearized unsteady Bernoulli equation.

Geometrical and mathematical discretization

The boundary of the gaseous plug was discretized by the cubic spline elements, whereas the internal surface of the cylindrical vein was discretized by the linear segments.

Fig. 3A illustrates the discretized geometry of the problem, when the cylindrical vein was assumed as a rigid cylinder. The z-axis along the cylinder length was the axis of symmetry and discretization on the internal wall of the rigid cylinder along the z-axis was continued from both left- and right sides of the gaseous plug until the points in the far field. The points in the far field are wherever that the explosive growth of the gaseous plug has negligible effects on the hydrodynamic behavior of the liquid domain.

Fig. 3

Schematic discretization of the geometry in the case of (A) gaseous plug inside a vertical rigid vein, (B) gaseous plug inside a vertical rigid vein with a compliant coating (hollow circles indicate the end points of the elements and black circles represent the colocation points).

×

Fig. 3.

Fig. 3B illustrates the discretized geometry of the problem when the vein has been assumed to be a rigid cylinder with an internal compliant coating. As it is stated in the case of Fig. 3A, the boundary of the gaseous plug has been discretized by the cubic spline elements and the internal wall of the cylindrical vein, which has a compliant coating around the explosive growing gaseous plug, was discretized by the linear segments. The height of the cylindrical compliant coating was 2.5 times of the maximum diameter of the gaseous plug and its outer edge was attached to the rigid wall from its both left- and right sides. The surface of the gaseous plug is discretized by the cubic spline elements and the internal wall of the cylindrical vein is discretized by the linear segment. The boundary conditions for the cubic spline elements over the gaseous plug surface were clamped boundary conditions and the first and last cubic spline elements over the half profile of the gaseous plug must be tangential to any line parallel to the axis of r. The geometrical characteristics of the problem have been non-dimensionalized by the equivalent maximum radius of the plug. Due to the significant importance of the discretization in the numerical simulations, the gird independency was considered.

As it can be seen in Fig. 4, the numerical procedure has been done for four different number of nodes in order to gain grid independence. Fig. 4 shows that the numerical results for the case of N=300 and N=380 are almost the same. On this base, discretization with N=300 nodes has been chosen as standard discretization and the problem has been solved with the discretization of 300 nodes. In the case of the rigid cylindrical vein, the grid independence, which has been discussed in Fig. 4, is enough for considering the stability of the problem. But in the case of having a compliant surface the size of the element on the compliant surface in the direction of z-axis becomes a more critical problem due to the Tolliman-Schilichting instability. In this case, Δz must be as small as possible in a manner that discretized compliant surface could follow the wavelength of the disturbed wavy compliant wall.

Fig. 4

Grid independency of the numerical simulation .

×

Fig. 4.

The hydrodynamic equations represent the mathematical models for the dynamic behaviors of the gaseous plug and the liquid domain in it’s around have been discretized as:

\[\begin{gathered} c({p_i})\varphi ({p_i}) + \hfill \\ \sum\nolimits_{j = 1}^{MN} {\left\{ {\varphi ({q_j})\int\limits_{{S_{gp}} + {S_{cw}}} {\frac{\partial }{{\partial n}}\left( {\frac{1}{{\left| {{p_i} - {q_j}} \right|}}} \right)} dS} \right\} = \sum\nolimits_{j = 1}^{MN} {\left\{ {\psi ({q_j})\int\limits_{{S_{gp}} + {S_{cw}}} {\left( {\frac{1}{{\left| {{p_i} - {q_j}} \right|}}} \right)} dS} \right\}} } , \hfill \\ \end{gathered} \] (Eq. 8)

\[\frac{{\varphi _j^{n + 1} - \varphi _j^n}}{{{t^{n + 1}} - {t^n}}} = \frac{{{P_W} - {P_\infty }}}{\rho } + \frac{1}{2}{\left| {\nabla {\varphi _j}} \right|^2}\] (Eq. 9)

\[\frac{{\varphi _j^{n + 1} - \varphi _j^n}}{{{t^{n + 1}} - {t^n}}} = \frac{{{P_W} - {P_\infty }}}{\rho } + \frac{1}{2}{\left| {\nabla {\varphi _j}} \right|^2}\] (Eq. 10)

Where equation (8) represents discretization of the Green’s integral formula and M is the number of cubic splines on the bubble surface and MN is the number of total segments on the bubble surface and the internal wall of the cylindrical vein. Equations (9) and (10) represent discretization of the Eulerian and Lagrangian forms of the unsteady Bernoulli equation, respectively.

The mathematical model of the dynamic response of the compliant coating to its nearby liquid domain behavior has been discretized as:

\[\frac{{\psi _j^{n + 1} - \psi _j^n}}{{{t^{n + 1}} - {t^n}}} = \frac{T}{m}\frac{{\eta _{j + 1}^n - 2\eta _j^n + \eta _{j - 1}^n}}{{{{(\Delta z)}^2}}} - \frac{K}{m}\eta + ({P_W} - {P_\infty })\] (Eq. 11)

Since equation (11) represents an explicit discretized equation, then attempts have been done for obtaining the stability conditions of this scheme. Therefore the time step for the evolution of the problem and ∆z along the internal surface of the cylindrical vein, have been chosen in a manner that has satisfied the stability conditions of this explicit scheme.

The boundary conditions for the numerical solution of equation (11) \[\frac{{\partial \eta }}{{\partial z}} = 0\] are at z=0 and η=0 at z=±z_cs. The former boundary condition is because of the symmetry, while the latter boundary conditions are because of the attachment of the compliant coating to the rigid wall.

Numerical implementation

By having the distribution of the normal velocity on the internal surface of the cylindrical vein and by having the distribution of the velocity potential on the boundary of the gaseous plug, the Green’s integral formula provides the solution for obtaining normal velocity on the boundary of the gaseous plug at any instance as well as velocity potential on the internal surface of the cylindrical vein. In the case of having a rigid cylindrical vein, the normal velocity on its internal surface is known and is equal to zero. Also according to Best ²¹ the velocity potential on the boundary of an explosive growing bubble at the instant of its initial minimum volume is known and is equal to zero as well. Then, as it is stated above, the Green’s integral formula provides the solution for obtaining the velocity potential on the internal surface of the rigid cylindrical vein and the normal velocity on the boundary of the gaseous plug. For the next time step, the velocity potential on the boundary of the gaseous plug can be obtained from the discretized unsteady Bernoulli equation. Furthermore, the normal velocity on the internal surface of the rigid cylindrical vein is known at every time step and is equal to zero.

Therefore by definition of a variable time step as:

\[\Delta t = {\min _i}\left( {\frac{{\Delta \varphi }}{{\frac{{{P_{gp}} + {P_\infty }}}{\rho } + {{\left| {\nabla {\varphi _j}} \right|}^2}}}} \right)\] (Eq. 12)

The numerical simulation of the explosive growth of a gaseous plug inside a rigid cylindrical vein can proceed through its elapsed time from the instant of bubble generation to the instant when the bubble reaches its maximum volume. In equation (12) that comes from the discretized unsteady Bernoulli equation and has been modified by choosing the value of P_gp + P_∞ instead of the pressure difference of P_gp – P_∞, ∆φ is some constant and is equal to the maximum difference of the velocity potential on the boundary of the gaseous plug between two successive time steps.

The discretized forms of the equations which provide the boundary conditions at the liquid domain-compliant surface interface are given as:

\[\left[ {\frac{{\partial \varphi }}{{\partial r}}} \right]_j^n = \left[ {\frac{{\partial \eta }}{{\partial t}}} \right]_j^n = \left[ \psi \right]_n^j\](Eq. 13)

and

\[\left[ {\frac{{\partial \varphi }}{{\partial r}}} \right]_j^n = \left[ {\frac{{\partial \eta }}{{\partial t}}} \right]_j^n = \left[ \psi \right]_j^n\](Eq. 14)

Solution strategy for the case of a rigid cylindrical vein with an internal compliant coating

If the vein has been assumed as a rigid cylinder with an internal compliant coating, the velocity and displacement of any point on the surface of the compliant coating have to be obtained from equation (5) which represents the mathematical model of a thin spring-backed compliant surface. On the other hand equation (5) needs to have the initial distribution of the pressure on the compliant surface at the initial instant of the explosive growth of the gaseous plug for starting of the solution. But the initial pressure distribution on the compliant surface at the initial instant of the explosive growth of the gaseous plug is not known. The initial pressure distribution on the compliant surface has to be found from the acceleration of the surface of the compliant coating at the initial instant of the problem. Furthermore, it should be noted that the acceleration of the surface of the compliant coating can be calculated from equation (5) which needs to have the distribution of the pressure on its surface. Therefore, the unknown physical parameters of the problem are one more than the known physical parameters. Consequently, the numerical simulation of the problem under investigation in the case of the explosive growth of a gaseous plug inside a rigid cylinder with an internal compliant coating is not possible unless an iteration scheme has been employed for obtaining the initial distribution of the pressure on the surface of the compliant coating.

In this study, an iteration scheme which has been proposed by Duncan and Zhang²² for a collapsing cavitation bubble from its maximum volume and adapted by Shervani-Tabar²³ for an explosion bubble growing from its initial minimum volume has been employed for the numerical simulation of interaction of the dynamic behavior of the liquid domain around the gaseous plug and the dynamic reaction of the internal compliant coating of the cylindrical vein.

Non-dimensionalizing of the numerical solution

The numerical solution of the problem has been non-dimensionalized by R_m, the maximum radius of the gaseous plug, as a scale for any length, \[{R_m}{\left( {\frac{\rho }{{{P_{gp}} - {P_\infty }}}} \right)^{\frac{1}{2}}}\] as a scale for time, \[{R_m}{\left( {\frac{{{P_{gp}} - {P_\infty }}}{\rho }} \right)^{\frac{1}{2}}}\] as a scale for velocity potential and \[{\left( {\frac{{{P_{gp}} - {P_\infty }}}{\rho }} \right)^{\frac{1}{2}}}\]as a scale for velocity.

The non-dimensional parameters which identify the elastic characteristics of the compliant coatings are defined as:

\[M* = \frac{m}{{\rho {R_m}}}\]

\[K* = \frac{{K{R_m}}}{{{P_{gp}} - {P_\infty }}}\]