Fluid structure interaction for small deformations in Hookean vessel
The governing equations for the Hookean vessel
The averaged Cauchy equations
The Cauchy equations were derived in the section ref{sec:CauchyEqnMotion} in cylindrical coordinates (see equation \eqref{eq:185}. We will in this section average the equations for the z-direction and r-direction over the vessel wall, assuming azimuthal symmetry (i.e., all \( \partial(\cdot)/\partial \theta-\mathrm{terms}=0 \)). Based on this assumption the Cauchy equation \eqref{eq:190} in the \( z \)-direction may be integrated over the vessel wall to yield: $$ \begin{equation} \tag{144} \int_{A_w} \rho_w a_z \, \dA = \int_{A_w} \left ( \frac{1}{r} \partd{}{r} \left ( r \tau_{zr} \right ) + % \partd{\sigma_z}{z} + \rho_w b_z \right ) \;\dA \end{equation} $$ where \( \rho_w \) denotes the density of the vessel wall, which is assumed to be constant the following derivation.
By introducing the common notation that a bar-ed quantity denote the cross-sectional averaged of the same quantity without a bar, i.e., \( \bar{(\cdot)}= \int_A (\cdot) \dA/A \), the first term of equation (144) may be reformulated to: $$ \begin{equation} \tag{145} \int_{A_w} \rho_w a_z \dA = \rho_w \bar{a}_z A_w \approx \rho_w \bar{a}_z 2 \pi r_i h \end{equation} $$ where \( A_w \) is the vessel wall area, i.e., \( A_w = \pi (r_o^2 - r_i^2) \), and \( r_i \) and \( r_o \) denote inner and outer radius, respectively and \( h \) the wall thickness. For a thin walled vessels \( h/r_i \ll 1 \): $$ \begin{equation} \tag{146} A_w = \pi (r_o^2 - r_i^2) = \pi ( 2r_i h + h^2) \approx 2 \pi r_i h \end{equation} $$ and thus the approximation in equation (145) is valid for thin walled structures.
By expanding \( dA=r \, dr d\theta \), the first term in the rhs of equation (144) take the form: $$ \begin{equation} \tag{147} \Int_0^{2\pi} \Int_{r_i}^{r_o} \frac{1}{r} \partd{}{r} \left ( r \tau_{zr} \right ) \; r \, dr d\theta % = \left . 2 \pi r \tau_{rz} \right |^{r_o}_{r_i} = 2\pi r_i \tau_w \end{equation} $$ where we have assumed \( \tau_{rz} |_{r=r_o} = 0 \) and \( \tau_{rz} |_{r=r_i} = -\tau_w \).
The second term on the rhs of equation (144) may similarly be written: $$ \begin{equation} \tag{148} \int \partd{\sigma_z}{z} \dA = \overline{\partd{\sigma_z}{z}} \, A_w \approx \overline{\partd{\sigma_z}{z}} \, 2\pi r_i h \end{equation} $$
The last term of equation (144) may be treated in exactly the same manner as the lhs. Thus, by substitution of equations (145), (147), and (148) into equation (144) we get: $$ \begin{equation} \tag{149} \rho_w \bar{a}_z A_w = 2\pi r_i \tau_w + \overline{\partd{\sigma_z}{z}} \, A_w + \rho_w \bar{b}_z A_w \end{equation} $$ which by neglection of body forces and assumption of a thin walled structure reduces to: $$ \begin{equation} \tag{150} \rho_w \bar{a}_z \,h = \overline{\partd{\sigma_z}{z}} \, h + \tau_w \end{equation} $$ Thus, equations (149) and (150) represent cross-sectional averaged Cauchy equations in the axial direction.
In order to derive the cross-sectional averaged Cauchy equations in the radial, we proceed in the same way as above, for the Cauchy equation \eqref{eq:188} in the radial direction: $$ \begin{equation} \tag{151} \int \rho a_r \, \dA = \int \left ( \partd{\sigma_r}{r} + \frac{\sigma_r - \sigma_{\theta}}{r} % + \partd{\tau_{rz}}{z} + \rho b_r \right ) \; \dA \end{equation} $$ For the first and last terms of equation (151) we proceed in the same manner as for the axial Cauchy equation, whereas for the first term of the rhs we get: $$ \begin{equation} \tag{152} \int \partd{\sigma_r}{r} \; r \, dr \, d\theta = \left . 2\pi \sigma_r r \right |^{r_o}_{r_i} - \int \sigma_r \, dr \, d\theta \end{equation} $$ To evaluate the expression we need boundary conditions and assume \( \sigma_r|_{r=r_i} = -p \) and \( \sigma_r|_{r=r_o} = 0 \). Further, we introduce the hat-symbol for another averaging procedure: $$ \begin{equation} \tag{153} \hat{\sigma}_r = \frac{1}{2\pi h} = \Int_0^{2\pi} \Int_{r_i}^{r_o} \sigma_r \; dr \, d\theta \end{equation} $$ which by substitution into equation (152) yields: $$ \begin{equation} \tag{154} \int \partd{\sigma_r}{r} \; r \, dr \, d\theta = 2\pi p r_i - \hat{\sigma}_r \, 2 \pi h \end{equation} $$ By using the hat-convention for the second term of equation (151) we get: $$ \begin{equation} \tag{155} \int \frac{\sigma_r - \sigma_{\theta}}{r} \; r \, dr \, d\theta = \hat{\sigma}_r \; 2\pi h - \hat{\sigma}_{\theta} \; 2\pi h \end{equation} $$ while the third term of equation (151) may be expressed as: $$ \begin{equation} \tag{156} \int \partd{\tau_{rz}}{z} \; \dA = \overline{\partd{\tau_{rz}}{z}} \; A_w \approx \overline{\partd{\tau_{rz}}{z}} \; 2 \pi r_i h \end{equation} $$ Finally, we may substitute equations (154), (155), and (156) into (151) to get: $$ \begin{equation} \tag{157} \rho_w \bar{a}_r A_w = 2\pi r_i p - \hat{\sigma}_{\theta} \, 2 \pi h + \overline{\partd{\tau_{rz}}{z}} \; A_w + \rho_w \bar{b}_r A_w \end{equation} $$ which represents a cross-sectionally averaged Cauchy equation in the radial direction. For a thin-walled structure without body forces, equation (157) reduce to: $$ \begin{equation} \tag{158} \rho_w \bar{a}_r h = p - \frac{\hat{\sigma}_{\theta}\, h}{r_i} + \overline{\partd{\tau_{rz}}{z}} \, h \end{equation} $$
The constitutive equation for plane stress
In general the constitutive equations for a Hookean material in a plane stress situation have been provided previously in equation \eqref{eq:255}. For cylindrical coordinates \( (r,\theta, z) \) we have also derived previously equation (save for some notation) that the circumferential strain \( E_{\theta \theta} = u_r/r \), when \( u_r \) denote radial displacement. Thus, in engineering notation the constitutive equation for a thin walled vessel of a Hookean material, takes the form in cylindrical coordinates: $$ \begin{equation} \tag{159} \sigma_z = \frac{\eta}{1-\nu_p^2} \left ( \partd{u_z}{z} + \nu_p \frac{u_r}{r} \right ) \quad \text{and} \quad % \sigma_\theta = \frac{\eta}{1-\nu_p^2} \left ( \frac{u_r}{r} + \nu_p \partd{u_z}{z} \right ) \end{equation} $$ In the following we drop the symbols for averaging, and substitute equation (159) in the averaged Cauchy equations ((150) and (158)): $$ \begin{align} \rho_w a_z \,h &= \frac{\eta \, h}{1-\nu_p^2} \left ( \partdd{u_z}{z} + \frac{\nu_p}{r} \partd{u_r}{z} \right )+ \tau_w \\ \rho_w a_r h & = p - \frac{\eta \, h}{(1-\nu_p^2) \, r_i} \left ( \frac{u_r}{r} + \nu_p \partd{u_z}{z} \right ) \end{align} $$ while assuming equation (159) to be a valid constitutive equation also for averaged stress/strain relations (the term \( \overline{\partd{\tau_{rz}}{z}} \) was discarded as ???). Assume further, that \( a_z = \partial^2 u_z/\partial t^2 \) and \( a_r = \partial^2 u_r/\partial t^2 \), i.e., neglect convective terms, to obtain the averaged governing equations for the thin walled vessel: $$ \begin{align} \tag{160} \partdd{u_z}{t} &= \frac{\eta }{(1-\nu_p^2) \rho_w} \left ( \partdd{u_z}{z} + \frac{\nu_p}{r} \partd{u_r}{z} \right )+ \frac{\tau_w}{\rho_w \,h} \\ \partdd{u_r}{t} & = \frac{p}{\rho_w \, h} - \frac{\eta }{(1-\nu_p^2) \, \rho_w } \left ( \frac{u_r}{r^2} + \frac{\nu_p}{r} \partd{u_z}{z} \right ) \tag{161} \end{align} $$ where we implicitly assume that \( r\approx r_i \).
The governing equations for the fluid
The momentum equations without convective terms: $$ \begin{align} \tag{162} \partd{v_z}{t} & = -\frac{1}{\rho} \partd{p}{z} + \nu \left (\partdd{v_z}{r} + \frac{1}{r} \partd{v_z}{r} \right )\\ \tag{163} \partd{v_r}{t} & = -\frac{1}{\rho} \partd{p}{r} + \nu \left ( \frac{1}{r} \partd{}{r} \left (r \partd{v_r}{r} \right ) \right ) \end{align} $$
The equation for conservation of mass takes the form in cylindrical coordinates: $$ \begin{equation} \tag{164} \partd{v_r}{r} + \frac{v_r}{r} + \partd{v_z}{z} = 0 \end{equation} $$
The momentum equations
Assume solutions on the form: $$ \begin{equation} \tag{165} p = \hat{p} e^{i \omega \; (t-z/c)}, \quad% v_z = \hat{v}_z e^{i \omega \; (t-z/c)}, \quad % v_r = \hat{v}_r e^{i \omega \; (t-z/c)} \end{equation} $$ which by substitution into equation (162) yields: $$ \begin{align} \tag{166} i \omega \; \hat{v}_z= -\frac{1}{\rho} \left ( \frac{-i \omega}{c} \right ) \; \hat{p}+ % \frac{\nu}{r_i^2} \left (\partdd{\hat{v}_z}{y} + \frac{1}{y} \partd{\hat{v}_z}{y} \right ) \end{align} $$ where we have introduced the nondimensional scale \( y = r/r_i \) and the Womersley parameter, previously introduced in equation \eqref{eq:147}, $$ \begin{align} \alpha = r_i \sqrt{\frac{\omega}{\nu}} \end{align} $$ Rearrangement of equation (166) yields: $$ \begin{align} \tag{167} \partdd{\hat{v}_z}{y} + \frac{1}{y} \partd{{v}_z}{y} + i^3 \alpha^2 \hat{v}_z = \frac{i^3 \alpha^2}{\rho c}\; \hat{p} \end{align} $$ Now, by introducing another scale: $$ \begin{align} \tag{168} s = i^{3/2} \alpha y = i^{3/2} \alpha r/r_i \end{align} $$ equation (167) may be transformed to an inhomogeneous Bessel equation of order zero (see equation \eqref{eq:211}): $$ \begin{align} \tag{169} \partdd{\hat{v}_z}{s} + \frac{1}{s} \partd{{v}_z}{s} + \hat{v}_z = \frac{1}{\rho c}\; \hat{p} \end{align} $$ By proceeding in the same manner for equation (163) we get: $$ \begin{align} i \omega \hat{v}_r & = -\frac{1}{\rho r_i} \partd{\hat{p}}{y} + \frac{\nu}{r_i^2} \left ( \partdd{\hat{v}_r}{y} + \frac{1}{y} \partd{\hat{v}_r}{y} % - \frac{\hat{v}_r}{y^2} \right ) \end{align} $$ which by rearrangement may be presented: $$ \begin{align} \tag{170} \partdd{\hat{v}_r}{y} + \frac{1}{y} \partd{\hat{v}_r}{y} + i^3 \alpha^2 \hat{v}_r - \frac{\hat{v}_r}{y^2} = \frac{r_i}{\mu} \partd{\hat{p}}{y} \end{align} $$ which also may be transformed to a Bessel equation by introducing the scale in equation (168), albeit of order one: $$ \begin{align} \tag{171} \partdd{\hat{v}_r}{s} + \frac{1}{s} \partd{\hat{v}_r}{s} + \left (1-\frac{1}{s^2} \right ) \hat{v}_r = \frac{i r_i}{\mu \alpha^2} \partd{\hat{p}}{y} = \frac{r_i}{\mu i^{3/2} \alpha} \partd{\hat{p}}{s} \end{align} $$ For the continuity equation (164) we get: $$ \begin{align} \tag{172} \partd{\hat{v}_r}{r} + \frac{1}{r} \hat{v}_r - \frac{i \omega}{c} \hat{v}_z = 0 \end{align} $$ which by introduction of the scale in equation (168) is transformed to: $$ \begin{align} \tag{173} \frac{1}{s} \partd{}{s} \left (s \hat{v}_r \right ) = - \frac{i^{3/2} \sqrt{\nu \omega}}{c} \hat{v}_z \end{align} $$ or alternatively: $$ \begin{align} \tag{174} \frac{1}{y} \partd{}{y} \left (y \hat{v}_r \right ) = - \frac{i r_i \omega}{c} \hat{v}_z \end{align} $$ The solutions of homogeneous Bessel equations of order zero and one, like equations (169) and (171), are given by their corresponding Bessel functions of order zero and one. Thus for general solutions of equations (169) and (171) we must provide particular solutions. Assume that a particular solution \( \hat{v}_z^p \) of equation (169) is given by: $$ \begin{equation} \tag{175} \hat{v}_z^p = B_1 J_o(ks) \quad \text{with} \quad \hat{p} = A J_0 (ks) \end{equation} $$ where \( k \) is to be determined. Substitution into equation (169) gives: $$ \begin{equation} \tag{176} \partdd{\hat{v}_z}{s} + \frac{1}{s} \partd{{v}_z}{s} + \hat{v}_z = \frac{1}{\rho c}\; \hat{p} % B_1 k^2 \frac{d^2}{dt^2} J_0(t) + \frac{B_1 k}{s} \frac{d}{dt} J_0(t) + B_1 J_0(ks) = \frac{A}{\rho c} J_0(ks) \end{equation} $$ where \( t=ks \). From the useful properties of the Bessel functions given in \eqref{eq:223} and \eqref{eq:226} we may deduce: $$ \begin{align} \frac{d}{dt} J_0(t) &= - J_1(t) \tag{177} \\ \frac{d}{dt} \left (t J_1(t) \right ) & = t J_0(t) = J_1(t) + t \frac{d}{dt} J_1(t) \tag{178}\\ \frac{d^2}{dt^2} J_0(t) & = -\frac{d}{dt} J_1(t) = \frac{J_1(t)}{t} - J_0(t) \tag{179} \end{align} $$ Equations (177) and (179) may subsequently be substituted into equation (176) to yield: $$ \begin{equation} \begin{aligned} B_1 k^2 \left ( \frac{J_1(ks)}{ks} - J_0(ks) \right ) & - \frac{B_1 k}{s} J_1(ks) + B_1J_0(ks) \\ & = B_1 (1-k^2) J_0(ks) = \frac{A}{\rho c} J_0(ks) \end{aligned} \end{equation} $$ Thus, equation \( \hat{v}_z^p \) in (175) is a valid particular solution provided that: $$ \begin{align} \tag{180} B_1 = \frac{A}{(1-k^2) \rho c} \end{align} $$ and a general solution of equation (169) is given by: $$ \begin{align} \tag{181} \hat{v}_z &= \hat{v}_z^h + \hat{v}_z^p = \frac{C_1}{J_0(i^{3/2} \alpha)} \; J_0(s) + \frac{A}{ \rho c \, (1-k^2)} \; J_0(ks) \end{align} $$ where \( C_1 \) and \( k \) are constants to be determined (\( C_1 \) is scaled with \( J_0(i^{3/2} \alpha) \) for analogy with rigid pipe solution).
Similarly, we propose a particular solution for equation (171) of the form: $$ \begin{align} \tag{182} \hat{v}_r^p = B_2 J_1 (ks) \end{align} $$ and we deduce from equation (177) and (179): $$ \begin{align} \tag{183} \frac{d^2}{dt^2} J_1 = \frac{d}{dt} J_0 - \frac{d}{dt} \left (\frac{J_1}{t} \right ) = % \left ( \frac{1}{t^2} -1 \right ) J_ 1 -\frac{1}{t} \frac{d}{dt} J_1 \end{align} $$ Substitution of equations (182) and (183) into equation (171) yields: $$ \begin{equation} \tag{184} \begin{aligned} B_2 k^2 \left ( \frac{1}{(sk)^2} -1 \right ) J_ 1 & -\frac{B_2 k^2}{sk} \frac{d J_1}{dt} +\frac{B_2 k}{s} \frac{d J_1}{dt} + \left (1 - \frac{1}{s^2} \right ) B_2 J_1 \\ &= B_2 \; (1-k^2) J_1 (ks) = \frac{r_i A}{\mu i^{3/2} \alpha} \partd{\hat{p}}{s}= \frac{- r_i A k}{\mu i^{3/2} \alpha} \;J_1(ks) \end{aligned} \end{equation} $$ where we have used: $$ \begin{equation} \tag{185} \partd{\hat{p}}{s} = -Ak \; J_1(ks) \end{equation} $$ From equation (184) we get that \( \hat{v}_r^p \) in equation (182) is a particular provided that: $$ \begin{equation} B_2 = \frac{r_i A k }{\mu i^{3/2} \alpha (k^2 -1)} \end{equation} $$ And consequently, a general solution of equation (171) is: $$ \begin{equation} \tag{186} \hat{v}_r = \hat{v}_r^h + \hat{v}_r^p = \frac{C_2}{J_0(i^{3/2} \alpha)} \; J_1(s) + \frac{r_i A }{\mu i^{3/2} \alpha}\, \frac{k}{k^2-1} \; J_1(ks) \end{equation} $$
Fulfillment of the continuity equation
The two general solutions given by equations (181) and (186) must also satisfy the continuity equation (173). Observe first that: $$ \begin{equation} \tag{187} \begin{aligned} \frac{1}{s} \left ( \frac{d}{ds} ( \hat{v}_r s ) \right ) & = \frac{1}{s} \left ( \frac{d}{ds} ( \hat{v}_r^h s ) + \frac{d}{ds} ( \hat{v}_r^p s ) \right ) \\ & = \frac{1}{s} \frac{d}{ds} ( \hat{v}_r^h s ) + \frac{1}{s} \frac{d}{dt} ( \hat{v}_r^p(t) t ) \end{aligned} \end{equation} $$ From equation (186) we get: $$ \begin{align} \tag{188} \frac{1}{s} \frac{d}{ds} ( \hat{v}_r^h s ) = \frac{C_2}{J_0(i^{3/2} \alpha)} J_0 (s) \end{align} $$ by the use of equation (178) and similarly: $$ \begin{align} \tag{189} \hat{v}_r^p (t) t = \frac{r_i A k}{\mu i^{3/2} \alpha (k^2-1)} \; J_1(t) t \end{align} $$ and thus: $$ \begin{equation} \tag{190} \begin{aligned} \frac{1}{s} \frac{d}{dt} ( \hat{v}_r^p(t) t ) & = \frac{r_i A k}{\mu i^{3/2} \alpha (k^2-1)} \; \frac{sk}{s} J_0(sk) \\ & = \frac{r_i A k^2}{\mu i^{3/2} \alpha (k^2-1)} J_0(sk) \end{aligned} \end{equation} $$ Substitution of equation (188) and (190) in equation (186) results in the following expression: $$ \begin{equation} \tag{191} \frac{1}{s} \left ( \frac{d}{ds} ( \hat{v}_r s ) \right ) = \frac{C_2}{J_0(i^{3/2} \alpha)} J_0 (s) + \frac{r_i A k^2}{\mu i^{3/2} \alpha (k^2-1)} J_0(ks) \end{equation} $$ which is the left hand side of the continuity equation (173). For the right hand side of the continuity equation we get by substitution of equation (181): $$ \begin{equation} \tag{192} \begin{aligned} - \frac{i^{3/2} \sqrt{\nu \omega}}{c} \hat{v}_z & = \frac{i r_i \omega}{i^{3/2} \alpha c} \hat{v}_z \\ & = \frac{i r_i \omega}{i^{3/2} \alpha c J_0(i^{3/2} \alpha)} \, C_1 \; J_0(s) + \frac{i r_i \omega}{i^{3/2} \alpha c} \; \frac{A}{ \rho c \, (1-k^2)} \; J_0(ks) \end{aligned} \end{equation} $$ By comparing equation (191) and equation (192) we find that the ratio of \( C_2 \) and \( C_1 \) must satisfy: $$ \begin{align} \tag{193} \frac{C_2}{C_1} = \frac{i r_i \omega}{i^{3/2} \alpha c} \end{align} $$ and \( k \) must fulfill the condition: $$ \begin{align} \tag{194} \frac{k^2}{\mu} = \frac{i^3 \omega}{\rho c^2} \quad \Rightarrow \quad k=\pm \frac{i^{3/2} \sqrt{\nu \omega}}{c} = \pm \frac{i^{3/2} r_i \omega}{\alpha c} \end{align} $$ Thus, one may argue that for physiological values of \( \nu \), \( \omega \), and \( c \) we will have \( |k|\ll 1 \). Further, by introducing the definition of the scale in equation (168) we get: $$ \begin{align} \tag{195} ks = \pm \frac{i^{3/2} r_i \omega}{\alpha c} i^{3/2} \alpha y = \mp \frac{i r_i \omega}{c} y =\mp \frac{i\omega}{c} r \end{align} $$ For terms involving \( J_0(ks) \) and \( J_1(ks) \) one may approximate: $$ \begin{align} \tag{196} J_0(ks) = I_0(\mp \frac{r_i \omega}{c} \, y) \approx 1 \quad \text{and} \quad J_1(ks) = I_1(\mp \frac{r_i \omega}{c} \, y) \approx \mp \frac{i r_i \omega}{2c} \, y \end{align} $$ where \( I_n \) denotes the modified Bessel function of order \( n \) (see equation \eqref{eq:243} in the section ref{sec:prop-bess-funct}). From equation (194) and (196) we get: $$ \begin{align} \tag{197} k J_1(ks) &\approx \pm \frac{i^{3/2} r_i \omega}{\alpha c} \mp \frac{i r_i \omega}{2c} \, y = -\frac{i^{3/2} r_i \omega}{\alpha c} \frac{i r_i \omega}{2c} \, y \end{align} $$
This is a simple second order algebraic equation in \( x \) which is easy to solve. The complex solutions of equation \eqref{eq:279} will be functions of \( \nu_p \), \( k \), and \( \alpha \), i.e., \( x=x(\nu_pk,\alpha) \), which may readily be found with e.g., Maple. From equation \eqref{eq:275} one may deduce: $$ \begin{equation} \tag{198} x = k \frac{\eta}{1-\nu_p^2} \frac{1}{\rho_w} \frac{1}{c^2} = \frac{h}{\rho r_i} \frac{\eta}{1-\nu_p^2} \frac{1}{c^2} \end{equation} $$
The Moens-Korteweg formula for the wave speed \( c_0 \) under inviscid conditions is given by equation (94), which by substitution and rearrangement in equation (198) results in: $$ \begin{equation} \tag{199} \gamma = \frac{c}{c_0} = \left (\frac{(1-\nu_p^2) \; x}{2} \right )^{-1/2} \end{equation} $$ The symbol \( \gamma \) has been introduced for the ratio of the wave speed \( c \) for the coupled problem to \( c_0 \). As \( x \) is complex, so will \( \gamma \) be, and we write \( \gamma = \gamma_r + i \gamma_i \) for the real and complex parts. Thus, having found the solutions of equation \eqref{eq:279}, we may plot \( \gamma_r \) as a function of \( \nu_p \), \( k \), and \( \alpha \) (see figure ref{fig:c_c0}).
The dependence of \( k \) and \( \nu_p \) may also be illustrated in 3D, albeit some less quantitative accuracy (see figure ref{fig:c_c03d}).
Figure 10: \( \gamma_r = \Re(c/c_0) \) as function of \( \alpha \) and \( \nu_p \) for \( k = 0.1 \) (left) and \( k \) for \( \nu_p = 0.25 \) (right).