Flow Regulation Using Biphasic-FSI Domains

Overview

In biomechanics, we may often need to regulate flow within a conduit using structures that act as valves. Typical examples include the heart chambers or lymphatic vessels. Modeling valves explicitly may require advanced features in a finite element code, such as the immersed boundary method, not currently available in FEBio. However, with the introduction of the biphasic-FSI material in the fluid-FSI solver (starting with FEBio 3.4)^Shim-2021Shim JJ, Maas SA, Weiss JA, Ateshian GA. Finite Element Implementation of Biphasic-Fluid Structure Interactions in FEBio. J Biomech Eng. 2021 Mar 25. doi: 10.1115/1.4050646. Epub ahead of print. PMID: 33764435., one can easily simulate a domain that opens and closes to fluid flow, to simulate such flow regulation. This case study illustrates this type of flow regulation using two representative, and conceptually similar problems.

A biphasic-FSI material represents a generalization of the fluid-FSI material. It models the dynamic flow of a viscous fluid through a deformable, porous and permeable (biphasic) domain. The fluid viscosity determines frictional interactions within the fluid, whereas the hydraulic permeability determines the frictional interactions between the fluid and the solid, similar to the phenomena modeled in the Darcy-Brinkman equation ^{Brinkman-1949}Brinkman, H.C., 1949. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Flow, Turbulence and Combustion, 1(1), pp.27-34.. The behavior of a biphasic-FSI material reduces to that of a fluid-FSI material when the hydraulic permeability is infinite (its inverse is zero), the solid mass density is set to zero, and the solid stiffness is set to a very small, but non-zero value. In practice, the limiting case of infinite hydraulic permeability is reproduced in FEBio by setting the permeability to zero.

Axisymmetric Expanding Chamber

Geometry and Mesh

In this model, a fluid-filled chamber with deformable elastic wall is supplied continuously with fluid through a bottleneck inlet on the negative $Z-$face (Figure 1). The flow out of the chamber is regulated by a biphasic-FSI material located at the end of a bottleneck on the positive $Z-$face. Since the geometry and flow conditions are axisymmetric, only a 3° wedge geometry is modeled. The entire fluid domain is 12 cm long, whereas each bottleneck is 2 cm long. The inner radius of the chamber is 3 cm and that of the bottlenecks is 1.5 cm. The elastic wall thickness is 0.5 cm. The actual model dimensions are in SI units (length units of meter).

The mesh has 131 elements along the axial $Z-$ direction. The fluid chamber has 19 elements along the radial $Y-$direction, with a mesh bias generated using FEBioStudio’s boundary layer tool at the interface with the wall. The wall has 4 elements along the radial direction. The total number of elements was 3013 (2882 hex8 and 131 penta6 elements), producing 32337 equations.

Material Models

The fluid-FSI domain models a Newtonian viscous fluid with a mass density $\rho=1000\,\text{kg}/\text{m}^3$, shear viscosity $\mu=0.004\,\text{Pa}\cdot\text{s}$ and bulk modulus $K=2.2\,\text{GPa}$. The solid material that regularizes the mesh deformation in the fluid-FSI domain is neo-Hookean, with Young’s modulus $E=10\,\text{Pa}$ and Poisson’s ratio $\nu=0$. The wall material is neo-Hookean with $E=100\,\text{kPa}$ and $\nu=0.3$, with mass density $\rho=1000\,\text{kg}/\text{m}^3$. The biphasic-FSI domain that regulates the flow at the outlet has the same properties as the fluid-FSI domain. In addition, its hydraulic permeability model is const-perm-iso, with a hydraulic permeability value that steps from $k=1\,\text{m}^4/\text{N}\cdot\text{s}$ for open flow to $k=10^{-9}\,\text{m}^4/\text{N}\cdot\text{s}$ for closed flow.

Boundary Conditions

A fluid normal velocity load was prescribed at the inlet, using a parabolic profile with a peak velocity of 0.2 m/s. The parabolic distribution was prescribed using the surface_data attribute to reference a corresponding SurfaceData element in the MeshData section. The velocity profile was associated with a load curve that smoothly raised the velocity from time $t=0$ to $t=1\,\text{s}$ using the cosine function $\frac{1}{2}\left(1-\cos{\pi t}\right)$. No-slip conditions were enforced between the fluid-FSI/biphasic-FSI domains and the wall domain. Zero fluid dilatation was prescribed on the outlet face of the biphasic-FSI domain to allow free efflux when the permeability was set to a high value, along with backflow and tangential stabilization (with $\beta=1$). Solid displacements were prescribed on all three domains to enforce axisymmetric conditions for the deformations, also requiring the use of a symmetry plane on the inclined wedge face. Finally, the $Z-$displacement was fixed on the inlet and outlet boundary surfaces.

To regulate the flow at the outlet, a load curve was associated with the biphasic-FSI material’s hydraulic permeability element perm, flipping its value from $1$ to $10^{-9}$ and back, using a STEP interpolation. The low value of $10^{-9}$ was set for time intervals $1.01\le t \le 2$, $3.01\le t \le 4$, and $5.01\le t \le 6$ (Figure 2).

Analysis Settings

For this model, the analysis setting requested $300$ time steps with step size of $0.02\,\text{s}$, for a total analysis time of $6\,\text{s}$. The solver employed the BROYDEN method with 5 maximum reformations and 50 maximum updates, turning off the diverge_reform and reform_each_time_step parameters for computational efficiency. For the time stepper, the maximum step size was set to $0.02\,\text{s}$ and the minimum step size was set to $2\times 10^{-8}\,\text{s}$, while the Cutback parameter was set to aggressive to accommodate the sudden “closing” of the outlet boundary. A maximum of 20 retries was allowed, and the optimal number of iterations was set to 50.

Expanding Chamber Results

With the settings given in the previous section, a desktop computer with 18 cores could complete the finite element analysis of the expanding chamber in less than 3 minutes (actual total elapsed time = 160 s). Results are presented in Figure 3, with the flow from right to left. During the first second ($0\le t \le 1)$, the parabolic flow profile at the left inlet rises and the flow propagates through the expanding tube until reaching the outlet. During this time frame, the hydraulic permeability of the outlet is set to its “open” value. Then, at $t=1.01$, the outlet “closes” (Figure 2) and the flow must suddenly turn back on itself in the outlet bottleneck. This becomes apparent in the movie of Figure 3 by visualizing the arrow glyphs (which represent the relative fluid velocity vector) and the streamlines. During the time span $1.01\le t \lt 2$ the chamber starts to expand due to the continuous fluid inflow from the left side. The pressure in the chamber concomitanty increases, to approximately $2\,\text{kPa}$. Then, at time $t=2$ the outlet suddenly “reopens” and the pressurized fluid in the chamber rushes out at velocities that reach a peak value of $1.6\,\text{m}/\text{s}$. Over the time span $2\le t \le 3$, there is a short period of oscillations in the wall and fluid domains, though nearly steady-state conditions are recovered before this cycle of closing and opening repeats itself.

Quarter-Symmetry Contracting Tube

Geometry and Mesh

This second example is conceptually similar to the previous one. It illustrates flow regulation during active contraction of a fluid-filled tube, to pump fluid continously in one direction. The geometry consists of a tube, 4 cm long, with an inner radius of 1 cm and a wall thickness of 2 mm. The actual model is in SI units. The inner tube domain is a fluid-FSI domain, capped on both ends by biphasic-FSI domains to represent the flow-regulating “valves”. The tube wall is an elastic solid that can also contract under a prescribed active stress oriented along the circumferential direction. There are 48 elements along the axial direction, including biased meshes in the biphasic-FSI domains to capture large gradients in fluid pressure when the flow is “turned off”. There are 16 elements in the radial direction, including a biased mesh at the fluid-wall interface to capture gradients in tangential fluid velocity. There are also 4 elements through the wall thickness. The total number of elements is 15360 (14592 hex8 and 768 penta6 elements), and the model has a total of 96282 equations.

Material Models

The fluid-FSI domain models a Newtonian viscous fluid with a mass density $\rho=1000\,\text{kg}/\text{m}^3$, shear viscosity $\mu=1\,\text{Pa}\cdot\text{s}$ and bulk modulus $K=1\,\text{GPa}$. The solid material that regularizes the mesh deformation in the fluid-FSI domain is neo-Hookean, with Young’s modulus $E=10\,\text{Pa}$ and Poisson’s ratio $\nu=0$. The wall material is a solid mixture, with mass density $\rho=1000\,\text{kg}/\text{m}^3$. It includes a neo-Hookean solid with $E=100\,\text{kPa}$ and $\nu=0$, and a prescribed uniaxial active contraction material with a peak stress of $T_0=200\,\text{kPa}$. The biphasic-FSI domain that regulates the flow at the inlet and outlet has the same fluid properties as the fluid-FSI domain. In addition, its solid matrix is neo-Hookean with $E=10\,\text{kPa}$ and $\nu=0$, and its hydraulic permeability model is const-perm-iso, with a hydraulic permeability value that steps from $k=1\,\text{m}^4/\text{N}\cdot\text{s}$ for open flow to $k=10^{-9}\,\text{m}^4/\text{N}\cdot\text{s}$ for closed flow.

Boundary Conditions

No-slip conditions were enforced between the fluid-FSI/biphasic-FSI domains and the wall domain. Zero fluid dilatation was prescribed on the inlet and outlet faces of the biphasic-FSI domains to allow free flow when the permeability was set to a high value. Solid displacements were prescribed on all three domains to enforce axisymmetric conditions for the deformations. The $Z-$displacement was fixed on the inlet and outlet boundary surfaces. The contractile stress $T_0$ was associated with a cosine function $\frac{1}{2}\left(1-\cos\left(2 \pi t\right)\right)$ with a period of $1\,\text{s}$ and repeated (Figure 5).

The hydraulic permeability of the inlet and outlet biphasic-FSI domains alternated their “closed” and “opened” configurations in synchronization (inlet) or out of synchronization (outlet) with the contractile stress, as shown in Figure 6.

Figure 6. Inlet and outlet boundary permeability load curves for tube contraction problem. The inlet is closed when the wall is contracting ($0\le t \le 0.5$, then opens while the wall contraction returns to zero ($0.5\le t \le 1$). The outlet boundary follows the opposite trend, staying open during contraction and closing during expansion.

Analysis Settings

For this model, the analysis setting requested $200$ time steps with step size of $0.01\,\text{s}$, for a total analysis time of $2\,\text{s}$. The solver employed the BROYDEN method with 5 maximum reformations and 50 maximum updates, turning off the diverge_reform and reform_each_time_step parameters for computational efficiency. For the time stepper, the maximum step size was set to $0.01\,\text{s}$ and the minimum step size was set to $10^{-8}\,\text{s}$, while the Cutback parameter was set to aggressive to accommodate the sudden “closing” and “opening” of the inlet and outlet boundaries. A maximum of 20 retries was allowed, and the optimal number of iterations was set to 50.

Results

With the settings given in the previous section, a desktop computer with 18 cores could complete the finite element analysis of the contracting tube in 14 minutes. Results are presented in Figure 7, with the flow from left to right. During the first half second ($0\le t \le 0.5)$, the wall contracts while the inlet boundary is closed and the outlet boundary is opened, causing the fluid to exit the tube through the outlet. During the next half second, the tube rebounds due to its elasticity and the decreasing contractile stress (Figure 5), while the inlet boundary is opened and the outlet boundary is closed (Figure 6), causing fluid to enter the tube from the left. This cycle repeats itself one more time before the analysis ends.

Discussion

The two case studies presented here illustrate how biphasic-FSI domains may be used to regulate the flow of a viscous fluid using valve-like functions. No actual valves are modeled explicitly here, instead the hydraulic permeability of the biphasic-FSI domain is decreased or increased in a stepwise fashion, to “close” or “open” the valve. In a strict sense, decreasing the permeability to a non-zero value does not shut the flow completely. However, in the two examples illustrated here, the reduction in flow velocity is so significant that, for all intents and purposes, the response is equivalent to reducing the flow to zero through those biphasic-FSI domains.

To achieve these computational results efficiently, there were two critical solution settings that needed to be adjusted: (1) The elasticity of the fluid-FSI domain was set to $10\, \text{Pa}$ instead of the arbitrarily small value of $10^{-9}\, \text{Pa}$ that we have previously employed in our FSI simulations. This value was still several orders of magnitude smaller than the wall elasticity ($100\, \text{kPa}$), thus having little influence on the observed flow responses. However, it helped regularize the deformation of the fluid-FSI mesh with minimal occurrences of negative jacobians. (2) The time settings Cutback parameter was set to aggressive, the minimum step size was set to a very small value, and the number of max retries was increased to 20, all of which helped the analysis sail through the sudden closing and opening of “valves”.

In both illustrated models, the opening and closing of valves was controlled by load curves. In physiological systems, one-way valves open and close in response to changes in the pressure gradient. Therefore, the illustrations provided here do not necessarily replicate the true behavior of heart valves or lymph valves. In particular, if the load curves are not timed properly with the rest of the response, it is possible to get non-physiological behavior. For example, if we reduce the fluid shear viscosity from $\mu=1\,\text{Pa}\cdot\text{s}$ to $\mu=0.004\,\text{Pa}\cdot\text{s}$ in the tube contraction problem, the fluid entering the tube during the rebound half of the cycle reaches the outlet valve before the rebound has completed, thus briefly encountering a closed outlet, which causes fluid recirculation until the time that the outlet valve opens (Figure 8). Clearly, in this case, the interaction of factors such as fluid viscosity, tube wall stiffness, contractile stress profile, and tube length causes undesirable results that are difficult to predict a priori. Therefore, users should be aware of these potential limitations when using load curves to open and close these biphasic-FSI valve-like domains.

In summary, this case study illustrates how biphasic-FSI materials may be used to regulate flow, using stepwise changes in hydraulic permeability.

Footnotes

• Shim-2021
  Shim JJ, Maas SA, Weiss JA, Ateshian GA. Finite Element Implementation of Biphasic-Fluid Structure Interactions in FEBio. J Biomech Eng. 2021 Mar 25. doi: 10.1115/1.4050646. Epub ahead of print. PMID: 33764435.
• Brinkman-1949
  Brinkman, H.C., 1949. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Flow, Turbulence and Combustion, 1(1), pp.27-34.