4 Finite element based approach for: Fluid-deformable structure interactionssolidt = t1solidThe method we want to develop has to be efficient. For a situation where the fluid and solid interacts, typically, the interface between the fluid and the solid domains needs to be traced and the mesh needs to be updated accordingly based on the actual dynamics of the solid or the fluid. As we all know, mesh updating can be a very time consuming procedure, and it’s not a viable approach especially when there is a relatively large movement or deformation involved in the process. In our IFEM method, we will de-couple the fluid and solid into two separate domains where the two sets of governing equations are solved independently. By doing it, there’s certainly complications involved, but we’ll discuss in just a minute how the interface is handled between the two domains.In the IFEM method, there are two assumptions must be made before we proceed further on developing the equations. 1. No-slip boundary condition is assumed at the fluid-solid interface; 2. the solid must remain immersed in the fluid domain completely the whole timet=0Assumptions:No-slip boundary condition at the fluid-solid interfaceSolid is completely immersed in the fluidFluid is everywhere in the domain

6 Principle of virtual work:Equations of motionPrinciple of virtual work:123s123

7 Solid: in s fluid: in  Overlapping s In our IFEM approach, the solid actually occupies a volume in the system. Therefore, that creates an overlapping domain when we assume fluid is everywhere in the domain.Using the Principle of virtual work. We come up with two sets of independent equations, one for solid and one for fluid. The derivations for these two equations are rather complicated. We’ll omit discussing on the details. However, notice that in the solid equation, the effects from the artificial fluid in the overlapping domain are taken care of by subtracting its influence in all force terms: the inertial, internal, and external forces.How do we link these two domain together? In other words, how will the fluid monitor the solid movement and react to the movement or vice versa? We will consider the following two terms: the velocity and the interacting force. By using the first assumption we made which says that there is no-slip at the interface, we can interpolate the fluid velocity onto the solid domain because the velocity at the overlapping domain should be the same. Another variable we’ll consider is the interaction force which is calculated in the solid domain. To understand this variable in a more physical sense, you can think of it as the force that is generated by the artificial fluid which is governed by the constitutive law for the solid. This force is been interpolated from the solid to the fluid and will be applied as an external force onto the fluid.

8 Interpolations at the interfaceForce distributionVelocity interpolationInfluence domainSurrounding fluid nodessolid nodeUniform spacingAs we said earlier, we have two independent domains for the solid and the fluid. As the solid moves, it’s almost impossible to keep all the nodes overlap on top of each other. If we do that, then we lose the advantage of having two independent meshes. So, how do we do the interpolations? We need to define the neighbors for every solid domain. The interpolation process should not lose any energy. Therefore, the interpolation function must satisfy a very important “consistency condition”, at least to the first order.

19 Use numerical methods to understand and study cardiovascular diseases.Find non-invasive means to predict physical behaviors and seek remedies for diseasesSimulate the responses of blood flow (pressure and velocities) under different physiologic conditions.Compare our results (qualitatively) with published clinical data and analyze the results.

21 Why heart?Cardiovascular diseases are one of the leading causes of death in the western world.Cardiovascular diseases (CVD) accounted for 38.0 percent of all deaths or 1 of every 2.6 deaths in the United States in It accounts for nearly 25% of the deaths in the word.In 2005 the estimated direct and indirect cost of CVD is$393.5 billion.

22 Cardiovascular systemA: The oxygen-rich blood (red) from the pulmonary vein fills the left atrium.B: The oxygen-rich blood in the left atrium fills the left ventricle via the mitra valve.C: The left ventricle contracts and sends the oxygen-rich blood via aortic valve and aorta to the systemic circulation.FAD: The oxygen-poor blood (blue) from the superior vena cava and inferior vena cava fills the right atrium.E: The oxygen-poor blood in the right atrium fills the right ventricle via tricuspid valve.F: The right ventricle contracts and sends the oxygen-poor blood via pulmonary valve and pulmonary artery to the pulmonary circulation.CDBE

23 Atrial fibrillation and blood flowDuring Atrial Fibrillation (a particular form of an irregular or abnormal heartbeat):The left atrium does not contract effectively and is not able to empty efficiently.Sluggish blood flow may come inside the atrium.Blood clots may form inside the atrium. Blood clots may break upResult in embolism.Result in stroke.Without blood clotswith a blood clotLeft atrial appendage

25 Left atrium geometry Left atrium Left atrial appendage Pulmonary veinsFrom Schwartzman D., Lacomis J., and Wigginton W.G., Characterization of left atrium and distal pulmonary vein morphology using multidimensional computed tomography. Journal of the American College of Cardiology, (8): pThe detailed characterization of left atrium and distal pulmonary vein are published by Schwartzman in The geometry includes the 4 pulmonary veins entering the atrium.It is only 10 years ago, that physicians used modern visualization device such as TTE and TEE and found that beside atrium and ventricle that we have known of, there is something identified as the left atrium appendage.The left atrium appendage is a small tube attached to the atrium wall. It may or may not be present depending of the patients. It has been shown that the presence of the appendage is more likely to reduce the blood flow and then increase the probability of thrombus formation. Different geometrical shapes and sizes of the appendage have important impact on flow properties. In 50% of the cases of AF, there are two lobes in the LAA and in 25% there are multi-lobes. The geometrical data for appendage is collected by Ernst in his 1995 paper on “morphology of the left atrial appendage”.We used the geometrical data and constructed a simplified geometry for the left atrium and its appendage. will conduct the simulations with and without the appendage in order to show its influence.Label dimensions for the veins, mitrial valve, and the atriumErnst G., et al., Morphology of the left atrial appendage. The Anatomical Record, : pLeft atrium77mm28mm20mm17mm56mmLeft atrial appendagePulmonary veins

26 Left atrium with pulmonary veinsThe pulmonary venous flow is prescribed at four locations on the left atrial wall where the pulmonary veins enter the chamber. This pulmonary venous flow is given and is coming from experimental data collected by Klein and TajikThe physilogical Reynolds number is about 2400 inside the pulmonary veins and 1000 inside the left atriumThe mitral valve is situated at the bottom of the atrium. In our model, an actual opening and closing mitral valve is not incorporated. Instead, the valve is modeled as a fluid boundary condition. During systole (ventricle systole), the atrium is filling of blood, then there is no flow through the mitral valve (velocities are set to zero) and during diastole (ventricle diastole) blood is ejected from the atrium through the mitral valve to the ventricle. Thus, blood flow is allowed through the mitral valve (free boundary condition).Modeling the valve in this manner provides a means for directing the flow into the left ventricle similar to the physiologic case, but does not account for leaflet flutterDuring diastole (relaxes, 0.06s < t < 0.43s) , no flow through the mitral valve (v=0)During systole (contracts, 0.43s < t < 1.06s), blood flow is allowed through the mitral valve (free flow)Blood is assumed to be Newtonian fluid, homogenous and incompressible.Maximum inlet velocity: 45 cm/sBlood density: 1055 kg/m3Blood viscosity: 3.5X10-3 N/s.m2Fluid mesh: 28,212Nodes, 163,662 ElementsSolid mesh: 12,292 Nodes, 36,427 ElementsKlein AL and Tajik AJ. Doppler assessment of pulmonary venous flow in healthy subjects and in patients with heart disease. Journal of the American Society of Echocardiography, 1991, Vol.4, pp

27 Wall muscle constitutive equationPassive strain during diastoleStrain energyActive strain during systoleThis is done by implementing appropriate constitutive equations to closely mimic the behaviors of biomaterials.Material models have only recently been developed for cardiac muscle tissue. A major challenge in modeling muscle tissue is that in addition to passive nonlinear properties, muscle tissue can generate force, termed activation force. To determine the overall or effective behavior of the myocardium, Xie and Perucchio assumed a strain energy function for the myocardial microstructure. This strain energy is divided in two parts, a passive (Wp) during diastole and an active (Wa) one during systole. The first term is similar to other soft tissues and represents the passive nonlinear properties of the myocardium muscle. The second term is new and accounts for the fact that muscle fibers can generate active stress. εij are components of the Green-Lagrange strain tensor. (using the first and second Piola Kirchooff, we can compute the stress along the muscle).WpWaa1=a2=a3=a4=a5=a6=a7=a8=a9=a10=a11=a12=a13=Green-Lagrange strainSecond Piola-Kirchhoff stressFirst Piola-Kirchhoff stressFrom W. Xie and R. Perucchio, “Computational procedures for the mechanical modeling of trabeculated embryonic myocardium”, Bioengineering Conference, ASME 2001, BED-Vol. 50, pp

28 Left atrium with appendagePressure distribution at the center of the atrium during a diastole and systole cycleWe implemented this constitutive equation and applied displacement along the walls to simulate the beating of the left atrium.The volume of the atrium at end of systole is around 2.5 times the volume of the atrium at the end of diastole. This simulates the contraction and relaxation of the heart muscles. Thus, we get this following regular left atrium heart beat.Here is the pattern of the blood flow profile inside the atrium during the cardiac cycles for the case of moving boundaries. As for the rigid case, on this movie, we plot the pressure distribution (top left plot) at the center of the atrium, the transmitral flow velocity (bottom left plot) and the blood pattern inside the atrium during one complete cardiac cycle. We can clearly identify on this movie both systole (filling of the atrium) and diastole (ejection of the blood). It must be pointed out that during diastole, the filling of the atrium occurs also from the pulmonary veins.Transmitral velocity during diastole

30 Left atrium (comparison with clinical data)Kuecherer H.F., Muhiudeen I.A., Kusumoto F.M., Lee E., Moulinier L.E., Cahalan M.K. and Schiller N.B., Estimation of mean left atrial pressure from transesophageal pulsed Doppler echocardiography of pulmonary venous flowCirculation, 1990, Vol 82,5Pressure (mm hg)2Time (s)1.51Pressure distribution at the center of the atrium during one cardiac cycleAs you can see, the pattern of the pressure distribution is realistic compared to clinical data. Same conclusion can be drawn for the transmitral velocity. Both E and A waves of the mitral valve during diastole are clearly identified.First, systole occurs. This is the filling of the atrium from the pulmonary veins. At this time, the mitral valve is closed (till 1.4s). The pressure increases quickly at the beginning of systole then decreases.Then, diastole starts at time 1.4s. Blood quickly ejects from the atrium and fills the ventricle, pressure declines sharply. This phase is known as the rapid ejection phase and is very short (around 100 ms) (E wave). Then, there is the reduce ejection phase. Flow across the mitral valve is greatly diminished. Then, pressure and velocity increases for 100 ms This is the atrial contraction phase. It results in the second filling of the ventricle (A wave).Now, we can introduce the appendage in our simulations in order to study its influence on the blood flow and blood clots formation.AETransmitral velocity during one cardiac cycle

31 Influence of the appendageIn order to show the influence of the appendage, we plot on the same graph the transmitral velocity with and without the appendage. With the appendage, the stroke volume decreases. The flow through the mitral valve reduces.The next plot shows the velocity profile of the blood inside the appendage. This velocity is very small compared to the pulmonary venous flow at the inlet and the transmitral flow (at the exit) of the atrium. Thus, there is a low shear rate at this location. Since, it has been proven experimentally that a low shear rate will increase the probability of thrombus formation; our results confirm that the appendage is an attractive location for thrombus formation. And as a consequence, in the case of atrial fibrillation, engenders stroke or embolism due to the breaking of a clot.Plot of a shear stress contour for the appendage????Transmitral velocity during one cardiac cycle (with and without the appendage)Velocity inside the appendage during one cardiac cycle

32 Influence of the appendageIn order to show the influence of the appendage, we plot on the same graph the transmitral velocity with and without the appendage. With the appendage, the stroke volume decreases. The flow through the mitral valve reduces.The next plot shows the velocity profile of the blood inside the appendage. This velocity is very small compared to the pulmonary venous flow at the inlet and the transmitral flow (at the exit) of the atrium. Thus, there is a low shear rate at this location. Since, it has been proven experimentally that a low shear rate will increase the probability of thrombus formation; our results confirm that the appendage is an attractive location for thrombus formation. And as a consequence, in the case of atrial fibrillation, engenders stroke or embolism due to the breaking of a clot.Plot of a shear stress contour for the appendage????Transmitral velocity during one cardiac cycle (with and without the appendage)Velocity inside the appendage during one cardiac cycle

33 Then what? Use realistic atrial geometry How?Medical School (Computed Tomography CT scan), but the device is ruined due to KatrinaHelp from Dr. A. Cristoforetti, University of Trento, ItalyA Computed Tomography (CT) scan shows the body’s organs in greater detail and more clearly than regular x-rays. CT uses x-rays and a computer to make a picture of sections of the body. “CTA” stands for computed tomography angiography. A CTA scan gives a view of specific blood vessels (arteries and veins).Computed Tomography is probably the most common source of three-dimensional data. CT scanners are relatively inexpensive. CT uses an X-ray radiation source to image the patient. The CT scanner consists of a couch upon which the patient is placed and a circular gantry through which the couch with patient is passed. Within the gantry is a rotating ring with an X-ray source opposed to a linear array of detectors. The X-ray source is collimated so that the X-rays form a flat fan beam with a thickness determined by the user. During the acquisition of a "slice" of data, the source-detector ring is rotated around the patient. The raw output from the detector array is back projected to form an image of the slice of the body. The couch is moved and then another slice is obtained. Computed tomography (CT) was the first non-invasive radiological method allowing the generation of tomographic images of every part of the human body without superimposition of adjacent structures.The output from a CT scanner is a series of transaxial slices of the patient. Each slice represents a slab of the patients' body with a thickness set by the collimation for the slice (typically 1-10mm). For most CT scanners each slab has 512 by 512 pixels. The size of a pixel can be varied within certain limits (generally 0.5 to 2 mm). Generally each slice is spaced such that they are either overlapping or contiguous, though some protocols call for gaps between the slices. Each pixel ideally represents the absorption characteristics of the small volume within its bounds.

34 atrial volume Atrial systole Atrial diastole 1 2 3 4 5 6 7 ECGECGLA VolumeMin LA volumeMax LA volume75% of max LA volume1 Atrial contraction2 Isovolumetric contraction3 Rapid ejection4 Reduced ejection5 Isovolumetric relaxation6 Rapid ventricular filling7 Reduced ventricular fillingThe period in which the heart has the least motion is usually (but not always) in diastole, near a phase between 55% and 75%.The maximum volume of the left atrium is at end atrial diastole. The minimum volume is at end atrial systole.

41 MAVs Loitering wings: high span and a large surface areaFast wings: a low wing span and a small areaFlying efficiently at high speed: small, perhaps, swept wingsFlying at slow speed for long periods: long narrow wingsFeatures:improved efficiency,more lift,high maneuverability,reduced noise.

50 IFEM: Fluid solving algorithmPetrov-Galerkin Weak Form and discretizationWith τm and τc as stabilization parameters depending on the grid sizeNewton Iteration: solve for the 4 unknowns per node: u, v, w, p (three velocity components + pressure)Matrix-free formulation is solved by the Generalized Minimum Residual Method (GMRES)Note that the force exerted from the structure is not updated during the Newton Iteration, therefore the coupling is explicit.

54 Red blood cell model RBCShear rate dependence of normal human blood viscoelasticity at 2 Hz and 22 °C (reproduced fromBulk aggregatesDiscrete cellsCell layersRBCThe density of RBC is about 15% higher than the fluid (density of blood in 45% of hematocrit is 1.07 g/ml)Thickness of RBC: 7.5 to 10nmHematocrite: the percent volume of red blood cellsNeo-Hookean or Mooney-Ravenlin material is used.The shape of blood cell is described by a empirical function from experimental measurement of RBC cellFrom Dennis Kunkel at

66 Multi-resolution analysisWindow function with a dilation parameter:a: dilation parameterProjection operator for the scale aWavelet function:To understand the ability of the RKPM to be used for multiresolution analysis, it is necessary to look at its behavior in the frequency or Fourier domain. The key to this ability is the dilation parameter a of the window function. The kernel function acts like a lowpass filter for the solution in this system. By changing the dilation parameter a, it is possible to construct a series of lowpass filter that provide different frequency parts of the solution. Substracting two low-frequency solution parts that are obtained for two different dilations a, a and 2a, we get the high frequency part, which is simply the difference between these two scales.My research was to take advantage of all the benefits that you get from the meshfree method, and to eliminate or reduce the inconvenience that are caused by the meshfree method. By accomplishing that, researchers can use the meshfree approach carefree without have to worry about boundary treatment or large amount of time they use for the calculation.Complementary projection operator:low scale + high scale