Computational characterization of vascular mechanics

Vascular mechanics plays a key role in both health and disease. Therefore, the mechanical properties of vessels have been under study for over a century. This thesis reports research with two computational models designed to better understand vessel mechanics in complex loading scenarios. Numerous methodologies have been utilized to evaluate the mechanical behavior of blood vessels, including distending arterial rings to investigate circumferential behavior, a configuration commonly used in wire myography. We previously used this configuration to experimentally characterize microstructural damage in cerebral arteries that may transpire in clinical procedures and due to trauma. However, due to the complexity of loading, we were not able to quantify strains throughout the vessel experimentally. As a consequence, we were not able to relate microstructural damage with vessel strains in all parts of the vessel. Thus, the aim of the current investigation was to quantify strains throughout the arterial ring by using a computational model. To achieve our goal, we created a finite element (FE) model of the experiment using FEBio. In the model, we observed complex vessel strain distributions along the circumference. Most vessel strains were observed to vary considerably through the wall thickness in regions near the needles, but circumferential strains remained largely constant throughout the ring. In this research, another computational model was constructed to understand the significance of perfusion in cerebral arteries' strain rate dependence. Although many investigators have attempted to characterize the strain rate dependence of arteries experimentally, there has been disagreement in the results. In our previous investigation, our lab observed strain rate dependence in dynamically-loaded middle cerebral arteries (MCAs) in rats. We hypothesized that perfusion was at least partly responsible for the observed behavior and designed a computational model using LS-DYNA to test our hypothesis. As expected, we observed a contribution of perfusion to strain rate dependence in the circumferential and the radial directions. However, it was not sufficient to influence experimentally witnessed axial strain rate dependence.

COMPUTATIONAL CHARACTERIZATION OF VASCULAR MECHANICS by Sameer Sanjay Nandikar A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical Engineering The University of Utah December 2017 Copyright © Sameer Sanjay Nandikar 2017 All Rights Reserved The University of Utah Graduate School STATEMENT OF THESIS APPROVAL The thesis of Sameer Sanjay Nandikar has been approved by the following supervisory committee members: , Chair Kenneth L. Monson 10/20/2017 Date Approved , Member Brittany Coats 10/20/2017 Date Approved , Member Ashley Spear 10/23/2017 Date Approved and by , Chair/Dean of Timothy Ameel the Department/College/School of Mechanical Engineering and by David B. Kieda, Dean of The Graduate School. ABSTRACT Vascular mechanics plays a key role in both health and disease. Therefore, the mechanical properties of vessels have been under study for over a century. This thesis reports research with two computational models designed to better understand vessel mechanics in complex loading scenarios. Numerous methodologies have been utilized to evaluate the mechanical behavior of blood vessels, including distending arterial rings to investigate circumferential behavior, a configuration commonly used in wire myography. We previously used this configuration to experimentally characterize microstructural damage in cerebral arteries that may transpire in clinical procedures and due to trauma. However, due to the complexity of loading, we were not able to quantify strains throughout the vessel experimentally. As a consequence, we were not able to relate microstructural damage with vessel strains in all parts of the vessel. Thus, the aim of the current investigation was to quantify strains throughout the arterial ring by using a computational model. To achieve our goal, we created a finite element (FE) model of the experiment using FEBio. In the model, we observed complex vessel strain distributions along the circumference. Most vessel strains were observed to vary considerably through the wall thickness in regions near the needles, but circumferential strains remained largely constant throughout the ring. In this research, another computational model was constructed to understand the significance of perfusion in cerebral arteries' strain rate dependence. Although many investigators have attempted to characterize the strain rate dependence of arteries experimentally, there has been disagreement in the results. In our previous investigation, our lab observed strain rate dependence in dynamically-loaded middle cerebral arteries (MCAs) in rats. We hypothesized that perfusion was at least partly responsible for the observed behavior and designed a computational model using LS-DYNA to test our hypothesis. As expected, we observed a contribution of perfusion to strain rate dependence in the circumferential and the radial directions. However, it was not sufficient to influence experimentally witnessed axial strain rate dependence. iv TABLE OF CONTENTS ABSTRACT ....................................................................................................................... iii Chapters 1. INTRODUCTION ......................................................................................................... 1 1.1 Consequences of traumatic brain injury .............................................................. 1 1.2 Vascular mechanics in TBI .................................................................................. 1 1.3 Isolated blood vessel experiments ....................................................................... 2 1.4 Computational modeling of experiments ............................................................. 3 1.5 Objective .............................................................................................................. 3 2. COMPUTATIONAL CHARACTERIZATION OF BLOOD VESSEL STRAIN DURING CIRCUMFERENTIAL STRETCHING BY WIRE MYOGRAPHY ................ 5 2.1 Background .......................................................................................................... 5 2.2 Methods................................................................................................................ 6 2.2.1 Model geometry ....................................................................................... 8 2.2.2 Mesh....................................................................................................... 11 2.2.3 Boundary conditions and contacts ......................................................... 11 2.2.4 Material model selection and verification.............................................. 14 2.2.5 Loading .................................................................................................. 16 2.2.6 Convergence study ................................................................................. 20 2.2.7 Model validation .................................................................................... 20 2.3 Results ................................................................................................................ 31 2.4 Discussion .......................................................................................................... 43 3. PERFUSION PROVIDES NEGLIGIBLE CONTRIBUTION TO STRAIN RATE DEPENDENCE IN BLOOD VESSELS .......................................................................... 47 3.1 Background ........................................................................................................ 47 3.2 Methods.............................................................................................................. 48 3.2.1 Model geometry ..................................................................................... 50 3.2.2 Mesh....................................................................................................... 53 3.2.3 Boundary conditions and contacts ......................................................... 59 3.2.4 Material model selection and verification.............................................. 65 3.2.5 Loading .................................................................................................. 71 3.2.6 Convergence study ................................................................................. 74 3.2.7 Model validation .................................................................................... 78 3.2.8 Postprocessing methods ......................................................................... 85 3.3 Results ................................................................................................................ 85 3.4 Discussion ........................................................................................................ 101 4. CONCLUSION .......................................................................................................... 117 APPENDIX:FEBIO INPUT FILES, LS-DYNA KEYWORD, AND MATLAB CODES ........................................................................................................................... 119 REFERENCES ............................................................................................................... 137 vi CHAPTER 1 INTRODUCTION 1.1 Consequences of traumatic brain injury According to the Centers for Disease Control and Prevention (CDC), traumatic brain injury (TBI) is a leading cause of death and disability in the United States. Thirty percent of all injury deaths can be attributed to TBI [1]. Survivors of TBI can suffer from its effects for an hour, a few days, or even for the rest of their lives. TBI impairs not only thinking and memory, but also sensations like vision and hearing. Patients suffering from TBI are also susceptible to personality changes and depression [2]. According to the CDC, these issues not only affect individuals, but can also cause permanent impacts on families and communities. 1.2 Vascular mechanics in TBI The grave consequences of TBI are often a result of damaged vasculature during the injury, but little is known about the mechanics of vessels during a TBI or about a vessel's response to such loading. It is well-known that overstretching arteries alters vessel mechanics and disrupts their functions [3]-[7], but the understanding of underlying structural alterations is incomplete. Blood vessels' strain rate dependence is another unresolved issue. Since TBI involves stretching arteries at high speeds, whether or not 2 they are strain rate dependent is another debated issue. A few investigators have tried to answer this question, but there is a lack of conclusive work [8]-[17]. Resolving these issues will help with more effiecient TBI diagnoses and will also lead to better safety equipment design. Surgical procedures like neuro angioplasty, in which blood vessels are distended beyond their physiological capacity, will also benefit from an improved understanding of vascular mechanics and damage. 1.3 Isolated blood vessel experiments Experiments to characterize vascular mechanics are commonly performed on isolated blood vessels from human and animal brains. In most of these experiments, the vessel segments are either clamped at the ends or mounted on hypodermic needles. Arterial segments are then stretched in the axial and/or circumferential directions. Some experiments are carried out on blood vessels perfused with a fluid to maintain pressure inside the artery segment during testing, and some tests have been performed without any perfusion. Imaging techniques are commonly used to track deformations, so as to characterize arteries' stress-strain behavior. While isolated vessel experiments have improved our understanding of vascular mechanics, these methods also have some limitations. For example, only average stresses and strains of vessel segments can be conveniently calculated by the typical experiments. It is thus difficult to characterize phenomena associated with small regions of the vessel, such as the generation of stress waves in high-rate testing. These methods are also limited by the materials and methods that can be used. For example, fiduciary markers for stretch evaluation, like microspheres, can only be placed on the outer surface of the vessel and 3 therefore cannot be used to quantify strains across the vessel wall thickness. 1.4 Computational modeling of the experiments Computational modeling of the isolated vessel experiments is an effective strategy for overcoming limitations with experimental approaches. It is comparatively easy to devise the necessary boundary conditions and loads in a computational model; due to its ability to represent a vessel segment as a cluster of tiny elements, the mechanics of any part of the vessel can be computed very efficiently. Thus, a validated model can not only resolve limitations of the experimental approach, but can also provide additional information not otherwise available in the experiments. As a consequence, computational models are commonly coupled with experimental approaches. 1.5 Objective The aim of the current study was to use computational models to more fully interpret isolated vessel experiments previously performed and investigated in our laboratory. Two models were created to explore two different experiments. In the first, middle cerebral arterial (MCA) rings from sheep were circumferentially stretched at a quasi-static speed using a modified wire myography technique. The purpose of this study and testing configuration was to understand the structural mechanics and damage of circumferentially overstretched arteries. In the previous experiments, a newly-created collagen hybridizing peptide (CHP) was utilized to identify arterial collagen damage during the circumferential stretch. However, since the full ring could not be visualized with the single camera used, we were not able to determine deformations along the entire 4 circumference or through the vessel wall thickness, in order to correlate the vessel strains with observed collagen unfolding. While wire myography is commonly used by vascular biologist to study vessel function, little is understood about the mechanics of this testing approach [18]-[23]. Thus, we created a computational model of the experiment to characterize the mechanics associated with a unique loading scenario. The second finite element model simulated axial stretching of rats' MCA. In one of our previous experiments, we observed arteries' strain rate dependence above strain rates of 500 s-1. We hypothesized that perfusion might be partly or wholly responsible for this behavior, owing to the inertia of water. The objective of this study was to investigate this hypothesis. CHAPTER 2 COMPUTATIONAL CHARACTERIZATION OF BLOOD VESSEL STRAIN DURING CIRCUMFERENTIAL STRETCHING BY WIRE MYOGRAPHY 2.1 Background Blood vessel mechanics play a vital role in both health and disease. Consequently, vessels' mechanical properties have been under study for over a century. A variety of approaches have been utilized to define blood vessel properties, including both wire and pressure myography. Wire myography is mostly used to characterize the active behavior of blood vessels through smooth muscle cells. This procedure involves intubating an arterial ring with hypodermic needles or wires. As a result, vessel rings subjected to wire myography exhibit complex strain distributions in the region around the needles. Through our literature study about wire myography, we concluded that vascular biologists use this method without considering the complex strain distributions in the vessel ring [18]-[23]. Conversely, in pressure myography, vessels are pressurized to an in-vivo state through perfusion, and the vessels' active response is studied. Due to the absence of wires or needles, pressure myography lacks the inadequacies of wire myography. However, using this method during passive response experiments, it is challenging to achieve sufficiently 6 high pressures to overstretch the vessel circumferentially. Previously, a passive response experiment was conducted in our lab to identify cerebral artery microstructural damage related to traumatic brain injury (TBI), as well as clinical procedures. It was not possible to use pressure myography, as explained above, so wire myography was utilized for that investigation. In the experiments, two hypodermic needles were inserted into a sheep's MCA ring (Figure 2.1). The ring was circumferentially distended from an unloaded state to failure at a rate of 0.1 mm/s by separating the needles quasi-statically. Strains between the needles were computed by tracking microspheres placed on the outer surface of the vessel. Collagen unfolding due to the overstretch was demonstrated using a recentlydeveloped collagen hybridizing peptide (CHP). One of the objectives of this experimental study was to correlate vessel strains with collagen damage. However, strains near the needles are complicated and differ through the vessel wall. Some portions of the vessel also move perpendicular to the imaging plane, while sliding around the needles. As a result, strains in these regions cannot be computed using microspheres alone. The goal of the computational investigation was to quantify strains throughout the vessel ring by accounting for vessel mechanics associated with wire myography. 2.2 Methods We chose a computational modeling approach to map strain distributions in vessel rings during wire myography. The FEBio software suite was used to create the model. The vessel ring was modeled as a cylinder, with symmetry utilized to reduce computational time. Loading was applied through a two-step process, first inducing 7 Figure 2.1 Experimental image of vessel ring cannulating needles and the vessel stretched circumferentially 8 residual strain in the vessel ring before stretching it circumferentially. The model was validated by relating nodal displacements and width reductions to comparable experimentally-measured values. Green-Lagrange strains of elements reported in the local cylindrical coordinate system by FEBio were plotted along vessel circumference and thickness to study strain distributions. 2.2.1 Model geometry The first step in creating a finite element model is to construct accurate geometry. The vessel ring was assumed to be a circular tube and was modeled with octant (1/8th) symmetry, such that its geometry was characterized as a quarter of an annular portion of the ring, and half its width (Figure 2.2). Experimentally-measured dimensions of a typical sheep MCA ring (outer diameter: 1.03 mm; wall thickness: 0.13 mm; width: 0.91 mm) were used to construct the geometry. The vessel ring was modeled as a homogeneous body without any distinct layers. It was initially modeled as a rectangular strip (Figure 2.3), which was then rotated into a quarter circle, to induce residual strain in the ring (Figure 2.4). With the application of symmetry, only one needle of half-length was modeled, rather than two full-length needles. The experimental setup used in the previous investigation (Figure 2.5) was measured to determine other geometry requirements of the computational model, such as the needle dimensions and the initial distance between the needles. Excluded from the model geometry were parts of the experimental setup not in contact with the vessel ring and not related to load and boundary conditions on the vessel ring. Thus, only needles (size: 28 gauge) in the experimental setup were included in the model geometry. Geometry was constructed in PreView 1.19 (University of Utah, UT). 9 Figure 2.2 Application of octant symmetry to vessel ring and needles; (a) complete geometry (highlighted portion in red was part of final geometry); (b) geometry with octant symmetry Figure 2.3 Initial geometry Figure 2.4 Final geometry 10 Figure 2.5. Experimental setup for circumferentially stretching the vessel ring, including X-Y stage, load cell, blocks for mounting needles, and voice coil actuator 11 2.2.2 Mesh PreView 1.19 was used for mesh generation and for the remaining preprocessing tasks, including the application of boundary conditions, contacts, and loads. A vessel ring was meshed using 8,400 linear hexahedral elements (Figure 2.6). Three element layers across the vessel wall thickness were formed, so as to observe strain distribution through the thickness. The needle was meshed with 65 linear hexahedral elements. Reasonable values of skewness (0), Jacobian (1), aspect ratio (10.98), and warpage angle (0) were maintained while meshing. The characteristic length of elements in the vessel ring was 4059 µm. The meshed model is as shown in Figure 2.7. 2.2.3 Boundary conditions and contacts Realistic boundary conditions and contacts were established to produce reliable results in the model (Figure 2.8). Octant symmetry boundary conditions were applied on the symmetry faces of the vessel ring, with details as shown in Figure 2.8. Rigid contact was created between the top surface of the needle and Symmetry Plane 2, preventing the vessel ring from slipping off the needle. A frictionless sliding contact was established between the top surface of the needle and the inner surface of the vessel ring, allowing the vessel ring to slide on the needle as it was stretched. A fixed contact was established between needle and the vessel ring, such that nodes of the needle and vessel ring were tied in X- and Y-directions. This contact ensured that the vessel ring didn't slip off the needle. However, the vessel ring was allowed to slide over the needle in an axial direction. All boundary conditions and contacts depicted in Figure 2.8 were maintained during both steps of the simulation (i.e., the residual stress generation step and the 12 Figure 2.6 Vessel ring mesh Figure 2.7 Model mesh 13 Fixed contact Figure 2.8 Boundary conditions to enforce octant symmetry and contact details 14 circumferential stretching step). For the first step, the needle was constrained in all degrees of freedom (DOF) except for rotation about the Z-axis, and for the second step, the needle was restricted in all DOF except for translation in the Y-direction. 2.2.4 Material model selection and verification Having a wide variety of biological material models was the main reason to choose FEBio for constructing this model. Owing to the anisotropic nonlinear nature of blood vessel tissue, the transversely isotropic Veronda-Westmann model was used for the vessel ring. It is an uncoupled material model with the following strain energy function [24], (Equation 1) (1) where Ĩ1 and Ĩ2 are the deviatoric invariants of right Cauchy-Green tensor and J is the Jacobian. C1 and C2 are material parameters which were established by fitting the material model to experimental force vs. needle displacement data (Figure 2.9). The parameter optimization function in FEBio was used for this purpose. Apart from C1 and C2, bulk modulus k for the material was also established (C1 = 0.3, C2 = 0.9, k= 4.99). The needle was modeled as a rigid body, as we were not interested in any deformations of the needle. We compared the theoretically-calculated Cauchy stress using Equations 2, 3, 4, and 5 with model prediction to verify the material model [24]. σ = pI + dev σ' σ' = (2 / J) [ (W1 + I1 W2) b' - W2 b'2 ] W1 = C1C2 ec2(I1 -3) (2) (3) (4) 15 Force on rigid body (N) 0.12 0.1 Experimental 0.08 Model 0.06 0.04 0.02 0 0 0.5 1 1.5 Needle displacement (mm) 2 Figure 2.9. Parameter optimization curve for fitting the transversely isotropic VerondaWestmann material model to experimental data 16 W2 = - (C1C2) / 2 (5) where b' is given by J-2/3 b. b is right Cauchy Green deformation tensor given by FTF. Since stresses generated in the model were complex, theoretically calculating these stresses was difficult. Thus, we created a simplified single element model and compared model predictions with theoretical calculations. This simple model included a rectangular block (width: 1 mm, length: 1 mm, thickness: 0.13367 mm) with boundary conditions as shown in Figure 2.10, and stretched in the X-direction. We compared the model prediction of Cauchy stresses and theoretically-calculated stresses in the Xdirection (Figure 2.11), and found that both matched with an average absolute error of 2.8%. 2.2.5 Loading To create a realistic representation of the blood vessel, we induced uniform circumferential residual strain in the vessel ring before stretching. Therefore, two steps of loading were applied. The presence of residual strain in blood vessels has been observed in vessels' no-load states (Figure 2.12) [25]. The opening angle ‘θ' is the measure of residual strain (Figure 2.12 b). While there was no reference in the literature about sheep MCA opening angles, we were able to find opening angles for human cerebral arteries. The opening angles for human cerebral arteries were witnessed to be variable and have been reported to vary in the range of 17-180 Deg. [26]. In some cases, angles more than 180 Deg. were also witnessed. Thus, for simplifying the model, we assumed, θ = 180 Deg. and modeled stress-free vessel ring as a flat rectangular strip (Figure 2.13). During the first step, the needle was subjected to 90 Deg. of rotation in the Z-direction to induce 17 Figure 2.10 Distribution of Cauchy stress in X direction for material verification model Stress in X-direction (MPa) 1.2 1 Theorotical Model 0.8 0.6 0.4 0.2 0 1 1.2 1.4 X-stretch Figure 2.11 Model verification comparison plot 1.6 18 Figure 2.12 Opening angle θ as a measure of residual strain in blood vessel adapted from [27] 19 (a) vessel ring needle (b) vessel ring needle Figure 2.13 Loading step 1 (legend represents Green-Lagrange strain in Y-direction), (a) vessel ring stress-free condition, (b) vessel ring in unloaded condition residual strain in the vessel ring (Figure 2.13 b) 20 residual strain. In the second step, the ring was subjected to circumferential stretch quasistatically, similar to the experiment itself, by applying 0.9 mm displacement in the Ydirection to the needle (Figure 2.14). As the model was created with symmetry about the XZ- plane, 0.9 mm displacement represents a total needle displacement of 1.8 mm. 2.2.6 Convergence study A convergence study was carried out to ensure proper mesh density. Because our chief objective was to observe hoop strain distribution along the circumference of the ring, the mesh size along the circumference was refined in the convergence study. The element size along the other dimensions of the vessel ring was adjusted such that a reasonable aspect ratio (<10) of elements was maintained. Peculiar strain distribution was observed in the portion of the vessel ring around the needle (Figure 2.15 a). Thus, average strains in that region (Figure 2.15 b) were measured for each mesh refinement, and were plotted against the number of elements (Figure 2.16). We found that results of the model mostly converged at a mesh size with 8465 elements. 2.2.7 Model validation We validated the model with experimental data using two methods. First, experimentally-measured width reduction of the ring specimen during needle displacement was compared with model predictions of width reduction at Symmetry Plane 3 (Figure 2.17). The width was measured at the top half of the circular ring. It was observed that experimental results matched with model predictions with an absolute average error of 0.94%. Second, the experimentally-measured distances between two 21 vessel ring needle Figure 2.14 Loading step 2 - circumferential stretch as a result of 1.8 mm needle displacement 22 (a) (b) Figure 2.15 Element selection for convergence study (a) hoop strain distribution along the circumference in the vessel ring portion near the needle, (b) elements selected for convergence study 23 (a) Average hoop strain 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 50000 100000 Number of elements in the model Average hoop strain 1.3 150000 (b) 1.25 1.2 1.15 1.1 1.05 1 0 5000 10000 15000 20000 Number of elements in the model Figure 2.16 Convergence plots - average strain vs. number of elements in the model (a) hoop, (b) hoop zoomed in, (c) radial, (d) radial zoomed in, (e) axial strain, (f) axial zoomed in, (g) radial-hoop shear, (h) radial-hoop shear zoomed in, (i) hoop-axial shear, (j) hoop-axial shear zoomed in, (k) radial-axial shear, and (l) radial-axial shear zoomed in 24 (c) Average radial strain -0.255 -0.26 0 50000 100000 150000 -0.265 -0.27 -0.275 -0.28 -0.285 -0.29 Average radial strain Number of elements in the model -0.28 -0.281 0 -0.282 -0.283 -0.284 -0.285 -0.286 -0.287 -0.288 -0.289 (d) 5000 10000 15000 Number of elements in the model Figure 2.16 continued 20000 25 Average axial strain 0.005 (e) 0 0 50000 100000 150000 -0.005 -0.01 -0.015 -0.02 Average axial strain Number of elements in the model 0.004 0.003 0.002 0.001 0 -0.001 0 -0.002 -0.003 -0.004 -0.005 (f) 5000 10000 15000 Number of elements in the model Figure 2.16 continued 20000 26 (g) Average rθ shear strain 0.0700 0.0600 0.0500 0.0400 0.0300 0.0200 0.0100 0.0000 0 50000 100000 150000 Number of elements in the model (h) Average rθ shear strain 0.0700 0.0650 0.0600 0.0550 0.0500 0.0450 0.0400 0 5000 10000 15000 20000 Number of elements in the model Figure 2.16 continued Average θz shear strain 27 (i) 0.0180 0.0160 0.0140 0.0120 0.0100 0.0080 0.0060 0.0040 0.0020 0.0000 0 50000 100000 150000 Number of elements in the model (j) Average θz shear strain 0.0160 0.0150 0.0140 0.0130 0.0120 0.0110 0.0100 0 5000 10000 15000 Number of elements in the model Figure 2.16 continued 20000 28 Average rz shear strain 0.0500 0.0300 0.0200 0.0100 0.0000 -0.0100 0 50000 100000 150000 -0.0200 -0.0300 Average rz shear strain (k) 0.0400 0.0000 -0.0020 0 -0.0040 -0.0060 -0.0080 -0.0100 -0.0120 -0.0140 -0.0160 -0.0180 -0.0200 Number of elements in the model 5000 10000 15000 Number of elements in the model Figure 2.16 continued 20000 (l) 29 (a) (b) Figure 2.17 Results of model validation method 1 (a) Model prediction of width at Symmetry Plane 3, (b) Experimental result of width, (c) Comparison of model predictions of width reduction and experimental results 30 (c) Vessel ring width at mid-section (mm) 1 0.9 0.8 0.7 0.6 0.5 0.4 model 0.3 experiment 0.2 0.1 0 0.11 0.61 1.11 1.61 Needle displacement in Y direction Figure 2.17 continued 2.11 31 microspheres and needles were compared with model predictions of the Y-direction distance between the needle and the two of the vessel nodes with the same initial position as those of the microspheres in the experiment (Figure 2.18). Trends in both the model and experimental results were similar. However, experimental results exhibited a small flat region initially. The reason for this could be the irregular shape of the vessel ring (i.e., not perfectly cylindrical), resulting in initial needle displacement when there was very small microsphere movement. We found that experimental results and model predictions varied with an absolute average error of 18.25% and 25.57% respectively for two nodes. The reason for this error might be our assumption about the homogeneity of the various layers of vessel material, whereas a blood vessel has three different layers (i.e., media, intima, and adventitia), which have different mechanical properties. We decided that this error was acceptable and proceeded to postprocessing. 2.3 Results The Green-Lagrange strains were recorded (Figure 2.19) and plotted to understand strain distributions in the vessel ring (Figure 2.20). As expected, modelpredicted hoop strains were the largest of all strains. We found that hoop strain was highest in the outer layer and was lowest in the inner layer, particularly where the vessel was adjacent to the needle (Figure 2.20 b). Hoop strain in the outer layer reduced along the circumference of the ring, whereas in the inner layer, hoop strain increased along the circumference. Thus, hoop strain distribution along the circumference of the vessel ring in the outer and inner layers displayed exactly opposite trends. Hoop strain distribution along the circumference of the middle layer was found to be relatively constant, with 32 (a) (b) Figure 2.18 Results of model validation method 2 (a) Model prediction of Y-direction distance between node and needle, (b) Experimental image for calculating distance between needle and microsphere, (c) Comparison of model predictions of distance between node and needle and experimental results of distance between microsphere and needle 33 Perpandicular distance between the needle and the microsphere / node (mm) 1.3 1.2 node 1 1.1 microsphere 1 1 0.9 (c) node 2 microsphere 2 0.8 0.7 0.6 0.5 0.4 0.3 0.11 0.61 1.11 1.61 Needle displacement in Y direction Figure 2.18 continued 2.11 34 Figure 2.19 Hoop strain distribution along circumference of vessel ring and vessel wall thickness on Symmetry Plane 3 35 -0.15 -0.17 (a) 0 0.2 0.4 0.6 0.8 Radial Lagrange strain -0.19 -0.21 -0.23 inner -0.25 middle -0.27 outer -0.29 -0.31 -0.33 Distance from symmetry plane 2(mm) (b) 1.7 1.6 Hoop strain 1.5 1.4 1.3 1.2 1.1 outer middle 1 inner 0.9 0.8 0 0.2 0.4 0.6 0.8 Distance from Symmetry Plane 2 Figure 2.20 Strain distributions on Symmetry Plane 3 along circumference from Symmetry Plane 2 (Figure 2.17) as a function of thickness for needle displacement of 1.8 mm (a) radial, (b) hoop, (c) axial, (d) radial-hoop shear, (e) hoop-axial shear, (f) radial-axial shear 36 (c) 0 0 0.2 0.4 0.6 0.8 -0.02 Axial strain -0.04 inner -0.06 middle -0.08 outer -0.1 -0.12 -0.14 -0.16 Distance from symmetry plane 2(mm) 0.2 (d) 0.18 inner 0.16 middle rθ shear strain 0.14 outer 0.12 0.1 0.08 0.06 0.04 0.02 0 0 0.2 0.4 0.6 Distance from symmetry plane 2(mm) Figure 2.20 continued 0.8 37 (e) θz shear strain 3.00E-03 2.50E-03 inner 2.00E-03 middle outer 1.50E-03 1.00E-03 5.00E-04 0.00E+00 0 0.2 0.4 0.6 Distance from symmetry plane 2 (mm) 0.8 1.50E-04 1.00E-04 5.00E-05 rz shear strain (f) inner middle outer 0.00E+00 0 0.2 0.4 0.6 -5.00E-05 -1.00E-04 -1.50E-04 -2.00E-04 Distance from symmetry plane 2 (mm) Figure 2.20 continued 0.8 38 little apparent influence from the model (Figure 2.20 b). Unlike hoop strain distribution, strain distributions in the other directions showed the same trend along the circumference in all three layers (i.e. radial strains, axial strains, radial-hoop shear strains, hoop-axial shear strains, and radial-axial shear strains; Figures 2.20, 2.21, and 2.22). In all three layers, strain distributions varied along the circumference in a similar manner. The radial strain was observed to be compressive. It was lowest in the inner layer and highest in an outer layer, near the needle (Figure 2.20 a). On the other hand, the radial-hoop shear strain was greatest in the middle layer and lowest in the outer layer in the region of the vessel ring near the needle (2.20 d). Radial, hoop, and radial-hoop shear strain distributions were observed to be relatively constant along the width of the vessel ring (Figure 2.20 a, b, d). Axial, hoop-axial shear and radialaxial shear strains (as shown in Figure 2.20 c, e, and f) were also uniform along the width for the most part. However, these strains had peculiar distribution along the width of the vessel ring in the region marked as Region A (Figure 2.21 c, e, f). Axial, axial-hoop shear, and radial-axial shear strains in Region A are shown in Figure 2.22. The axial strain was approximately constant through the thickness of the vessel (Figure 2.20 c), except in Region A (Figure 2.21 c), i.e., the area around the needle, near its end surface. In that region, the axial strain was highest in the inner layer and lowest in the outer layer (Figure 2.21 a). Hoop-axial and radial-axial shear strains were insignificant almost everywhere in the vessel ring (Figure 2.20 e, f) except in Region A (Figure 2.21 e, f). The hoop-axial shear strain was constant in all three layers (Figure 2.20 e), except in Region A (Figure 2.21 e). In Region A, it was highest in the inner layer and lowest in the outer layer (Figure 2.22 b). Radial-axial shear strain showed 39 Z Y X (b) Z Y X (c) Region A Z X Y X Figure 2.21 Strain distributions along width of vessel ring (inner layer) 40 Z Y X (e) Region A Z Y X (f) Region A Z Y X Figure 2.21 continued 41 0.15 (a) inner 0.1 middle Axial strain 0.05 outer 0 0 0.2 0.4 0.6 0.8 -0.05 -0.1 -0.15 -0.2 Distance from symmetry plane 2 (mm) θz shear strain 0.06 (b) 0.05 inner 0.04 middle 0.03 outer 0.02 0.01 2E-17 0 0.2 0.4 0.6 0.8 -0.01 Distance from symmetry plane 2 (mm) Figure 2.22 Strain distribution in Region A (Figure 2.20 c, e, f) along circumference as a function of thickness for needle displacement of 1.8 mm (a) Axial strain distribution, (b) Hoop-axial shear strain distribution, (c) Radial axial strain distribution 42 0.01 (c) -3E-17 -0.01 0 0.1 0.2 0.3 0.4 0.5 0.6 rz shear strain -0.02 -0.03 inner -0.04 middle -0.05 outer -0.06 -0.07 -0.08 -0.09 Distance from symmetry plane 2 (mm) Figure 2.22 continued 0.7 0.8 43 trend along the circumference of the vessel in the region around the needle; it briefly increased before becoming uniform (Figure 2.20 f). Conversely, radial-axial shear strain distribution in Region A (Figure 2.21 f) showed an increasing trend along the circumference of the vessel in the region around the needle (Figure 2.22 c). 2.4 Discussion In this investigation, we computationally modeled wire myography and achieved the primary objective of this model, i.e., to quantify strains in the region around the needles and through vessel wall thickness. The model demonstrated trends in strain distribution along the ring circumference and across the vessel wall thickness, furthering our understanding of experimentally-observed collagen damage. Since collagen fibers are oriented circumferentially in the medial layer of the blood vessel and the ring was stretched circumferentially in wire myography, collagen damage in the media was observed in the experiments. This model was mainly designed to quantify hoop or circumferential strains along the circumference to correlate with the observed collagen damage. It was interesting that circumferential strains reported by the computational model were relatively constant along the circumference of the middle layer of the ring. Remarkably, experimentally-observed hoop strains in the media were also constant along the circumference of the ring. However, it should be noted that the correspondence between the middle layer here and the media is not strict in the model. A more detailed model with layer-specific geometry, material properties, and residual strains is expected to confirm the trend of circumferential strain distribution in media, which can then be correlated with medial collagen fiber damage. 44 Interesting trends in axial and radial directions were also observed. Axial and radial strains were compressive, which was theoretically correct, as the ring was expected to shrink in radial and axial directions while being stretched in the circumferential direction. It was sensible that radial strain was lowest in an inner layer near the needle as the inner layer was expected to become more compressed by the needle. Conversely, axial strain distribution through the wall thickness was largely uniform. The absolute value of axial strains was lower near the needle and higher in the portion between the needles, which is the same as the experiment. It was interesting that the model could capture this phenomenon with frictionless contact between the needle and the ring. This trend of axial deformations was clearly visible in Region A (2.21 a and 2.22 c). Region A underwent compression in the radial direction due to the needle, which induced positive axial deformation. The plane of Region A (the surface opposite to Symmetry Plane 1) had no constrains; therefore, Region A deformations were peculiar. Radial and axial deformations in Region A were most extreme in the inner layer and least extreme in the outer layer. This difference in deformations induced extreme radial-axial and hoop-axial shear stresses in Region A. This difference in deformations could be further used to justify why we ignored diverging shear and axial strains. It was observed that of all the strains, the average hoop strain and radial strain were greatest in the region around the needle, which was selected for convergence study. Averaged hoop strain and radial strain of the finest mesh diverged from that of converged mesh with 0.82% and 0.35% absolute error, respectively. On the other hand, averaged axial strain of the finest mesh diverged from the strain of converged mesh with 63% absolute error. However, as compared to averaged hoop and radial 45 strains, values of averaged axial strain in the selected region were negligible (~ 1100 times less). The averaged axial strain in the midsection was relevant, and it diverged with only 1.35 % absolute error. Radial-hoop shear strains, hoop-axial shear strains, and radial-axial shear strains diverged with 16%, 5%, and 28% absolute error, respectively. However, these strains were almost 1/10th of the other relevant strains in the selected region. Thus, based on our observations, we concluded that even though strains slightly diverged, the mesh converged at 8465 elements. Alike normal strains, finite element model reporting shear strain distributions in the middle layer were also remarkable. Of all shear strains, the radial-circumferential shear strain was the highest, and it was maximized in the middle layer in the region of the ring around the needle. This result was in contrast with the expectation that the inner layer should have maximum radial-circumferential shear strain. Remarkably, this strain distribution through the thickness was uniform along the width of the ring. Interestingly, shear strains displayed at extremes in the region around the needles and were approximately zero in the region between the needles. It was understandable because needles induced shear deformation in the ring, whereas the portion between the needles was relatively free of shear strains. Wire myography has often been used to study the active behavior of arteries, so arterial rings are usually not circumferentially stretched to the extent that we studied in our model. However, arterial rings in such experiments are usually subjected to an initial stretch to get maximum contractile response [18]. Therefore, the finding that almost all the strains vary considerably through the vessel wall thickness near the needle suggests that experiments involving circumferential stretching, such as wire myography, should 46 consider these differences in deformations when drawing any conclusions. In the present study, the vessel ring material was assumed to be homogenous. Strain distributions observed through the thickness of the model may change notably if the model is created with layer-specific material properties, geometry, and residual strains. Thus, in the future, this model can be modified to include layer-specific details as explained above. A layer-specific model will improve correspondence between the model and a blood vessel. In this computational investigation, the ring was modeled as transversely isotropic material. This assumption may have affected the values of the strains. However, we believe that trends of strain distribution along the circumference and width of the ring would remain the same, as the trends are a function of the loading and the geometry. Thus, correspondence between the model and a vessel can be further improved by modeling the ring as nonlinear anisotropic, preferably orthotropic, material. inclusion of smooth muscle cell-related geometry and function in the model would further improve the efficacy of this model for the vascular biologist. CHAPTER 3 PERFUSION PROVIDES NEGLIGIBLE CONTRIBUTION TO STRAIN RATE DEPENDENCE IN BLOOD VESSELS 3.1 Background Traumatic brain injury (TBI) accounted for 2.5 million emergency room visits in 2010, and 2% of those cases resulted in death [1]. TBI most often leads to vascular damage due to the loading of cerebral blood vessels at high strain rates during impact. However, there is no definitive study regarding strain rate dependence of cerebral blood vessels. Initial tests on cerebral bridging veins proposed a strong influence of rates between 1 and 1000 s-1 [8], but later studies concluded that strain rate does not play a significant role in these vessels [9]-[12]. However, these subsequent studies did not include strain rates more than 250 s-1. Similarly, three more reports on strain rate dependence in cerebral arteries established an insignificance of strain rate in these vessels [13]-[15]. Remarkably, significant rate dependence in the aorta has been observed in other studies [16]-[17]. Thus, there is a need for a better understanding of strain rate significance. To resolve this issue, we conducted experiments on middle cerebral arteries (MCAs) of rats to characterize their mechanical properties during axial and stretch at different strain rates ranging from 0.05 s-1 to more than 700 s-1. 48 These experiments, which motivated our present computational investigation, were carried out on vessel segments perfused with saline. Vessel segments were mounted on two needles and were axially stretched (Figure 3.1) by moving one needle at a constant speed, while keeping the other needle stationary. A fluid column was connected to the moving needle for perfusion while the flow of fluid out of the stationary needle was restricted. Vessel segments were pressurized at different pressures (13.3 kPa and 6.7 kPa), by maintaining the height of the fluid column. Stress-stretch responses were recorded and compared to study rate dependence. No change in the axial stress-stretch response was observed in response to changes in pressure or at strain rates below 500 s -1. However, strain rate dependence in the axial direction was witnessed for strain rates above 500 s -1 (Figure 3.2). It has been proposed that rate dependence observed in cerebral vessels may be attributed to perfusion of the vessel rather than to viscoelasticity of the vessel wall (Lee and Haut, 1989). Our hypothesis was that perfusion would primarily offer rate-dependent resistance to circumferential deformation of the vessel, which in turn would induce axial strain rate dependence in the vessel. Simillarly, we also hypothesised that perfusion plays a role in the threshold of 500 s-1. The objective of the current investigation was to use a computational model to test these hypotheses. 3.2 Methods A computational model of the experiment was created to study the effect of perfusion on strain rate dependence observed in our previous investigation. LS-DYNA was used to construct and solve the model. A blood vessel segment perfused with fluid 49 Figure 3.1 Axially-stretched rat MCA segment in our previous investigation Figure 3.2 Representative axial Cauchy stress-stretch curves for all groups, including high (>700 s-1, HR), medium (400-500 s-1; MR), and low (100-200 s-1; LR) strain rate cases. The internal pressure was fixed at 6.7 kPa. These representative cases suggest trends toward higher stresses and lower stretches with higher strain rates; there was no apparent effect of pressure (from an unpublished investigation in our lab) 50 was modeled. The simulation consisted of two steps: first, pressurization of the blood vessel; second, axial stretch. An arrangement representing a fluid column was built to pressurize the vessel at 13.3 kPa pressure and to keep the vessel perfused during the axial stretch. Arbitrary Lagrange Eulerian (ALE) formulation was utilized to simulate the fluid-structure interactions in the model. Axial stretch was applied at different strain rates, and corresponding stress-strain responses were compared to quantify the effect of perfusion on stress-strain responses. 3.2.1 Model geometry Model geometry consisted of vessel segment, void, fluid reservoir, and needles (Figure 3.3). Averages of cross-sectional dimensions (outer diameter: 0.25 mm; wall thickness: 0.04 mm) and lengths (1 mm) of rat MCA segments tested in the experiments, were used to create the geometry of the blood vessel. The blood vessel was assumed to be a circular tube, and quadrant symmetry was utilized to reduce computational time (Figure 3.4). LS-PrePost was used to create the geometry, which was required to accommodate saline (fluid) and structure (vessel segment) interaction. Since LS-DYNA requires an empty mesh in which fluid may flow, a void was included in the geometry (Figure 3.5). As with the blood vessel, the void was modeled as a quarter cylinder (diameter: 0.4811 mm; length: 1.8 mm). The diameter of the void was established based on the maximum diameter of the blood vessel during pressurization. The void was extended on both sides of the vessel segment to ensure the presence of fluid in the vessel during stretch and to facilitate fluid flow similar to the experiment. 51 Figure 3.3 Line diagram of entire model (isometric view) 52 (a) (b) Figure 3.4 Blood vessel segment geometry, (a) isometric view and (b) top view (a) (b) Figure 3.5 Void (a) isometric view and (b) top view 53 To simulate pressurization by the fluid column, a reservoir of square crosssections (side: 0.028 mm; length: 0.1mm) was included in the model geometry. The reservoir was located at the end of the void on the moving end of the vessel segment. The shape, size, and location of the fluid reservoir (Figure 3.6) in the model were selected using trial and error methods. Shape, size, and location of the fluid reservoir in the model were varied, and the time for pressurizing the vessel to 13.3 kPa was recorded. Dimensions and location of the reservoir with minimum computational cost were selected to generate accurate initial conditions as well as to create a realistic representation of the fluid flow in the experimental setup. Rigid needles connected to either end of the vessel segments simulated needles, allowing fluid to flow in and out of the blood vessel (Figure 3.7). To accommodate fluid flow during axial stretch, void was needed to be longer than the vessel segment and larger in the cross-section to allow the vessel to expand during pressurization. Hence, without needles, fluid would have gone below the vessel segment. The purpose of the needles was to keep the fluid over the vessel surface or inside the lumen. Cross-sectional dimensions of the needles were the same as that of the blood vessel. The length of the needle on the stationary end of the vessel was arbitrarily selected as 0.2 mm to simulate blocked passage to the fluid, as in the experiment. The span of the needle on the moving end of the vessel was 0.6 mm to accommodate the axial displacement of 0.5 mm, as well as the reservoir of 0.1 mm length. 3.2.2 Mesh LS-PrePost was used to mesh the geometry. The vessel segment was meshed using 1260 hexahedral solid elements, as they can fully capture three-dimensional 54 (b) (a) Figure 3.6 Geometry (a) reservoir and (b) location of reservoir in the void (a) (b) Figure 3.7 Needles on either side of the vessel segment (a) Isometric view and (b) top view 55 states of stress and are ideal for modeling thick parts [28]. The constant stress solid element form (ELEFORM 1) in LS-DYNA was used. This formulation is underintegrated, yet efficient and accurate [28]. Since vessel segments underwent severe deformations, this choice was ideal. Belytschok-Bindmen strain co-rotational stiffness hourglass control (type 6) with larger values (0.9) was suggested to work better for anisotropic material; therefore, it was utilized in this model for a blood vessel segment [29]. Needles were also modeled with hexahedral solid elements (type 1; Figure 3.8 and 3.9). They were modeled as rigid bodies, so no hourglass control was used. As for the vessel segment, the void and fluid reservoir were also modeled with hexahedral solid elements (Figures 3.10, 3.11, and 3.12) because the fluid-structure interaction mechanism in LS-DYNA requires fluid to be modeled with solid elements. Overlapping nodes of the fluid reservoir and the void were merged to allow fluid to flow into the void. The one point an ALE multimaterial element form (ELEFORM 11) was utilized, allowing both to have more than one material in them. Elements of the fluid reservoir were defined for constant pressure (AET=4), to enable it to act as a continuous source of fluid at the prescribed pressure. Elements of the void were assigned to have variable pressure (AET=0), the same as the rest of the model. Due to the absence of distortion, hourglass deformation is not an issue in the case of ALE elements [28]. Thus the standard LS-DYNA viscous hourglass formation (type 1) with very low hourglass coefficient (1.0E-6) was used for all ALE elements. In summary, under-integrated elements are vulnerable to nonphysical modes of deformation, but that can be limited by using hourglass stabilization. On the other hand, under-integrated elements are computationally less costly and robust. Hence, in this 56 (a) (b) Figure 3.8 Mesh of vessel segment (number of elements along length: 70; number of elements along circumference: 6; number of elements along vessel wall thickness: 3) (a) isometric view and (b) top view (a) (b) Figure 3.9 Mesh of stationary needle (number of elements along length: 14; number of elements along circumference: 6; number of elements along vessel wall thickness: 3) (a) isometric view and (b) top view 57 (a) (b) Figure 3.10 Mesh of stationary needle (number of elements along length: 42; number of elements along circumference: 6; number of elements along vessel wall thickness: 3) (a) isometric view and (b) top view (a) (b) Figure 3.11 Mesh of void (number of elements along length: 126; number of elements on cross-section: 93) (a) isometric view and (b) top view 58 (a) (b) Figure 3.12 Mesh of fluid reservoir (number of elements along length: 7; number of elements on cross-section: 4) 59 model, under-integrated elements with appropriate hourglass control were used. 3.2.3 Boundary conditions and contacts As our objective was to study the effects of perfusion, the primary challenge of this investigation was to realistically model fluid-structure interactions between the fluid and the surrounding vessel segment. The CONSTRAINED_LAGRANGE_IN_SOLID (CLIS) keyword was used to establish contact between the fluid and the vessel segment. The inside surfaces of the vessel segment and the two needles were specified as slave surfaces (Figure 3.13), and the void mesh was defined as the master, while the ALE material of the fluid reservoir was assigned to be in contact with slave surfaces. A penalty-type contact was selected. A penalty curve (Figure 3.14) was assigned to this coupling to counter leakage of fluid; other measures recommended in LS-DYNA literature were also employed to avoid leakage. Leakage sites were visually identified during the simulation, and maximum pressure at those locations was recorded, so as to assign the proper value of pressure in penalty curve. Next the "dbfsi" file outputted by LS-DYNA was studied to find out the amount of leakage, as well as the ratio of coupling forces to leakage control forces. The CLIS card was modified to keep leakage and the forces ratio minimum at optimum computational cost. Similar to the coupling mechanism, LS-DYNA has another card (CONTROL_ALE) to control advection of fluid material. Donor cell + HIS advection method with the new algorithm for continuum treatment were selected in the CONTROL_ALE card based on model requirements and recommendations in LS-DYNA 60 Figure 3.13 Inside surface of vessel segment and needles defined as slave for fluid structure interaction 61 Figure 3.14 Penalty curve assigned to coupling mechanism to apply coupling pressure proportional to depth of fluid leakage 62 literature [28],[29]. As we were modeling fluid flow, smoothing was turned off as recommended in LS-DYNA literature [28],[29], in order to save computational cost. While CLIS assumes a "slip" boundary condition between fluid and structure, the automatic Euler boundary condition (EBC) option, present in CONTROL_ALE , alternatively enabled us to apply a "no-slip" boundary condition. We conducted our entire investigation using the "slip" condition and did some additional tests with the "noslip" boundary condition to understand the significance of fluid flow on axial strain rate dependence of the vessel segment. Along with advection control, boundary conditions were specified to simulate the realistic movement of fluid in the model. As per quadrant symmetry, the flow of fluid out of the X-symmetry face (Figure 3.15 a) was restricted by constraining the nodes on that surface in the X-direction. Similarly, nodes on the Y-symmetry face were arrested in the Y-direction to prevent fluid flow (Figure 3.15 b). The nodes on the end of the void on the stationary needle were constrained in the Z-direction (Figure 3.15 c) to simulate blocked flow as in the experimental setup. Similarly, quadrant symmetry conditions were used on blood vessel symmetry faces (Figure 3.16 a and b). The stationary needle was constrained in the X-, Y-, and Z-direction displacements and rotations, while the moving needle was constrained in the X-, Y-, and Z-direction rotations as well as X- and Y-direction displacements. The moving needle was free to move in the Z-direction. A fixed rigid contact was established between the ends of the vessel segments and the needles by merging their nodes. Another challenge of the model was to simulate flow in and out of a blood vessel experimental setup, fluid flowed in and out of the moving end of the blood vessel. By 63 (a) (b) (c) Figure 3.15 Restrict the fluid flow out of (a) X symmetry face, (b) Y symmetry face, (c) bottom surface of void 64 (a) (b) 3.16 Symmetry boundary conditions on (a) X symmetry face and (b) Y symmetry face 65 trial and error, the inflow was defined by creating a reservoir with 13.5 kPa pressure as if the fluid column with approximately 13.3 kPa pressure was connected to it. At the moving end of the vessel-as explained in above sections-the outflow was defined by applying a pressure boundary condition of 13.8 kPa on the moving end of the void (Figure 3.17). Fluid flowed out of the reservoir and into the vessel continuously until the reservoir pressure was less than its surroundings and fluid flowed out of the void when the pressure of the fluid increased more than the specified value. Elements of the void inside the vessel segment and needles were filled up with fluid from the reservoir at the beginning of the simulation. INITIAL_VOLUME_FRACTION_GEOMETRY keyword was used for this purpose. 3.2.4 Material model selection and verification Material model selection involved selecting material models and parameter values for the blood vessel, fluid, void, and needles. In addition, LS-DYNA requires the Equation of State (EOS) to be defined for ALE materials. Appropriate materials were selected from available material models in LS-DYNA. LS-DYNA has a full range of anisotropic material models, including those for orthotropic materials, but most are linear, which is not consistent with our requirement for modeling a blood vessel. Instead, we required a nonlinear anisotropic material model to replicate the behavior of cerebral arteries. A compromise was made by selecting a transversely isotropic material model (MAT_SOFT_TISSUE), including an isotropic Mooney-Rivlin matrix reinforced by fibers that has a strain energy contribution with the qualitative material behavior of collagen [30]. The model also has a viscoelasticity 66 3.17 Pressure boundary condition (13.3 KPa) on top of void 67 option, which may be useful in the future, but this was not utilized in the present investigation. This model is based on the work of Weiss et al. [1996] and Puso and Weiss [1998] [29]. The overall "uncoupled" strain energy (W) function for the material model is given in Equation 6 [30]. W = C1 ( Ĩ1 - 3) + C1 ( Ĩ2 - 3) + F (λ) + 0.5 K [ln (J)] 2 (6) where Ĩ1 and Ĩ2 are the deviatoric invariants of the right Cauchy deformation tensor, and λ is the deviatoric part of the stretch along the current fiber direction. The straightening of fibers (i.e., before a critical stretch limit - λ < λ*) is described by an exponential function, whereas the behavior of straightened fibers past the critical stretch limit (λ ≥ λ*) is defined by a linear function. K is bulk modulus, and J =detF is the volume ratio. Shear strains were assumed to be insignificant. Incompressibility was enforced by including Lagrange multiplier p in stress equations. Corresponding average Cauchy stress values are expressed as given below [24],[31] tzz = 2 [ ( λzz2 - 1 / (λzz2 λcc2) ) C1 - ( (λzz2 λcc2) - 1 / λzz2 ) C2 ] + λWλ a a tcc = 2 [ ( λcc2 - 1 / (λzz2 λcc2) ) C1 - ( (λzz2 λcc2) - 1 / λcc2 ) C2 ] (7) (8) where tzz and tcc are the theoretical Cauchy stresses in the axial and circumferential directions. We assumed radial stresses to be negligible based on our previous work with human blood vessels [14],[15]. The term "a" is the unit vector in the fiber direction. For the present investigation, we assumed that the fibers of interest were oriented axially. Therefore, in the model, ‘a' was [0,0,1] and term a a was 0 for stresses in the circumferential and radial direction. As a consequence, the term λWλ was non-zero only for axial stresses. The term λWλ described fiber behaviors , i.e., unstretched, stretched up to critical stretch, and stretched beyond critical stretch λWλ delineated stretching of fibers at different stages of the stretch as given below [31]. 68 λWλ = 0, λ < 1 λWλ = C3 (exp ( C4 ( λ - 1) ) - 1 ), λ < λ* λWλ = C5 λ + C6 , λ ≥ λ* (9) We assumed λ* to be 1.2, based on our previous experiments. Exponential stresses were scaled by C3, and the collagen fiber uncrimping rate was determined by C4. C5 was the modulus in the linear function of straightened collagen fiber. C6 was characterized as shown below in Equation 10 to maintain continuity between the exponential and linear portions of the model30. C6 = C3 (exp ( C4 ( λ* - 1) ) - 1 ) - C5 λ* (10) Quasi-static experimental data were used to fit the above-defined stress equations because our hypothesis was that the experimentally-observed rate dependence in blood vessels was caused by luminal fluid inertia rather than by viscoelasticity inherent to the vessel material. In our experimental investigation, average curves were plotted, and data were grouped for constitutive model fitting. The group of data used for the present study included a quasi-static axial stretch test at ~13.3 KPa and ~6.7 KPa pressure, as well as quasi-static circumferential tests at medium and low pressures. Average experimental axial and circumferential stresses (Tzz and Tcc ) were compared with theoretical values. A constrained nonlinear regression routine (fminsearchbnd) in MATLAB (Mathworks, Natick, MA) was utilized to minimize the objective function f (Equation 11) in order to find material parameters C1 , C2, C3, C4, and C5 (units MPa) that best fit the experimental data (Figure 3.18). (11) where N is representative of a total number of data points, and the theoretical stresses are 69 (a) (b) Figure 3.18 Curve fitting for getting material parameters (a) circumferential stresses from circumferential tests and axial stresses from axial tests, (b) Axial stress from circumferential tests circumferential stress from axial tests 70 expressed in Equations 7 and 8. Based on the results of curve-fitting, we established parameters C1 =0.014953 MPa, C2 = 0.006494 MPa, C3 = 0.007929 MPa, C4 = 16.678829. However, C5 = 2.4146941 MPa was established by extrapolating experimental data linearly and fitting it to linear stress function. In addition, density (ρ) and bulk modulus (K) of blood vessels were assumed to be 1.075e-003 g/mm 3 and 22000 MPa, respectively. These values were selected based on the assumption that a blood vessel has a density of and bulk modulus the same as water. Similar to the vessel segment, appropriate material models were assigned to other parts in the model. MAT_NULL was assigned to the fluid reservoir with a fluid density 0.000998 g / mm3 [32], cavitation pressure =- 24 MPa [33], and viscosity = 1e-006 MPa.ms [34]. In addition, EOS_Gruneisen, which relates the change in pressure of the fluid to the change in its corresponding specific internal energy, was specified for fluid. The speed of sound in fluid was set as 1647 mm/ms, S 1 = 1.92, S2 = -0.092, gamma= 0.35 [28]. For no change in density, the pressure of the fluid can be calculated as a function of gamma and specific internal energy [28]. As we wanted to observe strain rate dependence of vessels with the 13.3 kPa pressure, we selected value of e 0 =0.038517 N-mm/g, to achieve pressure 13.5 kPa for the reservoir [28]. The void was assigned MAT_VACUUME, with a density of air 1.25e-006 g/mm 3. It is a dummy material which represents a vacuum in the multimaterial ALE model [30] (i.e., elements can have more than one material in them, and any material can flow in and out of MAT_VACCUME). MAT_RIGID_BODY was assigned to the needles. To verify our above-defined material parameters for the vessel segment, we conducted material model verification. As we were only verifying the blood vessel 71 material model, this simulation consisted of the vessel material only having the same shape and dimensions as specified in Section 3.2.1. Mesh and symmetry boundary conditions of the vessel segment were the same as explained in Sections 3.2.3 and 3.2.3, respectively. For verification, the one end of the vessel segment was fixed in the axial direction (Figure 3.19 a) and a quasi-static displacement of 0.5 mm at 0.02 m/s second rate was applied at the other end (Figure 3.19 b). Due to the quasi-static rate, the axial stress distribution (Figure 3.19 c) was uniform. The average axial stress of all the elements in the model was calculated. Axial and circumferential stretch values were extracted from the model, and axial stress was theoretically calculated using Equations 2 and 3. Model and theoretical stress-stretch responses were compared by plotting both model and theoretical axial stresses against axial stretch (Figure 3.20 a). We found that they both match with an insignificant error. To confirm the strain rate independence of the vessel material, we also plotted the average axial stress-strain response at three strain rates (1000 s-1, 500 s-1, 100 s-1). As expected, the material model was strain-rateindependent (Figure 3.20 b) 3.2.5 Loading The aim of this investigation was to replicate high-rate axial stretch experiments in a computational model, but it was important to apply the axial stretch from the same initial state as that of the experiment. In that case, the initial state of the vessel was predominantly defined by the pressure of the luminal fluid. As explained in Section 3.2.3, pressurization of the vessel was carried out using the fluid reservoir and pressure boundary condition. Zero displacement was applied to the 72 (a) (b) (c) Figure 3.19 Verification model (a) boundary condition to axially fixed vessel segment, (b) axial stretch applied at the end of the vessel segment, (c) axial stress distribution 73 (a) (b) Figure 3.20 Response of the vessel segment (a) comparison of theoretical calculations and model response of axial stress-stretch; (b) strain rate independent response of vessel segment 74 needle and pressure at the center of the void was calculated and plotted against time to determine the time step at which 13.3 kPa pressure was attained (Figure 3.21). Elements in the "slip" model were pressurized in 4.1125 ms, while the "no-slip" model required 22.8 ms. Once the desired pressure of the fluid was reached (Figure 3.22), 0.5 mm displacement in the axial direction was applied to the moving needle while the stationary needle was kept fixed, thereby axially stretching (λmax = 1.5) the blood vessel (Figure 3.23). This loading was applied at 100 s-1 strain rate (velocity = 0.1 m/s), 500 s-1 strain rate (velocity = 0.5 m/s), and 1000 s-1 strain rate (velocity = 1 m/s) (Figure 3.24). Tests with fluid having density 100 times that of water (0.0998 gm / mm 3) were carried out to determine the effect of perfusion fluid density on strain-rate dependence as well as on overall stress-strain response. This model was run at the three strain rates and was referred as the high-density model here onwards. 3.2.6 Convergence study A convergence study was carried out to determine the mesh size beyond which further mesh refinement would produce the same results to confirm consistency of model response, irrespective of change in mesh density. The mesh size along the axis of the blood vessel was varied from coarsest to finest. For mesh size along the circumference and through the vessel, the thickness was maintained such that the aspect ratio was less than 5. The mesh size of the void along the axis and circumference was kept the same as that of the vessel, to maintain suitable coupling between fluid and structure [30],[31]. First, a simulation was run with just a pressurization step, without any axial stretch, to get the time step at which 13.3 kPa pressure was achieved. The simulation was later rerun, as 75 (a) (b) Figure 3.21 Pressurization illustration (a) elements at the center of the void selected for pressure calculations; (b) average pressure of selected elements plotted against time (for 1000 s-1 strain rate model) (a) (b) Figure 3.22 Pressurized geometries (a) vessel segment and (b) column of fluid at the end of the pressurization stage 76 (a) Figure 3.23 Axially-stretched vessel at 1000 s-1 strain rate (1.3 stretch) (a) isometric view and (b) back view (b) 77 (a) (b) Figure 3.24 Axial displacement load curve for model for (a) ‘slip' model and (b) ‘no-slip' model 78 explained in Section 3.2.5, with a modified load curve of the moving needle so as to apply axial stretch after that time step. The strain rate for all simulations in the convergence study was 1000 s-1 to save computational time. Average axial stresses and strains, as well as mean circumferential and radial stresses and strains at 1.3 stretch over the flat portion of the pressurized artery (Figure 3.25), were evaluated and plotted against a number of elements to study convergence (Figure 3.26). We considered the flat portion of the artery to eliminate end effects arising from connection to the rigid needles. Mesh converged at 13986 elements with a vessel segment having 70 elements along the length, 6 elements along the circumferenc, and 3 elements through the thickness. Averaged axial, hoop, and radial stress of the finest mesh diverged from that of converged mesh with 2.09%, 3.13%, and 2.34% absolute error, respectively. Averaged axial, hoop and radial strain of the finest mesh diverged from that of converged mesh with 1.3 %, 3.11%, and 2% absolute error, respectively. Averaged axial strain didn't diverge. 3.2.7 Model validation In order to have confidence in the results of a finite element model, validation with experimental results is required. Since material model parameters of the artery were found by fitting the material model to quasi-static test data, the simulation for validation was run at relatively low 20 s-1 strain rate (0.02 m/s); slower speeds were not feasible given the large required computation times. The average axial and circumferential stress over the midsection of the vessel (as explained in Section 3.2.6) was plotted against axial and circumferential stretch, respectively. Experimental data from our previous investigation with rat MCA were used for this comparson. We used quasi-static test data. 79 (a) (b) Figure 3.25 Elements selected for convergence study in the flat portion of the vessel (a) back view and (b) isometric view Average axial stress (MPa) 80 (a) 0.60 0.50 0.40 0.30 0.20 0.10 0.00 0 10000 20000 30000 40000 50000 60000 70000 Average circumferential stress (MPa) Number of elements in the model (b) 9.00E-02 8.00E-02 7.00E-02 6.00E-02 5.00E-02 4.00E-02 3.00E-02 2.00E-02 1.00E-02 0.00E+00 0 20000 40000 60000 80000 Number of elements in the model Figure 3.26 Convergence plots at 1.3 stretch of average stress / strain vs. number of elements in the model (a) axial stress, (b) circumferential stress, (c) radial stress, (d) axial strain, (e) circumferential strain, (f) radial strain 81 Average radial stress (MPa) 0.000 0 20000 40000 60000 80000 (c) -0.005 -0.010 -0.015 -0.020 Number of elments in the model (d) 0.40 Average axial strain 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0 20000 40000 Number of elments in the model Figure 3.26 continued 60000 80000 82 (e) Average circumferential strain 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 0 20000 40000 60000 80000 Number of elments in the model 0.00 Average radial strain -0.10 0 20000 40000 60000 -0.20 -0.30 -0.40 -0.50 -0.60 -0.70 Number of elments in the model Figure 3.26 continued (f) 80000 83 The model response was compared against the the experimental axial and circumferential stress-stretch response (Figure 3.27). In an earlier investigation, results from 12 quasistatic experiments were reported. We designed the model to predict the response of the vessel only up to the yield point. Thus, results of experiments up to yield were used to validate the model. As vessel segments reached yield point at different stretch values, the stress-stretch response of all of the experiments was plotted, instead of using an average response of all experiments. It was observed in the experiments that vessels yielded mostly around the 1.3 axial stretch. Thus, we validated the model only up to that point. Our observations of the model were also based on an axial stretch up to 1.3. In general, the model predicted a stiffer axial response than that of the experiments. While the model exhibited maximum absolute errors up to 53% in comparison to experimental data, the overall behavior was deemed to provide a good match of the data. For example, at 1.35 axial stretch, the minimum absolute error in axial response was 13% and maximum absolute error was 59%. Similarly, circumferential response of the model (3.27 b) deviated from the experimental response, but it was in the same range as that of the 12 experiments. At 1.3 axial stretch for the circuferential stretch minimum, absolute error was 7% and maximum absolute error was 19%. Similarly for circumferential stress, minimum and maximum absolute errors were respectively 15% and 39%. The material parameters of the computational model were determined based on an average of the 12 experiments; thus, model predictions were close to average exeperimental response. In addition, the diameter of the vessel segment at 13.3 kPa pressure (0.34 mm) was approximately the same as the experimentally-measured diameters (0.33-0.35 mm) of the vessel segments at the same pressure. 84 (a) (b) Figure 3.27 Validation of model with experimental data 85 3.2.8 Postprocessing methods As expected, dynamic loading produced a strain wave in the vessel segment (Figure 3.28). Therefore, we studied the results of individual elements in the flat portion of the vessel segment (Figure 3.29) instead of averaging the stress-strain response. Elements were selected such that they were representative of axial, radial, and circumferential distribution of vessel segment characteristics. However, the axial stressstrain response of all the elements was found to be similar. Therefore, for comparisons between different models, axial stress-strain response of element-Axial 1 was used. Moreover, in case of circumferential and radial stress-strain responses, individual element data were noisy. Hence, for all those cases, average stress-strain response was calculated. We were interested in the axial and circumferential stress-strain response of the vessel segment. Since the vessel segment was assumed to be cylindrical and was axially stretched, we were able to establish 2nd principal stress-strain as circumferential stressstrain and 3rd principal stress-strain as radial strain. 3.3 Results In this investigation, the stress-strain response of the vessel was plotted and compared with the three strain rates to study strain-rate dependence in the axial, circumferential, and radial directions. We observed axial deformations of the vessel segment to be similar (Figure 3.30) for all three strain rates (1000 s-1, 500 s-1, 100 s-1). The axial stress-strain response of all the elements for all the strain rates was witnessed to be nearly uniform (Figure 3.31). The highest absolute difference between stresses corresponding to highest stretch was 86 Figure 3.28 Elements selected for postprocessing Figure 3.29 Strain wave in the dynamically stretched vessel segment for model with 1000 s-1 strain rate 87 (b) (a) (c) Figure 3.30 Deformation of vessel segment at 1.3 stretch with (a) 1000 s -1 strain rate; (b) 500 s-1 strain rate; (c) 100 s-1 strain rate 88 (a) (b) Figure 3.31 Axial stress - strain response of model (with slip boundary condition) for strain rate dependence study (a) element-a, (b) element-b, (c) element-g, (d) element-c, (e) element-d, (f) element-f, (g) element-e 89 (c) (d) Figure 3.31 continued 90 (e) (f) Figure 3.31 continued 91 (g) Figure 3.31 continued 92 observed at element-e (Figure 3.31 e) and it was approximately 8%. However, this difference was observed between 100 s-1 and 1000 s-1 strain rates; on the other hand, absolute difference between responses at 1000 s -1 and 500 s-1 was relatively negligible. Similar trends were witnessed in all the other elements. Also, a small strain-rate dependence in the averaged axial stress-strain response was witnessed in the high-density model (Figure 3.32). At 1.3 stretch, the difference in stresses at 1000 s-1 and 100 s-1 was approximately 2 %. Furthermore, similar to the normal density model, difference between 1000 s-1 and 500 s-1 model was relatively negligible. Here we calculated averaged stress-strain response as individual response of the high density model was noisy. Similarly, a little change in axial stress-strain response was observed when results of the model with perfusion were compared with the results of the model without perfusion (Figure 3.33). The stresses recorded in the perfused model and the model without perfusion differed by 2% at the end of the axial stretch and they were higher in the perfused model. Perfusion didn't change the trend of the stress-strain response. Axial strain vs. time was plotted for element "a," for the perfused model, the model without perfusion, and the high-density model, to understand the trend of axial deformation during axial stretch at 1000 s-1 strain rate (Figure 3.34). It was observed that axial deformation was tensile, as expected, and was the same for all three models. Unlike axial stress-strain response, circumferential stress-strain response varied considerably with strain rate (Figure 3.35). The difference in circumferential deformations at different strain rates could also be seen in Figure 3.30. In Figure 3.30, it was observed that for 1000 s-1 and 500 s-1 strain rates, the circumferential deformation was not uniform along the length of the vessel segment, while it was almost uniform for 93 (a) (b) Figure 3.32 Axial stress-strain of high density model averaged over midsection (a) overall respone, (b) zoomed in response at the end of the stretch 94 Figure 3.33 Element-a comparison of stress-strain response of the model with perfusion and the model without perfusion Figure 3.34 Axial deformation trend for element-a at 1000 s -1 strain rate 95 Figure 3.35 Average circumferential stress-strain response of the vessel for strain-rate dependence study, 96 the 100 s-1 strain rate case. This phenomenon was more prominently witnessed in deformations of the high-density model (Figure 3.36). The average circumferential stressstrain response of the high-density model at the three strain rates was as shown in Figure 3.37. In addition, results of models with and without perfusion were compared to examine the effect of perfusion on the circumferential stress-strain response. It was observed that without perfusion, the circumferential strain was negative with the almost constant trend of circumferential stress-strain response (Figure 3.38). Conversely, with perfusion, circumferential strain was observed to be positive with an increasing trend of circumferential stress-strain response. To verify the circumferential stress trend, we plotted circumferential stress versus circumferential stretch of one of the axial-stretch experiments carried out at 987 s-1 strain rate (Figure 3.39). Except for the initial response (up to circumferential strain 0.1), the circumferential stress-strain trend reported by the model was largely similar to experimental observations. Circumferential strain versus time was plotted for element "Axial 1," for the perfused model, the model without perfusion, and the high-density model, to understand the trend of circumferential deformation during axial stretch at 1000 s-1 strain rate (Figure 3.40). Radial stress was assumed to be zero while determining material model parameters (Section 3.2.4). Hence, we calculated average radial stresses and strains of the midsection of the vessel to check the validity of our assumption. Radial stresses were compressive and negligible as compared to the axial and the circumferential stresses (Figure 3.41). The radial stress-strain response was calculated at three strain rates, and strain rate dependence was witnessed in the radial direction (Figure 3.42). The absolute value of radial stress and strain increased with an increase in strain rates. A similar trend 97 (b) (a) (c) Figure 3.36 Deformation of vessel segment perfused with fluid having density 0.0998 g / mm 3 at 1.3 stretch with (a) 1000 s-1 strain rate; (b) 500 s-1 strain rate; (c) 100 s-1 strain rate 98 Figure 3.37 Average circumferential stress-strain response of the vessel in the high-density model for strain rate dependence study Figure 3.38 Comparison of average circumferential stress-strain response of the model with perfusion and model without perfusion 99 Figure 3.39 Circumferential stress vs. stretch response from an axial stretch experiment at 987 s-1 strain rate Figure 3.40 Circumferential deformation trend for element - "a" at 1000 s -1 strain rate 100 Figure 3.41 Element a - comparison of axial, circumferential, and radial stress-strain responses (strain rate : 1000 s-1) Figure 3.42 Radial stress-strain response at the three strain rates to determine strain rate dependance in radial direction 101 was observed in the radial stress-strain response of the high-density model (Figure 3.43). Radial strain vs. time was plotted for element "a," for the perfused model, the model without perfusion, and the high-density model to understand the trend of radial deformation during axial stretch at 1000 s-1 strain rate (Figure 3.44). We also computed axial reaction forces on the needle for the three strain rates and plotted those against axial stretch. We observed that like axial stress, axial reaction force was also independent of strain rate influence. Like stress-strain responses, pressure traces of the model for the three strain rates were calculated as explained in Section 3.2.5 and were plotted against axial stretch. It was found that the pressure of the fluid inside the lumen was slightly higher at 1000 s -1 and 500 s-1 strain rates than at 100 s-1 (Figure 3.45). We also compared axial stress-strain responses of the model with the ‘slip' boundary condition and the model with the ‘no-slip' boundary condition. Axial stressstrain responses of the seven elements (Figure 3.28) for both ‘slip' and ‘no-slip' models at different strain rates were compared (Figure 3.46-3.53). We found that the axial stressstrain response was almost the same for both boundary conditions at all strain rates and all over the vessel geometry. However, the ‘slip' boundary condition case was computationally less costly than that of the ‘no-slip' boundary condition. Thus, results regarding strain rate dependence presented in this investigation were with ‘slip' boundary condition. 3.4 Discussion In this research, we used a computational model to test our hypothesis that perfusion directly resists circumferential deformation and the resulting circumferential 102 Figure 3.43 Average radial stress-strain response of the vessel in the high-density model for strain rate dependence study Figure 3.44 Radial deformation trend for element-a at 1000 s -1 strain rate 103 Figure 3.45 Effect of strain rate on pressure trace Figure 3.46 Axial reaction force vs. axial stretch 104 (a) (b) (c) Figure 3.47 Axial stress-strain response of element-a (Figure 3.28) at strain rate (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 105 (a) (b) (c) Figure 3.48 Axial stress-strain response of element-b (Figure 3.28) at strain rate (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 106 (a) (b) (c) Figure 3.49 Axial stress-strain response of element-g (Figure 3.28) at strain rate (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 107 (a) (b) (c) Figure 3.50 Axial stress-strain response of element-c (Figure 3.28) at strain rate (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 108 (a) (b) (c) Figure 3.51 Axial stress-strain response of element-d (Figure 3.28) at strain rates (a) 1000 s -1, (b) 500 s-1, (c) 100 s-1 109 (a) (b) (c) Figure 3.52 Axial stress-strain response of element-f (Figure 3.28) at strain rates (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 110 (a) (b) (c) Figure 3.53 Axial stress-strain response of element-e (Figure 3.28) at strain rates (a) 1000 s-1, (b) 500 s-1, (c) 100 s-1 111 strain rate dependence influences axial strain rate dependence. The model shows that this effect has negligible influence on results, suggesting that fluid inertia is not the source of rate dependence observed in vessel stretch experiments. According to our hypothesis, we observed the effect of perfusion and strain rate on the circumferential stress-strain response of the vessel segment (Figure 3.34). One cause for such behavior can be the inertia of fluid during the dynamic axial stretch. The effect of inertia could be understood when the behavior of the model with perfusion was compared with the behavior of the model without perfusion. In the model without perfusion, the vessel segment stretched in the axial direction and shrank in the circumferential and radial directions when axial stretch was applied at any rate. Thus, in the model without perfusion, circumferential strain was negative and reduced as the vessel was stretched even at the high strain rates. However, in the perfused model at the high strain rates, circumferential strains increased as the fluid in the lumen could not shrink in diameter as fast as the axial deformation of the artery. Particularly, the fluid in the lumen of the midsection of the vessel didn't deform as fast as the vessel segment. The increase in the circumferential strains may have been result of displacement of the fluid into the midsection from adjacent sections of the lumen, which deformed relatively fast. Fluid present in the midsection may have offered resistance to circumferential deformation. Thus, circumferential stress also increased along with circumferential strain. Consequently, more resistance was offered to circumferential deformation due to inertia of the fluid at high strain rates. This behavior was witnessed more prominently in the high-density Model (Figure 3.35), but its stress-strain response was largely similar to the stress-strain response of the normal density model. Nevertheless, that perfusion was 112 responsible for this phenomenon was confirmed by comparing corresponding results of the perfused model and the model without perfusion. Circumferential stress-strain in the model without perfusion was very negligible and compressive in nature. On the other hand, circumferential stress-strain in the perfused vessel model was tensile and significant. Only at 100 s-1 strain rate, circumferential stress-strain started to reduce at the end of the axial stretch. For 1000 s-1 strain rate, circumferential stress-strain response had an increasing trend. This observation supported our earlier explanation about less resistance to circumferential deformation at low strain rates. Therefore, perfusion was responsible for strain rate dependence in the circumferential direction. However, from our results, clearly perfusion did not play a significant role in axial strain rate dependence as observed in our previous experimental investigation. Even with an unrealistically high-density fluid in the high-density model, we were not able to see any significant strain rate dependence. A small shift in stress-strain response was witnessed when axial stress-strain responses of models with and without perfusion were compared. This shift was due to pressurized lumen in the perfused model as opposed to a pressure-free blood vessel in the model without perfusion. Axial stress increased by approximately 2% at the end of the stretch as result of pressure in the lumen. Thus, for our investigation, perfusion offered a small resistance to axial deformation, but this resistance was not strain rate dependent. According to the Generalized Hook's law, strain rate dependence in a circumferential direction would be expected to indirectly influence the strain rate dependence in the axial direction. However, we were not able to determine the exact cause of the absence of axial strain rate dependence in our model. It might be because 113 circumferential stresses were comparatively much smaller than axial stresses (approximately 10 times smaller). The amount of strain rate induced by circumferential stresses was even smaller than that. Thus, axial strain induced due to strain rate dependent circumferential stresses might be negligible, which resulted in no axial strain rate dependence. However, more study is required to prove this or to find better explanations. Like circumferential direction, the radial stress-strain response was also influenced by perfusion. An increase in compressive radial strain could be seen in Figure 3.44. However, strains in the high-density model were largely similar to strains in the normal density model. Moreover, strain rate dependence was observed in radial directions as well. Again, the reason for that may be inertia of fluid, as explained above. At higher strain rates, the bulk of water present in lumen compressed the vessel wall, thus absolute values of both radial strains and radial stresses were the highest at 1000 s -1. Absolute values of radial strains and stresses were lower for lower strain rates as less fluid was present in the midsection lumen at low strain rates. Thus, there was less compression of the vessel wall. Our investigation with slip and no-slip boundary conditions between the vessel wall and the fluid lead us to believe that flow characteristics do not have any effect on strain rate dependence. Thus, if there is any strain rate dependence due to perfusion, it can only be because of the inertia of fluid. We observed strain rate dependence beyond 500 s-1 in our earlier experimental investigation. We hypothesized that at strain rates below 500 s-1, the resistance offered by perfusion to circumferential deformations, and in turn axial deformations, was low, 114 but that the resistance offered by perfusion at strain rates above 500 s-1 was considerable enough to have a notable influence. However, in the computational model, we observed that perfusion plays no role. Another explanation for the experimentally-observed rate dependence could be due to viscoelasticity of the vessel wall, but more investigation will be required to understand the mechanism of such influence. In conclusion, our hypothesis was true about the contribution of perfusion in circumferential strain rate dependence. However, the circumferential stresses were insignificant; thus, axial stresses influenced by circumferential stresses were too small to cause axial strain rate dependence. While, this is one possible explanation, further study is required to confirm this hypothesis or to find a better explanation. In addition, the density of the fluid did not influence stress-strain behavior in the axial direction, which eliminated the density of the fluid as a factor contributing to axial strain rate dependence. Nevertheless, a parametric study with larger diameters of vessel segments will help exclude the presence of fluid in the lumen as a cause of axial strain rate dependence. In the present computational investigation, the vessel segment was modeled as a transversely isotropic material. In the future, the vessel can be modeled as an orthotropic material to improve correspondence of the model with vessel tissue. In modeling the vessel segment as a homogenous, transversely isotropic material, we assumed axial direction as the principal direction of fibers for the entire vessel segment. This is true for adventitia. However, principal fiber direction in media is circumferential. Hence, layerspecific modeling of the vessel segment will further improve correspondence between a blood vessel and the model. In addition, a model with viscoelastic material properties of the vessel segment will be required to further understand strain rate dependence. 115 Moreover, modeling the fluid as blood, instead of as saline, may also improve model relevance. Relative viscosity of plasma at 37 Deg. Celsius is about 1.8 [35]. Thus, for a model with blood, a no-slip boundary condition would be established between the vessel segment and the fluid. Based on our observation, we believe that this would not change the results of the simulation regarding axial strain rate dependence. However, such study would be required to establish the effects of the change in viscosity of the fluid. Our future investigations may include other factors that might have played a role in the different conclusions regarding strain rate dependency of blood vessels. Apart from perfusion, other variables that might have contributed to different conclusions could be the range of strain rates at which the experiments were conducted, different animals that were used in different experiments, and different organs from which blood vessels were dissected. For example, initial tests on human cerebral bridging veins suggested a strong effect of rates between 1 and 1000 s-1 [8], but in later studies on human cerebral bridging veins conducted at rates below 250 s-1, no strain rate dependence was observed [9]-[12]. In another example, studies performed on human and pig thoracic aortas showed remarkable influence of strain rate [16],[35] while experiments performed on human cerebral arteries didn't report any strain rate dependence [13]-[15]. The reason for different conclusions in these cases might be the different physical constitutions of these blood vessels in different organs of different animals. Apart from the factors explained above, the method of specimen preparation could also have contributed to different conclusions. For example, experiments performed on cerebral blood vessels were conducted on round blood vessels [13]-[15] whereas experiments performed on thoracic 116 aortas were carried out on I-shaped specimens [16]-[17]. Determining which of these factors influence strain rate dependence would further improve our understanding of TBI. CHAPTER 4 CONCLUSION We created two computational models of isolated blood vessel experiments in this investigation. In the first model, we wanted to quanitfy vessel strains throughout the vessel to correlate with microstructural damage observed in an arterial ring subjected to circuferential stretching by wire myography. We found out that vessel strains vary considerably through the vessel wall and along the circumference. Near the needles, vessel strain distributions were particularly complex. However, hoop strains in the middle layer of the ring were largely uniform along the circumference. Interestingly, experimentally-observed collagen damage was also uniform in the middle layer, i.e., media. However, we could not correlate the hoop strains in the middle layer with microstructural damage observed in media, as the ring was modeled as homogenous material. Thus, a model with layer-specific material proporties, geometry, and residual strain will further improve the correspondance of the model and a blood vessel. In this model, the arterial ring was modeled as transversly isotropic material; orthotropic representation of the ring would make the model more realistic and its results more relatable. This finite element model will be more effective to a vascular biologist if smooth muscle cells were included in the arterial ring. The second computational model was designed to determine the contribution of 118 perfusion to the strain rate dependence of blood vessels. The model showed that perfusion didn't contribute to axial strain rate dependence. We identified insignificant values of circumferential stresses as one possible cause of a negligible axial strain rate dependence. Even though perfusion played a role in circumferential and radial strain rate dependence, insignificant values of both the stresses influenced the insignifcant axial strain rate dependent response. However, more study and additional tests are required prove this hypothesis. The density of the fluid was not a contributing factor in axial strain rate dependence. However, a parametric study with large diameters of the vessel segments would be required to deterimine the role played by the vessel size in perfusioninduced strain rate dependence. In the model, the vessel segment was fashioned as transversly isotropic and homogenous. Both of these assumptions require further improvement. An orthotropic vessel segment with layer-specific principal fiber directions will improve the correspondance between the model and the vessel tissue. Relevance of the model would be further improved by modeling blood as the fluid. In the future, a study with a viscoelastic vessel segment will further our understanding of strain rate dependence of blood vessels. In conclusion, computational modeling of experiments furthered our understanding about vascular mechanics. However, limitations of our computational models, as well as interesting findings from the models themselves, created new questions. Therefore, further efforts are required to better understand vascular mechanics in isolated blood vessel experiments. APPENDIX FEBIO INPUT FILE, LS-DYNA KEYWORD, AND MATLAB CODES 120 FEBio input file <?xml version="1.0" encoding="ISO-8859-1"?> <febio_spec version="2.0"> <Globals> <Constants> <T>0</T> <R>0</R> <Fc>0</Fc> </Constants> </Globals> <Material> <material id="1" name="blood vessel" type="trans iso VerondaWestmann"> <density>1</density> <c1>0.3</c1> <c2>0.9</c2> <k>4.99997</k> <c3>0</c3> <c4>0</c4> <c5>0</c5> <lam_max>0</lam_max> <fiber type="local"> 0, 0</fiber> </material> <material id="2" name="needle" type="rigid body"> <density>1</density> <center_of_mass>0,0,0</center_of_mass> </material> <material id="3" name="square bar" type="rigid body"> <density>1</density> <center_of_mass>0.247835,-0.878595,0</center_of_mass> </material> </Material> <Geometry> <Contact> <contact type="facet-to-facet sliding" name="sliding cotact between inner surface of blood vessel and outer surface of needle"> <laugon>1</laugon> <tolerance>0.2</tolerance> <penalty>400</penalty> <two_pass>1</two_pass> <auto_penalty>0</auto_penalty> <fric_coeff>0</fric_coeff> <fric_penalty>0</fric_penalty> <search_tol>0.01</search_tol> <minaug>0</minaug> <maxaug>10</maxaug> <gaptol>0</gaptol> <seg_up>0</seg_up> 121 <surface type="master"> </Contact> <Constraints> </Constraints> <LoadData> <loadcurve id="1" type="smooth"> <point>0,0</point> <point>1,1</point> </loadcurve> <loadcurve id="2" type="smooth"> <point>1,1</point> <point>2,1</point> </loadcurve> <loadcurve id="3" type="linear"> <point>0,0</point> <point>1,0</point> <point>2,1</point> </loadcurve> </LoadData> <Output> <plotfile type="febio"> <var type="contact gap"/> <var type="contact pressure"/> <var type="contact traction"/> <var type="displacement"/> <var type="reaction forces"/> <var type="stress"/> </plotfile> </Output> <Step name="Step01"> <Module type="solid"/> <Control> <time_steps>10</time_steps> <step_size>0.1</step_size> <max_refs>15</max_refs> <max_ups>10</max_ups> <dtol>0.001</dtol> <etol>0.01</etol> <rtol>0</rtol> <lstol>0.9</lstol> <time_stepper> <dtmin>0.01</dtmin> <dtmax>0.1</dtmax> <max_retries>5</max_retries> <opt_iter>10</opt_iter> </time_stepper> <analysis type="static"/> </Control> <Boundary> </Boundary> 122 </contact> <Contact> </Contact> <Constraints> <rigid_body mat="2"> <fixed bc="x"/> <fixed bc="y"/> <fixed bc="z"/> <fixed bc="Rx"/> <fixed bc="Ry"/> </rigid_body> <rigid_body mat="3"> <fixed bc="z"/> <fixed bc="Rx"/> <fixed bc="Ry"/> <fixed bc="Rz"/> </rigid_body> <rigid_body mat="2"> <prescribed bc="Rz" lc="1">1.57143</prescribed> </rigid_body> </Constraints> </Step> <Step name="Step02"> <Module type="solid"/> <Control> <time_steps>10</time_steps> <step_size>0.1</step_size> <max_refs>15</max_refs> <max_ups>10</max_ups> <dtol>0.001</dtol> <etol>0.01</etol> <rtol>0</rtol> <lstol>0.9</lstol> <time_stepper> <dtmin>1e-009</dtmin> <dtmax>0.1</dtmax> <max_retries>5</max_retries> <opt_iter>10</opt_iter> </time_stepper> <analysis type="static"/> </Control> <Boundary> <fix bc="uvw"> </Boundary> <Constraints> <rigid_body mat="2"> <fixed bc="x"/> <fixed bc="z"/> <fixed bc="Rx"/> 123 <fixed bc="Ry"/> </rigid_body> <rigid_body mat="2"> <prescribed bc="Rz" lc="2">1.57143</prescribed> </rigid_body> <rigid_body mat="2"> <prescribed bc="y" lc="3">0.8982</prescribed> </rigid_body> </Constraints> </Step> </febio_spec> 124 LS-DYNA Keyword $# LS-DYNA Keyword file created by LS-PrePost(R) V4.3 - 22Dec2016(09:00) $# Created on Jul-10-2017 (16:35:48) *KEYWORD *TITLE $# title LS-DYNA keyword deck by LS-PrePost *CONTROL_ALE $# dct nadv meth afac bfac cfac dfac efac -1 1 1 -1.0 0.0 0.0 0.0 0.0 $# start end aafac vfact prit ebc pref nsidebc 0.01.00000E20 1.01.00000E-6 0 0 0.0 0 $# ncpl nbkt imascl checkr 1 50 0 0.0 *CONTROL_ENERGY $# hgen rwen slnten rylen 2 2 2 2 *CONTROL_MPP_DECOMPOSITION_DISTRIBUTE_ALE_ELEMENTS *CONTROL_STRUCTURED *CONTROL_TERMINATION $# endtim endcyc dtmin endeng endmas 4.6125 0 0.0 0.01.000000E8 *CONTROL_TIMESTEP $# dtinit tssfac isdo tslimt dt2ms lctm erode ms1st 0.0 0.9 0 0.0 0.0 0 0 0 $# dt2msf dt2mslc imscl unused unused rmscl 0.0 0 0 0.0 *DATABASE_ELOUT $# dt binary lcur ioopt option1 option2 option3 option4 0.0025 0 0 1 0 0 0 0 *DATABASE_GLSTAT $# dt binary lcur ioopt 0.0025 0 0 1 *DATABASE_MATSUM $# dt binary lcur ioopt 0.0025 0 0 1 *DATABASE_BINARY_D3PLOT $# dt lcdt beam npltc psetid 0.0025 0 0 0 0 125 $# ioopt 0 *DATABASE_EXTENT_BINARY $# neiph neips maxint strflg engflg 0 0 3 1 1 $# cmpflg ieverp beamip dcomp ialemat 0 0 0 1 1 $# nintsld pkp_sen sclp hydro nodout 0 0 1.0 0 $# dtdt resplt neipb 0 0 0 *DATABASE_FSI $# dt 0.0025 $#dbsfi_id sid stype swid 1 2 2 0 *BOUNDARY_PRESCRIBED_MOTION_RIGID $# pid dof vad lcid birth 5 3 2 3 0.0 *BOUNDARY_SPC_SET $# nsid cid dofx dofy dofrz 1 0 1 0 1 *SET_NODE_LIST_TITLE NODESET(SPC) 1 $# sid da1 da2 da3 1 0.0 0.0 0.0 *BOUNDARY_SPC_SET $# nsid cid dofrz 2 0 1 *SET_NODE_LIST_TITLE NODESET(SPC) 2 $# sid da1 2 0.0 sigflg epsflg rltflg 1 1 1 shge stssz n3thdt 1 1 2 msscl therm intout 0 0 convid 0 ndsetid 0 cid 0 sf vid death 1.0 01.00000E28 dofz dofrx dofry 0 1 1 da4 solver 0.0MECH dofx dofy dofz dofrx dofry 0 1 0 1 1 da2 0.0 da3 0.0 da4 solver 0.0MECH 126 MATLAB code to fit material model to experimental data %This file was originally written by E. David Bell as an example of how %Matlab can be used to fit experimental data using a user defined strain %energy function to fit, as well as a user defined objective function. clc; clear; close all global lamC lamZ Tcc Tzz LB UB % Load sampled data load QS_CircData_avg load QS_AxialData_avg % Group data into typical groupings (LamC, LamZ, Tcc, Tzz) S7 = cat(2,S7_LamC_avg,S7_LamZ_avg,S7_Tcc_avg,S7_Tzz_avg); S13 = cat(2,S13_LamC_avg,S13_LamZ_avg,S13_Tcc_avg,S13_Tzz_avg); % S20 = S20_sampled; Plow = cat(2,Plow_LamC_avg,Plow_LamZ_avg,Plow_Tcc_avg,Plow_Tzz_avg); Pmed = cat(2,Pmed_LamC_avg,Pmed_LamZ_avg,Pmed_Tcc_avg,Pmed_Tzz_avg); % Phi = Phi_sampled; % Pressure Tests (1=Plow; 2=Pmed; 3=Phigh;) % Axial Tests (4=S7; 5=S13; 6=S20;) % Circumferential Stretch Data lamC1 = Plow(:,1); lamC2 = Pmed(:,1); lamC3 = []; lamC4 = S7(:,1); lamC5 = S13(:,1); lamC6 = []; % Circumferential Stress Data Tcc1 = Plow(:,3); Tcc2 = Pmed(:,3); Tcc3 = []; Tcc4 = S7(:,3); Tcc5 = S13(:,3); Tcc6 = []; % Axial Stretch Data lamZ1 = Plow(:,2); lamZ2 = Pmed(:,2); lamZ3 = []; lamZ4 = S7(:,2); lamZ5 = S13(:,2); lamZ6 = []; % Axial Stress Data Tzz1 = Plow(:,4); Tzz2 = Pmed(:,4); Tzz3 = []; Tzz4 = S7(:,4); Tzz5 = S13(:,4); Tzz6 = []; %% Fit pressure and axial test data 127 lamC = cat(1,lamC1,lamC2,lamC3,lamC4,lamC5,lamC6); % circumferential stretch lamZ = cat(1,lamZ1,lamZ2,lamZ3,lamZ4,lamZ5,lamZ6); % axial stretch Tcc = cat(1,Tcc1,Tcc2,Tcc3,Tcc4,Tcc5,Tcc6); % circumferential stress Tzz = cat(1,Tzz1,Tzz2,Tzz3,Tzz4,Tzz5,Tzz6); % axial stress x0=[1,1,1,1]; LB=[0.0149530,0.0064940,-inf,-inf]; UB=[inf,inf,inf,inf]; options options options options options options options options = = = = = = = = optimset('Largescale','off'); optimset(options,'Algorithm','levenberg-marquardt'); optimset(options,'Display','iter'); optimset(options,'FunValCheck','on'); optimset(options,'MaxFunEvals',500000); optimset(options,'MaxIter',500000); optimset(options,'TolFun',1e-12); optimset(options,'TolX',1e-12); [x1,fval1]=fminsearchbnd(@Objective_Exponential_Function_Fung3,x0,LB,UB,optio ns);; % Create fit data coeffs = x1; % Stress fit data from pressure tests (1-3) [YfitC1 YfitZ1] = CalcYfitsFung(coeffs,lamC1,lamZ1); [YfitC2 YfitZ2] = CalcYfitsFung(coeffs,lamC2,lamZ2); [YfitC3 YfitZ3] = CalcYfitsFung(coeffs,lamC3,lamZ3); % Stress fit data from axial tests (4-6) [YfitC4 YfitZ4] = CalcYfitsFung(coeffs,lamC4,lamZ4); [YfitC5 YfitZ5] = CalcYfitsFung(coeffs,lamC5,lamZ5); [YfitC6 YfitZ6] = CalcYfitsFung(coeffs,lamC6,lamZ6); % Plot 1: Circum Stresses from P-Tests (1-3)and Axial Stresses from STests % (4-6) subplot(2,2,1) % h1 = Figure(1); plot(lamC1,Tcc1,'ro',lamC1,YfitC1,'r.-'); hold on plot(lamC2,Tcc2,'bo',lamC2,YfitC2,'b.-'); hold on plot(lamC3,Tcc3,'ko',lamC3,YfitC3,'k.-'); hold on plot(lamZ4,Tzz4,'rx',lamZ4,YfitZ4,'r.-'); hold on plot(lamZ5,Tzz5,'bx',lamZ5,YfitZ5,'b.-'); hold on plot(lamZ6,Tzz6,'kx',lamZ6,YfitZ6,'k.-'); legend('Plow Exp','Plow Fit',... 'Pmed Exp','Pmed Fit',... 'S7 Exp','S7 Fit',... 'S13 Exp','S13 Fit',... 'Location','NorthEast'); ylabel('Cauchy Stress (MPa)') xlabel('Stretch') title('Model 3A: Circum Stresses from Circum Tests and Axial Stresses from Axial Tests') % saveas(h1,'Fung Plot 1.png') 128 % Plot 2: Axial Stress from P-Tests (1-3) and Tests % (4-6) subplot(2,2,2) % h2 = Figure(2); plot(lamZ1,Tzz1,'ro',lamZ1,YfitZ1,'r.-'); plot(lamZ2,Tzz2,'bo',lamZ2,YfitZ2,'b.-'); plot(lamZ3,Tzz3,'ko',lamZ3,YfitZ3,'k.-'); plot(lamC4,Tcc4,'rx',lamC4,YfitC4,'r.-'); plot(lamC5,Tcc5,'bx',lamC5,YfitC5,'b.-'); plot(lamC6,Tcc6,'kx',lamC6,YfitC6,'k.-'); legend('Plow Exp','Plow Fit',... 'Pmed Exp','Pmed Fit',... 'S7 Exp','S7 Fit',... 'S13 Exp','S13 Fit',... 'Location','NorthEast'); ylabel('Cauchy Stress (MPa)') xlabel('Stretch') title('Model 3A: Axial Stress from Circum from Axial Tests') % saveas(h2,'Fung Plot 2.png') Circum Stresses from S- hold hold hold hold hold on on on on on Tests and Circum Stresses x1 fval1 % %% Fit only pressure test data % % lamC = cat(1,lamC1); % circumferential stretch % lamZ = cat(1,lamZ1); % axial stretch % Tcc = cat(1,Tcc1); % circumferential stress % Tzz = cat(1,Tzz1); % axial stress % % x0=[1,1,1,1]; % % options = optimset('Largescale','off'); % options = optimset(options,'Algorithm','levenberg-marquardt'); % options = optimset(options,'Display','iter'); % options = optimset(options,'FunValCheck','on'); % options = optimset(options,'MaxFunEvals',500000); % options = optimset(options,'MaxIter',500000); % options = optimset(options,'TolFun',1e-12); % options = optimset(options,'TolX',1e-12); % % [x2,fval2]=fminsearch(@Objective_Exponential_Function_Fung3,x0,options); % % % Create fit data % coeffs = x2; % % Stress fit data from pressure tests (1-3) % [YfitC1 YfitZ1] = CalcYfitsFung(coeffs,lamC1,lamZ1); % % % Plot 3: Circumferential Stress from Pressure Tests (1-3) % h3 = Figure(3); % subplot(1,2,2) % plot(lamC1,Tcc1,'ko',lamC1,YfitC1,'r-'); hold on % legend('Piv Exp','Piv Fit','Location','NorthWest'); 129 % % % % % % % % % % % % % % % % % % ylabel('Cauchy Stress (MPa)') xlabel('Stretch') title('Fung Model - Circum Stress from P-Tests') saveas(h3,'Fung Plot 3.png') % Plot 4: Axial Stress from Pressure Tests (1-3) h4 = Figure(3); subplot(1,2,1) plot(lamZ1,Tzz1,'ko',lamZ1,YfitZ1,'r-'); hold on legend('Piv Exp','Piv Fit','Location','NorthWest'); ylabel('Cauchy Stress (MPa)') xlabel('Stretch') title('Fung Model - Axial Stress from P-Tests') % saveas(h4,'Fung Plot 4.png') x2 fval2 %% Fit only stretch test data lamC = cat(1,lamC4,lamC5); % circumferential stretch lamZ = cat(1,lamZ4,lamZ5); % axial stretch Tcc = cat(1,Tcc4,Tcc5); % circumferential stress Tzz = cat(1,Tzz4,Tzz5); % axial stress x0=[1,1,1,1]; options options options options options options options options = = = = = = = = optimset('Largescale','off'); optimset(options,'Algorithm','levenberg-marquardt'); optimset(options,'Display','iter'); optimset(options,'FunValCheck','on'); optimset(options,'MaxFunEvals',500000); optimset(options,'MaxIter',500000); optimset(options,'TolFun',1e-12); optimset(options,'TolX',1e-12); [x3,fval3]=fminsearch(@Objective_Exponential_Function_Fung3,x0,options); % Create fit data coeffs = x3; % Stress fit data from stretch tests (4-6) [YfitC4 YfitZ4] = CalcYfitsFung(coeffs,lamC4,lamZ4); [YfitC5 YfitZ5] = CalcYfitsFung(coeffs,lamC5,lamZ5); % % Plot 5: Axial Stress from Stretch Tests (4-6) h5 = Figure(5); subplot(2,2,3) plot(lamZ4,Tzz4,'rx',lamZ4,YfitZ4,'r.-'); hold on plot(lamZ5,Tzz5,'bx',lamZ5,YfitZ5,'b.-'); legend('S7 Exp','S7 Fit','S13 Exp','S13 Fit','Location','NorthWest'); ylabel('Axial Stress (MPa)') xlabel('Axial Stretch') title('Model 3B: Axial Stress from Axial Tests') % saveas(h5,'Fung Plot 5.png') 130 % Plot 6: Circumferential Stress from Axial Tests (4-6) subplot(2,2,4) plot(lamC4,Tcc4,'rx',lamC4,YfitC4,'r.-'); hold on plot(lamC5,Tcc5,'bx',lamC5,YfitC5,'b.-'); legend('S7 Exp','S7 Fit','S13 Exp','S13 Fit','Location','NorthWest'); ylabel('Circumferential Stress (MPa)') xlabel('Circumferential Stretch') title('Model 3B: Circumferential Stress from Axial Tests') % saveas(h6,'Fung Plot 6.png') x3 fval3 function [YfitC, YfitZ] = CalcYfitsFung(coeffs,lamC,lamZ) c = coeffs(1); c1 = coeffs(2); c2 = coeffs(3); c3 = coeffs(4); YfitC = ((2.*((lamC.^2) 1./((lamC.^2).*(lamZ.^2)))).*c)+((2.*(((lamC.^2).*(lamZ.^2)) 1./(lamC.^2))).*c1); YfitZ = ((2*(lamZ.^2 1./((lamC.^2).*(lamZ.^2))))*c)+((2*(((lamC.^2).*(lamZ.^2)) 1./(lamZ.^2)))*c1)+((exp(c2*(lamZ-1))-1)*c3); function [x,fval,exitflag,output]=fminsearchbnd(fun,x0,LB,UB,options,varargin) % FMINSEARCHBNDNEW: FMINSEARCH, but with bound constraints by transformation % % Changes from fminsearchbnd: % 1) in options structure, user may pass an 'output function' and 'plot function' to fminsearch. % Original fminsearchbnd handled the output function via a nested wrapper function. I have extended % this to the plot function too. % 2) I have moved the 'intrafun' and 'xtransform' functions and wrappers to be nested functions % (INSIDE the fminsearchbnd function), so they do not need to pass the params structure around % (into fminsearch) - but have access to it directly. This maintains the integrity of the varargin, % which the user may be passing thru fminsearch to their optmization funciton (fminsearchbnd had % passed the params structure to fminsearch, thus ruining any varargin that the user passed in). % This also obviates the params.(whatever) structure the author had, so I've eliminated it so things % are simpler. % 3) I have created a test example so the user can see not only how fminseachbnd works, but also how % the OutputFn and PrintFns functions work, which were heretofore poorly documented by MathWorks. % Many thanks to the original author, John D'Errico, for excellent work very useful! % 131 % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % Modifications by: Ken Purchase Email: kpurchase at yahoo Date: 2007-Nov-29 usage: usage: usage: usage: usage: usage: x=FMINSEARCHBND(fun,x0) x=FMINSEARCHBND(fun,x0,LB) x=FMINSEARCHBND(fun,x0,LB,UB) x=FMINSEARCHBND(fun,x0,LB,UB,options) x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...) [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...) arguments: fun, x0, options - see the help for FMINSEARCH LB - lower bound vector or array, must be the same size as x0 If no lower bounds exist for one of the variables, then supply -inf for that variable. If no lower bounds at all, then LB may be left empty. Variables may be fixed in value by setting the corresponding lower and upper bounds to exactly the same value. UB - upper bound vector or array, must be the same size as x0 If no upper bounds exist for one of the variables, then supply +inf for that variable. If no upper bounds at all, then UB may be left empty. Variables may be fixed in value by setting the corresponding lower and upper bounds to exactly the same value. Notes: If options is supplied, then TolX will apply to the transformed variables. All other FMINSEARCH parameters should be unaffected. Variables which are constrained by both a lower and an upper bound will use a sin transformation. Those constrained by only a lower or an upper bound will use a quadratic transformation, and unconstrained variables will be left alone. Variables may be fixed by setting their respective bounds equal. In this case, the problem will be reduced in size for FMINSEARCH. The bounds are inclusive inequalities, which admit the boundary values themselves, but will not permit ANY function evaluations outside the bounds. These constraints are strictly followed. If your problem has an EXCLUSIVE (strict) constraint which will not admit evaluation at the bound itself, then you must provide a slightly offset bound. An example of this is a function which contains the log of one of its parameters. If you constrain the 132 % % % % % % % % % % % % % % % % % % % % % % % % % variable to have a lower bound of zero, then FMINSEARCHBND may try to evaluate the function exactly at zero. Example: rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; fminsearch(rosen,[3 3]) ans = 1.0000 1.0000 % unconstrained fminsearchbnd(rosen,[3 3],[2 2],[]) ans = 2.0000 4.0000 % constrained See test_main.m for other examples of use. See also: fminsearch, fminspleas Author: John D'Errico E-mail: [e-mail] Release: 4 Release date: 7/23/06 % size checks xsize = size(x0); x0 = x0(:); xLength=length(x0); if (nargin<3) || isempty(LB) LB = repmat(-inf,xLength,1); else LB = LB(:); end if (nargin<4) || isempty(UB) UB = repmat(inf,xLength,1); else UB = UB(:); end if (xLength~=length(LB)) || (xLength~=length(UB)) error 'x0 is incompatible in size with either LB or UB.' end % set default options if necessary if (nargin<5) || isempty(options) options = optimset('fminsearch'); end % % % % 0 1 2 3 --> --> --> --> unconstrained variable lower bound only upper bound only dual finite bounds 133 % 4 --> fixed variable BoundClass = zeros(xLength,1); for i=1:xLength k = isfinite(LB(i)) + 2*isfinite(UB(i)); BoundClass(i) = k; if (k==3) && (LB(i)==UB(i)) BoundClass(i) = 4; end end % transform starting values into their unconstrained % surrogates. Check for infeasible starting guesses. x0u = x0; k=1; for i = 1:xLength switch BoundClass(i) case 1 % lower bound only if x0(i)<=LB(i) % infeasible starting value. Use bound. x0u(k) = 0; else x0u(k) = sqrt(x0(i) - LB(i)); end % increment k k=k+1; case 2 % upper bound only if x0(i)>=UB(i) % infeasible starting value. use bound. x0u(k) = 0; else x0u(k) = sqrt(UB(i) - x0(i)); end % increment k k=k+1; case 3 % lower and upper bounds if x0(i)<=LB(i) % infeasible starting value x0u(k) = -pi/2; elseif x0(i)>=UB(i) % infeasible starting value x0u(k) = pi/2; else x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1; % shift by 2*pi to avoid problems at zero in fminsearch % otherwise, the initial simplex is vanishingly small x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k)))); end % increment k k=k+1; case 0 134 % unconstrained variable. x0u(i) is set. x0u(k) = x0(i); % increment k k=k+1; case 4 % fixed variable. drop it before fminsearch sees it. % k is not incremented for this variable. end end % if any of the unknowns were fixed, then we need to shorten % x0u now. if k<=xLength x0u(k:xLength) = []; end % were all the variables fixed? if isempty(x0u) % All variables were fixed. quit immediately, setting the % appropriate parameters, then return. % undo the variable transformations into the original space x = xtransform(x0u); % final reshape x = reshape(x,xsize); % stuff fval with the final value fval = feval(fun,x,varargin); % fminsearchbnd was not called exitflag = 0; output.iterations = 0; output.funcount = 1; output.algorithm = 'no call (all variables fixed)'; output.message = 'All variables were held fixed by the applied bounds'; % return with no call at all to fminsearch return end % Add the wrapper function to the user function right here inline: intrafun = @(x, varargin) fun(xtransform(x), varargin{:}); % Added code: Add wrappers to output function(s) and plot function(s) - you can specify multiple % output and/or print functions if you use a cell array of function handles. if ~isempty(options) % Add a wrapper to the output function(s) % fetch the output function and put it(them) into a cell array: 135 OutputFcn = createCellArrayOfFunctions(optimget(options,'OutputFcn',struct('OutputFcn',[] ),'fast'),'OutputFcn'); for ii = 1:length(OutputFcn) %stop = firstOutputFunction(OutStructure, optimValues, state, varargin) OutputFcn{ii} = @(x, varargin) OutputFcn{ii}(xtransform(x), varargin{:}); end % store the "wrapped" output function back into the options. options = optimset(options, 'OutputFcn', OutputFcn); % Add a wrapper to the plot function(s) % fetch the plot function and put it(them) into a cell array: PlotFcn = createCellArrayOfFunctions(optimget(options,'PlotFcns',struct('PlotFcns',[]), 'fast'),'PlotFcns'); for ii = 1:length(PlotFcn) %stop = firstOutputFunction(OutStructure, optimValues, state, varargin) PlotFcn{ii} = @(x, varargin) PlotFcn{ii}(xtransform(x), varargin{:}); end % store the "wrapped" output function back into the options. options = optimset(options, 'PlotFcns', PlotFcn); % Add a wrapper to the print function(s) end % now we can call fminsearch, but with our own % intra-objective function. [xu,fval,exitflag,output] = fminsearch(intrafun,x0u,options,varargin); output.algorithm = [output.algorithm ' bounded using fminsearchbnd']; % undo the variable transformations into the original space x = xtransform(xu); % final reshape x = reshape(x,xsize); % ====================================== % ========= begin NESTED subfunctions ========= % ====================================== function xtrans = xtransform(x) % converts unconstrained variables into their original domains xtrans = zeros(xsize); %zeros(xLength, 1); % I changed this to make it same dimension as the x in fminsearch % was zeros(1, params.xLength) % k allows some variables to be fixed, thus dropped from the % optimization. k=1; for i = 1:xLength 136 switch BoundClass(i) case 1 % lower bound only xtrans(i) = LB(i) + x(k).^2; k=k+1; case 2 % upper bound only xtrans(i) = UB(i) - x(k).^2; k=k+1; case 3 % lower and upper bounds xtrans(i) = (sin(x(k))+1)/2; xtrans(i) = xtrans(i)*(UB(i) - LB(i)) + LB(i); % just in case of any floating point problems xtrans(i) = max(LB(i),min(UB(i),xtrans(i))); k=k+1; case 4 % fixed variable, bounds are equal, set it at either bound xtrans(i) = LB(i); case 0 % unconstrained variable. xtrans(i) = x(k); k=k+1; end end end % sub function xtransform end end % mainline end function [error] = Objective_Exponential_Function_Fung3(x,varargin) % Bell 2012 Fung model but normalize by max stress global lamC lamZ Tcc Tzz coeffs = x; [YfitC, YfitZ] = CalcYfitsFung(coeffs,lamC,lamZ); error = sum(((YfitC - Tcc)).^2/max(Tcc) + ((YfitZ - Tzz)).^2/max(Tzz)); REFERENCES [1] C. A. Taylor, J. M. Bell, M J Breiding, and L.Xu, "Traumatic Brain Injury-Related Emergency Department Visits, Hospitalizations, and Deaths ," MMWR Surveillance Summaries, vol. 66, no. SS-9, pp.1-16, Mar. 1992. [2] https:// www.cdc.gov/ traumaticbraininjury/ get _the _facts.html [3] A. Jamal, M. Bendeck, and B. Langille, "Structural changes and recovery of function after arterial injury," Arteriosclerosis, Thrombosis, and Vascular Biology, vol. 12, no. 3, pp. 307-317, Mar. 1992. [4] M. Ohkawa, N. Fujiwara, M. Tanabe, H. Takashima, K. Satoh, K. Kojima, K. Irie, Y. Honjo, and S. Nagao, "Cerebral vasospastic vessels: histologic changes after percutaneous transluminal angioplasty," Radiology, vol. 198, no. 1, pp. 179184, Jan. 1996. [5] W. Castenada-Zuniga, A. Formanek, M. Tadavarthy, Z. Vlodaver, J. E. Edwards, C. Zollikofer, and K. Amplatz, "The mechanism of balloon angioplasty," Radiology, vol. 135, no. 3, pp. 565-571, Jun. 1980. [6] G.L. Wolf, R.F. LeVeen, and E.J. Ring, "Potential mechanisms of angioplasty," Cardiovascular and Interventional and Radiology, vol. 7, no. 1, pp. 1117, Jan. 1984. [7] P. D. Chan, J.M. Findlay, B. Vollrath, D.A. Cook, M. Grace, M. H. Chen, and R. A. Ashforth, "Pharmacological and morphological effects of in vitro transluminal balloon angioplasty on normal and vasospastic canine basilar arteries," Journal of Neurosurgery, vol. 83, no.3, pp. 522-530, Sep. 1995. [8] Löwenhielm, P., "Dynamic properties of the parasagittal bridging veins," Z Rechtsmedizin, vol. 74, no. 1, pp. 55-62, Mar 1974. [9] Lee, M. C., and R. C. Haut, "Insensitivity of tensile failure properties of human bridging veins to strain rate: implications in biomechanics of subdural hematoma," Journal of Biomechanics, vol. 22, no. 6/7, pp. 537-542, Dec. 1989. [10] D. F Meaney, "Biomechanics of acute subdural hematoma in the subhuman primate and man,"Ph.D dissertation, Dept. Bio Eng., University of Pennsylvania, Philadelphia, PA, 1991. 138 [11] H. Delye, J. Goffin, P. Verschueren, J. Vander Sloten, G. Van der Perre, H. Alaerts, I. Verpoest, and D. Berckmans, "Biomechanical properties of the superior sagittal sinus-bridging vein complex," Stapp Car Crash Journal, vol. 50, pp. 625-636, Nov. 2006. [12] A. G. Monea, K. Baeck, E. Verbeken, I. Verpoest, J. V. Sloten, J. Goffin, and B. Depreitere, "The biomechanical behaviour of the bridging vein-superior sagittal sinus complex with implications for the mechanopathology of acute subdural haematoma," Journal of the Mechanical Behavior of Biomedical Materials, vol. 32, no. c, pp. 155-165, Apr. 2014. [13] J. Chalupnik, C. Daly, and H. Merchant, "Material Properties of Cerebral Blood Vessels," University of Washington, Seattle, WA, USA, Rep No. ME 71-11; 1971. [14] K. L. Monson, W. Goldsmith, N. M. Barbaro, and G. T. Manley, "Axial mechanical properties of fresh human cerebral blood vessels," Journal of Biomechanical Engineering, vol. 125, no. 2, pp. 288-294, Apr. 2003. [15] K. L. Monson, W. Goldsmith, N. M. Barbaro, and G. T. Manley, "Significance of source and size in the mechanical response of human cerebral blood vessels," Journal of Biomechanics, vol. 38, no. 4, pp. 737-744, Apr. 2005. [16] D. Mohan and J. Melvin, "Failure properties of passive human aortic tissue I-uniaxial tension tests," Journal of Biomechanics, vol. 15, no. 11, pp. 887-902, 1982. [17] D. Mohan and J. Melvin, "Failure properties of passive human aortic tissue, II-biaxial tension tests," Journal of Biomechanics, vol. 16, no. 1, pp. 31-44, 1983. [18] Xue Xiao, Na-na Ping, Sen Li, Lei Cao, and Yong-xia Cao, "An optimal initial tension for rat basilar artery in wire myography," Microvascular Research, vol. 97, pp. 156- 158, Jan. 2015. [19] Elizabeth C. Raff, Ann Janice Brothers, and Rudolf A. Raff, "Mechanical properties of vascular smooth muscle cells in situ," Nature, vol. 260, pp. 617-619, Apr. 1976. [20] M. Kowalik and J. Pyrzanowska, "Determination of Mechanical Properties of rat aorta using ring-shaped specimen," Solid State Phenomena, vol. 240, pp. 255-260, 2016. [21] M. Stoiber , B. Messner, C. Grasl, and V. Gschlad, "A method for mechanical characterization of small blood vessels and vascular grafts," Experimental Mechanics, vol. 55, pp. 1591-1595, Oct. 2015. [22] Lakeesha E. Bridges, Cicely L. Williams, Mildred A. Pointer, and Emmanuel M. Awumey, "Mesenteric artery contraction and relaxation studies using automated wire myography," Journal of Visualized Experiments, Issue 55, pp. 3119-3123, Sept. 2011. 139 [23] W. Halpern, M. J Mulvany, and D. M. Warshaw, "Mechanical properties of smooth muscle cells in the walls of arterial resistance vessels," The Journal of Physiology, vol. 275, pp. 85-105, 1977. [24] S. A. Maas, B. J. Ellis, G. A. Ateshian, and J. A. Weiss. (2012). FEBio: Finite Elements for Biomechanics. Journal of Biomechanical Engineering. [Online]. Available: http://help.mrl.sci.utah.edu/help/index.jsp. [25] Fung YC, "What are the residual stresses doing in our blood vessels?" Annals Biomedical Engineering, vol. 19, no. 3, pp. 237-49, May 1991. [26] K. L. Monson, N. M. Brbaro, and G. T. Manley, "Biaxial response of passive human cerebral arteries," Annals of Biomedical Engineering, vol. 36, no. 12, pp. 2028-2041, Dec. 2008. [27] Y C Fung and S Q Liu, "Change of residual strains in arteries due to hypertrophy caused by aortic constriction," Circulation Research, vol. 65, no. 5, pp. 1340-1349, Nov. 1989. [28] LS-DYNA Aerospace working group modeling guidelines document, version 16-1 dated 30,2016, [Online]. Available: http://www.predictiveengineering.com/sites/default/files/awg_lsdyna_modeling_guidelings_document_v13-1.pdf. [29] LS-Dyna Keyword user's manual, volume 1, 05/24/16 (r:7643), LS-Dyna Dev, Livermore Software Technology Corporation (LSTC). [30] LS-DYNA Keyword user's manual, vol. 2, Material models, 07/22/16 (r:7782), LSDYNA Dev, Livermore Software Technology Corporation (LSTC). [31] J. A. Weiss and K. M. Quapp, "Material characterization of human medial collateral ligament," Journal of Biomechanical Engineering, vol. 120, no. 6, pp. 757-63, 1998. [32] https:// water.usgs.gov/ edu/ density.html [33] Herbert E, Balibar S, Caupin F, "Cavitation pressure in water," 2006 [Online]. https://pdfs.semanticscholar.org/0e50/5d19f5499a94981a61602d28992828ff95e2.p df [34] http:// www.engineersedge.com/ physics/ water__ density_ viscosity_ specific_ weight_ 13146.htm [35] http:// www.cvphysiology.com/ Hemodynamics/ H011 Reproduced with permission of copyright owner. Further reproduction prohibited without permission.