Coupled Solid–Fluid Problems

Abstract

We considered a variety of problems in Chaps. 2–10 that fall within the domain of either biosolid mechanics or biofluid mechanics, each of which is very important in its own right. Whether in the body (in vivo) or in the laboratory (in vitro), however, many “real-life” problems simultaneously involve solid–fluid interactions. For example, although we may seek to determine the stresses in the limbs of a pilot who has ejected from an aircraft, for purposes of identifying safety measures, it is the wind that induces the applied loads of importance; aneurysms may be considered as thin-walled, nearly spherical membranes that exhibit a solidlike character, but the applied loads are due to the internal flowing blood and the surrounding cerebrospinal fluid; mechanotransduction in bone, which exhibits a strong solidlike behavior, appears to be influenced directly by both loads due to weight bearing and those due to the flow of blood and bone fluid within the many different canals within the bone; and an atomic force microscopic examination of the mechanics of a cell may primarily reveal the properties of the cortical membrane and underlying solidlike cytoskeleton, but flow of the cytosol likely plays a key role as well. Hence, from these simple examples and many more like them, we see that solid–fluid interactions are important at the organism, organ, tissue, cellular, and molecular levels. Indeed, although it tends to be convenient to introduce students to a subject by focusing only on that subject, most research and clinical problems require interdisciplinary and multidisciplinary approaches (i.e., analysis and design of coupled problems). Such problems are typically complex and require advanced approaches, but here we consider a few introductory examples.

Appendix 11: Wave Equations

Recall that we studied the steady flow of a Newtonian fluid within a rigid circular cylinder in Sect. 9.2. Notwithstanding the unsteadiness of flow in actual arteries and airways, results for steady flows are often used to estimate the mean wall shear stress within arteries and airways as well as in the more appropriate applications of flows in veins, gravity-fed IV tubes, and so forth. Inasmuch as the unsteadiness can be important, we also considered a special case of pulsatile flow within a rigid circular cylinder in Sect. 9.5, which provided more insight into potential effects of pulsatility on wall shear stress. Yet, as noted eloquently by Zamir (2000), pulsatile flows differ fundamentally within rigid tubes versus either elastic or viscoelastic tubes (e.g., arteries, airways, and so forth). No matter how small the distensibility of the tube, flow propagates within a distensible tube as a wave having a finite wave speed. Hence, in the presence of any downstream material or geometric discontinuity (e.g., a bifurcation or a branch), the wave(s) can reflect, which can lead to forward and backward traveling waves that add either constructively or destructively. In the case of constructive interference, for example, a local systolic blood pressure due to a forward traveling wave can be augmented by a backward traveling wave, which thereby increases the local pulse pressure. Of course, an increased pulse pressure in the ascending aorta can increase the work-load on the heart whereas an increased pulse pressure within any segment of the arterial tree can also be a strong mechanobiological stimulus for cells within the wall of the artery.

In contrast, waves do not exist in a rigid tube. The so-called wave speed is infinite in a rigid tube (recall the Moens-Korteweg equation in Observation 7.2), which means that the pressure and flow are transmitted instantaneously along the length of the tube, hence resulting in an overall bulk motion. There is, therefore, strong motivation to study unsteady flows in distensible tubes, that is, fluid-solid-interactions. Indeed, as we also noted in Observation 7.2, “pulse wave velocity” (PWV) within large arteries is now recognized as an important indicator or initiator of diverse cardiovascular diseases and it merits increased scientific consideration. Although the complexities of analytical or computational solutions of wave motion within a distensible vasculature are beyond the scope of an introductory textbook, let us note some basics with regard to traveling waves.

A number of equations in mathematical physics exhibit common features and can be classified as elliptic, parabolic, or hyperbolic (partial) differential equations. They are,

$$ \begin{array}{l} Elliptic\kern2.64em 0={\nabla}^2\phi, \\ {} Parabolic\kern1.44em \frac{\partial \phi }{\partial t}={\alpha}^2{\nabla}^2\phi, \\ {} Hyperbolic\kern0.84em \frac{\partial^2\phi }{\partial {t}^2}={c}^2{\nabla}^2\phi, \end{array} $$

(A.11.1)

where α and c are parameters, ${\nabla}^2 \equiv \frac{\partial^2}{\partial {x}^2}+\frac{\partial^2}{\partial {y}^2}+\frac{\partial^2}{\partial {z}^2}$ in three spatial dimensions (3-D) in Cartesian coordinates, and t is time. These equations are known, respectively, as the Laplace equation, the diffusion equation, and the wave equation; general methods of solution can be found in textbooks on partial differential equations or applied mathematics.

Let us consider briefly the wave equation in one spatial dimension (1-D), namely

$$ \frac{\partial^2\phi }{\partial {t}^2}={c}^2\frac{\partial^2\phi }{\partial {x}^2}, $$

(A.11.2)

where c is a constant and we seek solutions of the form $\phi =\phi \left(x,t\right)$, where x is the spatial coordinate location and t is time. It was recognized by D’Alembert that functions of the form $\phi =\phi \left(x-ct\right)$ and $\phi =\phi \left(x+ct\right)$ satisfy this linear partial differential equation, hence these solutions can be superimposed to yield the general solution, $\phi =\phi \left(x-ct\right)+\phi \left(x+ct\right)$, which suggests that a wave could begin at x = 0 and travel simultaneously in both the positive and negative direction at speed c. That these functions satisfy the wave equation can be appreciated easily by noting that, if $\phi =\widehat{\phi}(u)$, with $u=x-ct$, then

$$ \begin{array}{l}\frac{\partial \phi }{\partial t}=\frac{\partial \widehat{\phi}}{\partial u}\frac{\partial u}{\partial t}=\frac{\partial \widehat{\phi}}{\partial u}\left(-c\right)\kern0.36em \mathrm{and}\kern0.36em \frac{\partial^2\phi }{\partial {t}^2}=\frac{\partial^2\widehat{\phi}}{\partial {u}^2}\left(-c\right)\left(-c\right)\\ {}\frac{\partial \phi }{\partial x}=\frac{\partial \widehat{\phi}}{\partial u}\frac{\partial u}{\partial x}=\frac{\partial \widehat{\phi}}{\partial u}(1)\kern0.36em \mathrm{and}\kern0.36em \frac{\partial^2\phi }{\partial {x}^2}=\frac{\partial^2\widehat{\phi}}{\partial {u}^2}\end{array} $$

(A.11.3)

whereby

$$ {c}^2\frac{\partial^2\widehat{\phi}}{\partial {u}^2} \equiv {c}^2\frac{\partial^2\widehat{\phi}}{\partial {u}^2}. $$

(A.11.4)

Moreover, that the parameter c has units of speed can be seen easily via the unit equation for A.11.2, namely

$$ \frac{\left[\phi \right]}{{\left[T\right]}^2}={\left[c\right]}^2\frac{\left[\phi \right]}{{\left[L\right]}^2}\kern0.36em \to \kern0.36em {\left[c\right]}^2=\frac{{\left[L\right]}^2}{{\left[T\right]}^2}\kern0.36em \to \left[c\right]=\frac{\left[L\right]}{\left[T\right]}. $$

(A.11.5)

Notwithstanding the interpretive advantage of D’Alembert’s solution, one often pursues solutions to the wave equation using the method of separation of variables. That is, we assume that the solution can be written in the form $\phi \left(x,t\right)=X(x)T(t)$, whereby we have

$$ X\frac{d^2T}{d{t}^2}={c}^2T\frac{d^2X}{d{x}^2}\kern0.36em \to \kern0.36em \frac{1}{T}\frac{d^2T}{d{t}^2}=\frac{c^2}{X}\frac{d^2X}{d{x}^2}. $$

(A.11.6)

Of course, the only way for a function of time (left hand side) to equal a function of position (right hand side) for all (x,t) is for both functions to equal the same constant, say A. Hence, we have two second order ordinary differential equations to solve, namely

$$ \frac{d^2T}{d{t}^2}=AT,\kern0.84em \frac{d^2X}{d{x}^2}=\frac{A}{c^2}X. $$

(A.11.7)

Three solutions are possible depending on whether A < 0, A = 0, or A > 0. It can be shown that A > 0 and A = 0 do not lead to periodic solutions and thus are not realistic. If A < 0, however, one finds solutions for both X(x) and T(t). In this case, it is convenient to let $A \equiv - {\lambda}^2\kern0.36em \left(\lambda >0\right)$, whereby the final solution is (Wylie and Barrett 1982)

$$ \phi \left(x,t\right)=X(x)T(t)=\left({C}_1 \cos \frac{\lambda }{c}x+{C}_2 \sin \frac{\lambda }{c}x\right)\left({C}_3 \cos \lambda t+{C}_4 \sin \lambda t\right) $$

(A.11.8)

where C ₁, C ₂, C ₃, C ₄ are material parameters to be determined from boundary and initial conditions. Clearly, this solution is periodic for it repeats itself every time t increases by 2π/λ. As we close this chapter, similar to Chap. 8, we are reminded that many fundamental approaches of applied mathematics prove essential in solving problems in biomechanics, hence one must commit to studying both.

Exercises

• 11.1 Repeat the nondimensional analysis of Sect. 11.2 using as length, time, and mass scales:

  $$ {L}_s=A,\kern1em {T}_s=\frac{A^2H}{Q},\kern1em {M}_s=\rho {A}^2H. $$

  Compare the results with those in Sect. 11.2.

• 11.2 Repeat Exercise 11.1 using

  $$ {L}_s=A,\kern1em {T}_s=\sqrt{\frac{\rho {A}^2}{\Delta p},}\kern1em {M}_s=\rho {A}^2H. $$

• 11.3 Recall from Observation 10.3 that the Buckingham Pi method can be used to nondimensionalize known equations. Show that the governing differential equation of motion for the aneurysm [Eq. (11.53)] can be written in nondimensional form as

  $$ \left(\frac{1}{x^2}+bx\right)\ddot{x}+\frac{3}{2}b{\dot{x}}^2+4m\frac{\dot{x}}{x}+2\frac{f(x)}{x}=F\left(\tau \right), $$

  where

  $$ \begin{array}{ll}x \equiv \lambda, \hfill & f=\frac{T}{c},\hfill \\ {}b=\frac{\rho A}{\rho_sH},\hfill & F=\frac{PA}{c},\hfill \\ {}m=\frac{\mu }{\sqrt{\rho_scH}},\hfill & \tau =\frac{t\sqrt{c}}{\sqrt{\rho_s{A}^2H}},\hfill \end{array} $$

  and c is a material parameter having units of force/length.

• 11.4 If the CSF surrounding a spherical aneurysm is assumed to be ideal, then the governing equation reduces to

  $$ \left(\frac{1}{x^2}+bx\right)\ddot{x}+\frac{3}{2}b{\dot{x}}^2+2\frac{f(x)}{x}=F\left(\tau \right). $$

  If F(τ) = F ₀, a constant, show that the equation can be integrated once in time to yield

  $$ \frac{1}{2}{\dot{x}}^2+\frac{1}{2}b{\dot{x}}^2{x}^3+2{\displaystyle \int xf(x)dx}-\frac{1}{3}{F}_0{x}^3=\mathrm{constant}. $$

  This form of the equation can be related to the first law of thermodynamics (e.g., ẋ ²/2 is a nondimensional kinetic energy term, and –F ₀ x ³/3 is a work-type term related to a pressure times volume). Hint: Note that

  $$ \dot{x}\ddot{x}=\frac{d}{dt}\left(\frac{1}{2}{\dot{x}}^2\right),\kern1em {x}^2\ddot{x}=\frac{d}{dt}\left(\frac{1}{3}{x}^3\right),\kern1em dx=\frac{dx}{dt}dt. $$

• 11.5 If F(τ) ≡ F ₀, a constant, if f(x) = f ₀, a constant surface tension, and if there is no external fluid, then the governing equation of motion of a spherical “soap bubble” is (cf. equation in Exercise 11.3)

  $$ \frac{1}{x^2}\ddot{x}+\frac{2{f}_0}{x}={F}_0. $$

  Solve this equation for x(τ) and comment on its interpretation.

• 11.6 In Sect. 11.3, we saw that the pressure field p(r, t) in a radial flow in a spherical domain for an incompressible, Newtonian fluid is independent of the viscosity; that is, we have found a special flow wherein the same pressure field satisfies both the Navier–Stokes and the Euler equations. Because the flow is radial, we can also define a radial streamline s ≡ r, where ds = dr. Show, therefore, that an unsteady Bernoulli equation can be written in the form

  $$ {p}_A+\rho g{z}_A+\frac{1}{2}\rho {v}_A^2={p}_B+\rho g{z}_B+\frac{1}{2}\rho {v}_B^2+\rho {\displaystyle {\int}_A^B\frac{\partial {v}_s}{\partial t}}ds. $$

  Hint: Integrate the appropriate Euler equation along a radial streamline.

• 11.7 Use the unsteady Bernoulli equation in Exercise 11.6, with point A at r = a, point B at r = ∞, and v _s  ≡ v _r  = g(t)/r ² from our mass balance relation, to show that one obtains the same pressure field p(r, t) from Bernoulli as obtained from Navier–Stokes.

• 11.8 Repeat Example 11.3 for η = 0.0625, ζ = −0.1, c = 1.0, α = 0.5, F(t) = 1.0sin t, x(0) = 0, and $\dot{x}(0)=1$. What does the negative value of ζ induce?

• 11.9 Repeat Exercise 11.8 with x(0) = 1 and $\dot{x}(0)=1$.

• 11.10 Viscoelastic characteristics include instantaneous elasticity, creep, stress relaxation, instantaneous recovery, delayed elasticity, permanent set, and hysteresis. Define and discuss each characteristic.

• 11.11 We found in Eq. (11.59) that the stress relaxation σ(t) in a Maxwell model is given by an exponential decay. Compute the rate of change of stress as a function of time. Observe that the relaxation is initially very rapid; indeed, show that only 37 % of the initial stress remains at time t = t _R , the relaxation time.

• 11.12 Similar to the previous exercise, investigate the creep response of a Kelvin–Voigt model. In particular, compute $\dot{\varepsilon}(t)$ and sketch it versus time. Show, too, that only 37 % of the asymptotic strains remains to be realized after t = t _c .

• 11.13 As we saw in Eq. (11.64), the Kelvin–Voigt model does not allow stress relaxation. Intuitively, we realize that a step change in ε from 0 to ε ₀ at time t = 0 can only be accomplished via an infinite stress (if that were possible) because a viscous dashpot cannot otherwise extend instantaneously. Thereafter, the stress in the viscous element drops to zero, for it requires a strain rate to produce a stress, thus the spring sustains the constant extension with a constant stress. In a similar way, qualitatively discuss the creep response of the Kelvin–Voigt model and, in particular, justify why this behavior is sometimes referred to as a delayed elasticity (cf. Findley et al. 1976, p. 56).

• 11.14 Recall again the stress relaxation of a Maxwell model:

  $$ \sigma (t)={\sigma}_0{e}^{-t/{t}_R}, $$

  where t _R is the so-called relaxation time. Because the model is linear, superposition holds (cf. Sect. 5.5 of Chap. 5). Thus, a set of Maxwell models in parallel (Fig. 11.26) and under a constant strain ε ₀ will stress relax according to

  

  Figure 11.26

  $$ \sigma (t)={\varepsilon}_0{\displaystyle \sum_{i=1}^n{E}_i{e}^{-t/{t}_R^i},} $$

  where ${t}_R^i={\mu}_i/{E}_i$ for all elements i = 1, 2,… n. For a continuous distribution of relaxation times, from 0 to ∞, we have

  $$ \sigma (t)={\varepsilon}_0{\displaystyle {\int}_0^{\infty }R\left({t}_R\right){e}^{-t/{t}_R}d{t}_R,} $$

  where R(t _R ) is called the relaxation spectrum—a distribution of relaxation times. Formulate a similar analysis of the creep response of a set of Kelvin–Voigt elements in series whereby

  $$ \varepsilon (t)={\sigma}_0{\displaystyle \sum_{i=1}^n{C}_i\left(1-{e}^{-t/{t}_c^i}\right)}, $$

  or

  $$ \varepsilon (t)={\sigma}_0{\displaystyle {\int}_0^{\infty }C\left({t}_c\right)\left(1-{e}^{-t/{t}_c}\right)d{t}_c}, $$

  where C(t _c ) is called the retardation spectrum.

• 11.15 A Maxwell model in series with a Kelvin–Voigt model is a four-parameter model sometimes called a Burgers model. If the spring stiffness and dashpot viscosity are given by (E ₁, μ ₁) and (E ₂, μ ₂) for the Maxwell and Kelvin–Voigt components, respectively, show that the creep function (or compliance) for the Burgers model is

  $$ J(t)=\frac{1}{E_1}+\frac{t}{\mu_1}+\frac{1}{E_2}\left(1-{e}^{-{E}_2t/{\mu}_2}\right). $$

• 11.16 Assume that a uniaxial member is subjected to a strain of the form ε(t) = ε _A  sin ωt. If the material behaves elastically, one would expect a stress response of the form σ(t) = σ _A  sin ωt (i.e., in phase). For a viscoelastic response, however, one would expect the stress response to be out of phase with the strain. Hence, let

  $$ \sigma (t)={\sigma}_A \sin \left(\omega t+\phi \right)={\sigma}_A\left( \sin \omega t \cos \phi + \cos \omega t \sin \phi \right), $$

  where ϕ is the phase angle. This form suggests that a complex representation may be useful, namely

  $$ \varepsilon (t)={\varepsilon}_A{e}^{i\omega t},\kern1em \sigma (t)={\sigma}_A{e}^{i\left(\omega t+\phi \right)}, $$

  where $i=\sqrt{-1}$. Show that

  $$ \frac{\sigma (t)}{\varepsilon (t)}=\frac{\sigma_A}{\varepsilon_A}{e}^{i\phi }={G}_1+i{G}_2, $$

  where G ₁ and G ₂ are called the storage modulus and the loss modulus, respectively. Show, too, that

  $$ {G}_1=\frac{\sigma_A}{\varepsilon_A} \cos \phi, \kern1em {G}_2=\frac{\sigma_A}{\varepsilon_A} \sin \phi, \kern1em \tan \phi =\frac{G_2}{G_1}. $$

  Note that, for example, G ₁ ~ 10⁹ Pa, G ₂ ~ 10⁷ Pa, and ϕ ~ 0.01 for a typical polymer.

• 11.17 Following up on Exercise 11.16, note that G* = G ₁ + iG ₂ is called the complex modulus. Show that for the sinusoidal straining in Exercise 11.16, the magnitude of G* = σ _A /ε _A . Note, too, that tanϕ is often called the mechanical loss.

• 11.18 If a Standard model consists of a spring in series with a Kelvin–Voigt element, and the spring has a stiffness E = 1 GPa, whereas the Kelvin–Voigt element has stiffness 10 kPa and viscosity 10⁷ P (poise), plot log J ₁, log J ₂, and log(tan ϕ) versus logω ∈ [−8, 8], where

  $$ {J}^{*}={J}_1-i{J}_2=\frac{\varepsilon_0}{\sigma_0}{e}^{-i\phi } $$

  for creep.

• 11.19 Show that for a Maxwell model subjected to an oscillatory motion,

  $$ {G}^{*}=\frac{\mu^2{\omega}^2/E}{1+{\mu}^2{\omega}^2/{E}^2}+i\left(\frac{\mu \omega }{1+{\mu}^2{\omega}^2/{E}^2}\right) $$

  and

  $$ \tan \phi =\frac{G_2}{G_1}=\frac{E}{\omega \mu }. $$

• 11.20 Consistent with the prior exercise, plot tanϕ and |G*| versus log ω ∈ [−4, 4] given values of E = l GPa and μ = 5 × 10⁹ P. Repeat for the same E but with μ = 5 × 10¹⁰ then 5 × 10⁸ P. Discuss the behavior in terms of changes in t _R . Note that these values are reasonable for a polymer that may be used in a biomedical device.

• 11.21 Show that the frictional drag coefficient C _f for the step-slider in Sect. 11.5.2 is correct as given. Hint: The total drag force F _D can be computed by integrating the shear stresses over the bearing surface; that is,

  $$ {F}_D=w{\displaystyle {\int}_0^{L\left(1-b\right)}{\tau}_1d{x}_1+}{\displaystyle {\int}_{L\left(1-b\right)}^{bL}{\tau}_2d{x}_2,} $$

  where, in general,

  $$ \tau \equiv {\left.\mu \frac{\partial {v}_x}{dy}\right|}_{y=0} $$

• 11.22 The so-called Reynolds’ equation governs general flows in hydrodynamic lubrication theory. In one dimension, show that it can be written as

  $$ \frac{d}{dx}\left(\frac{h^3(x)w}{\mu}\left(\frac{dp}{dx}\right)\right)+6U\frac{dh}{dx}=0. $$

  Hint: Use the result for the velocity v _x (y) in Sect. 11.5.2 and compute the volumetric flow rate

  $$ Q=w{\displaystyle {\int}_0^{h(x)}{v}_x(y)} dy. $$

  Exploit the fact that Q is constant with respect to x even though the gap distance is h = h(x).

• 11.23 Although the Moens-Korteweg equation is theoretically inappropriate for use in arterial mechanics, it provides correct order of magnitude results nonetheless. If ρ ~ 1,060 kg/m³ for blood and E ~ 1 MPa and h/a ~ 0.1 for an artery, estimate the speed of a pressure wave and compare to values measured clinically in a human aorta. Discuss in a 2-page report the implications of increased arterial stiffening with aging and associated increases in the so-called pulse wave velocity (PWV).

• 11.24 Because of the complexity of fluid-solid interactions in the vasculature, sophisticated computational methods must be used to study issues such as the effects of wall stiffness on pressure waves. Discuss in a 2-page report recent advances and associated computational findings regarding regional pulse wave velocities in humans and changes therein due to aging. As a start, see the paper by Xiao N. et al. (2013) Computer simulation of blood flow, pressure, and vessel wall dynamics in a full-body-scale three dimensional model of the human vasculature. J Comp Phys 244: 22-40.

• 11.25 Discuss in a 2-page report the utility of “dynamic similarity” in the design of medical devices that are deployed within the vasculature.