Délivré par l’Université de Montpellier

Préparée au sein de l’école doctorale I2S Et de l’unité de recherche IMAG Spécialité: Mathématiques et Modélisation Présentée par Étienne Gibaud

Numerical simulation of red blood cells flowing in a blood analyzer

Soutenue le 15 décembre 2015 devant le jury composé de Franck Nicoud

IMAG - Univ. de Montpellier

Directeur

Simon Mendez

IMAG CNRS - Univ. Montpellier

Co-encadrant

Anne-Virginie Salsac

BMBI CNRS - Univ. Tech. Compiègne

Rapporteur

Marc Jaeger

M2P2 Centrale Marseille

Rapporteur

Chaouqi Misbah

LIPhy - Univ. Joseph Fourier Grenoble

Examinateur

Vincent Moureau

CORIA CNRS - Rouen

Examinateur

Damien Isèbe

Horiba Medical SAS - Montpellier

Invité

Abstract The aim of this thesis is to improve the understanding of the phenomena involved in the measurement performed in a blood analyzer, namely the counting and sizing of red blood cells based on the Coulter effect. Numerical simulations are performed to predict the dynamics of red blood cells in the measurement regions, and to reproduce the associated electrical measurement used to count and size the cells. These numerical simulations are performed in industrial configurations using a numerical tool developed at IMAG, the YALES2BIO solver. Using the Front-Tracking Immersed Boundary Method, a deformable particle model for the red blood cell is introduced which takes the viscosity contrast as well as the mechanical effects of the curvature and elasticity on the membrane into account. The solver is validated against several test cases spreading over a large range of regimes and physical effects. The velocity field in the blood analyzer geometry is found to consist of an intense axial velocity gradient in the direction of the flow, resulting in a extensional flow at the micro-orifice, where the measurement is performed. The dynamics of the red blood cells is studied with numerical simulations with different initial conditions, such as its position or orientation. They are found to reorient along the main axis of the blood analyzer in all cases. In order to understand the phenomenon, the Keller & Skalak model is adapted to the case of extensional flows and is found to reproduce the observed trends. This thesis also presents the reproduction of the electrical measurement used to count red blood cells and measure their volume distribution. Numerous dynamics simulations are performed and used to generate the electrical pulse corresponding to the passage of a red blood cell inside the micro-orifice. The resulting electrical pulse amplitudes are used to characterize the electrical response depending on the initial parameters of the simulation by means of a statistical approach. A Monte-Carlo analysis is performed to quantify the errors on the measurement of cell volume depending on its orientation and position inside the micro-orifice. This allows to construct the volume distribution of a well defined population of red blood cells and the characterization of the associated measurement errors. Key words: Red Blood Cells, Blood Analyzer, Immersed Boundary Method, FluidStructure Interaction, Electrostatics.

iii

R´ esum´ e L’objectif de cette th`ese est d’am´eliorer la compr´ehension des ph´enom`enes jouant un rˆ ole dans la mesure effectu´ee dans un analyseur sanguin, en particulier le comptage et la volum´etrie d’une population de globules rouges reposant sur l’effet Coulter. Des simulations num´eriques sont effectu´ees dans le but de pr´edire la dynamique des globules rouges dans les zones de mesure et pour reproduire la mesure ´electrique associ´ee, servant au comptage et ` a la volum´etrie des cellules. Ces simulations sont effectu´ees `a l’int´erieur de configurations industrielles d’analyseurs sanguins, en utilisant un outil num´erique d´evelop´e ` a l’IMAG, le solveur YALES2BIO. En utilisant la m´ethode des fronti`eres immerg´ees avec suivi de front, un mod`ele de particule d´eformable prenant en compte le contraste de viscosit´e ainsi que les effets m´ecaniques de la courbure et de l’´elasticit´e sur la membrane est introduit. Le solveur est valid´e grˆace `a de nombreux cas tests pour ´evaluer diff´erents r´egimes et effets physiques. L’´ecoulement fluide dans cette g´eom´etrie d’analyseur sanguin est caract´eris´e par un fort gradient de vitesse axial dans la direction de l’´ecoulement, impliquant la pr´esence d’un ´ecoulement extensionel au niveau du micro-orifice, l` a o` u a lieu la mesure. La dynamique des globules rouges est ´etudi´ee par des simulations num´eriques pour diff´erentes conditions initiales, telles que sa position ou son orientation. Il est observ´e que les globules rouges se r´eorientent selon l’axe principal de l’analyseur sanguin dans tous les cas. Pour comprendre le ph´enom`ene, le mod`ele de Keller & Skalak est adapt´e au cas des ´ecoulements extensionels et reproduit correctement les tendances de r´eorientation. Cette th`ese pr´esente ´egalement la reproduction de la mesure ´electrique utilis´ee pour le comptage et la mesure de la distribution des volumes de globules rouges. De nombreuses simulations de la dynamique des globules rouges sont effectu´ees et utilis´ees pour g´en´erer l’impulsion ´electrique correspondant au passage du globule rouge dans le micro-orifice. Les amplitudes d’impulsions ´electriques r´esultantes permettent la caract´erisation de la r´eponse ´electrique en fonction des param`etres initiaux de la simulation par une approche statistique. Une analyse par simulation Monte-Carlo est utilis´ee pour la quantification des erreurs de mesure li´ees `a l’orientation et la position des globules rouges dans le micro-orifice. Ceci permet la g´en´eration d’une distribution de volume mesur´ee pour une population de globules rouges bien d´efinie et la caract´erisation des erreurs de mesure associ´ees. Mots clefs: Globules Rouges, Analyseur Sanguin, M´ethode des Fronti`eres Immerg´ees, ´ Int´eraction Fluide-Structure, Electrostatique.

v

”Le vide et le silence sont toujours en nous, ici et maintenant, dans le meilleur et dans le pire, sous le tumulte de nos chants d’amour et de nos cris de haine, de nos peurs et de nos pleurs, de nos rˆeves et de nos joies, de nos ambitions et de nos amertumes, de nos triomphes et de nos d´esarrois.” — Arnaud Desjardins

` mon p`ere... A

Acknowledgements Remerciements

Je souhaite tout d’abord remercier les organismes qui ont rendu cette th`ese possible, ` a travers un financement, une collaboration ou des moyens techniques. Merci ` a BPIFrance, au Labex Numev, Horibal Medical et les membres du projet DAT@DIAG, l’Agence Nationale de la Recherche, le CNRS, le CINES et ´evidemment le laboratoire I3M/IMAG et l’Universit´e de Montpellier. Je remercie ´egalement les membres du jury de ma th`ese: Anne-Virginie Salsac, Marc Jaeger, Chaouqi Misbah, Vincent Moureau et Damien Is`ebe pour leur disponibilit´e et leur professionnalisme. Je suis extrˆemement honor´e par la constitution de ce jury, qui r´eunit chercheurs et professionnels de haut niveau, que j’ai pu admirer tout au long de mes nombreuses recherches bibliographiques ou avec qui j’ai eu la chance de collaborer. Je souhaite remercier chaleureusement mes deux encadrants de th`ese, Franck Nicoud et Simon Mendez. Vous m’avez permis au fil des ann´ees d’aiguiser mes r´eflexions et d’affiner ma fa¸con d’appr´ehender le travail. Je vous remercie pour votre pr´esence permanente, votre intransigeance et votre comp´etence. Je suis heureux d’avoir pu effectuer ce travail avec vous, cel`a m’a permis de donner le meilleur de moi-mˆeme et d’aller au bout de ce travail ´eprouvant. Je vous remercie ´egalement de votre compr´ehension face aux probl`emes que j’ai pu rencontrer durant ma th`ese, qu’ils soient d’ordre personnel ou professionnel. J’esp`ere que le futur nous am`enera `a des opportunit´es de travail ensemble, car je pense en avoir encore beaucoup `a apprendre `a votre contact. Je remercie ´egalement les diff´erents collaborateurs qui ont contribu´e `a l’aboutissement de ce travail: Christophe Chnafa, Julien Sig¨ uenza, Julien Stoehr, Myriam Tami, Marco Martins Afonso, Manouk Abkarian, Damien Is`ebe, Daniele Di Pietro et Xavier Bry. Merci pour vos id´ees, nos d´ebats et toutes les discussions qui m’ont amen´e `a pousser la r´eflexion au maximum et ` a m’aventurer dans des domaines que je ne connaissais pas ou peu. Je souhaite ´egalement remercier mes co-bureaux principaux, ceux avec qui j’ai partag´e mon quotidien professionnel et personnel. Merci `a Christophe Chnafa, maˆıtre inconstest´e de l’investissement, qui m’a toujours impressionn´e par sa motivation mais aussi par ses intuitions et ses fulgurances. Je remercie aussi Julien Sig¨ uenza, ing´enieur et chercheur hors pair, qui m’a beaucoup appris de par sa rigueur et son acharnement. Merci pour les bons moments pass´es en votre compagnie, votre acceptation de mes travers et d´efauts, et pour votre disponibilit´e sans faille. Je n’aurais pas pu esp´erer ix

x meilleurs compagnons de gal`ere, et je vous suis reconnaissant d’avoir partag´e un bout de votre vie avec moi. Merci ´egalement `a Angelina Roche, Jonathan Ohayon et les autres, vous m’avez d´emontr´e votre ouverture d’esprit et votre force mentale en supportant mes nombreuses (et aga¸cantes) chansons improvis´ees et le non-sens dont je sais faire preuve en fin de journ´ee. Je remercie ´egalement l’ensemble des membres du laboratoire pour leur pr´esence et nos int´eractions. C’est grˆace `a vous que j’ai pu conserver ma sant´e mentale et sortir vivant de ces trois derni`eres ann´ees de labeur. Merci au personnel administratif, aux chercheurs, permanents et non-permanents, aux doctorants et aux stagiaires, vous m’avez appris un nombre incalculable de choses et m’avez soutenu sans repˆıt. Je remercie tout particuli`erement ma m`ere et mes soeurs, L´ea et Marie, pour leur soutien crucial au fil des mes ´etudes. Votre pr´esence m’a aid´e `a supporter les p´eriodes difficiles que j’ai pu traverser. Merci pour les longues nuits sans dormir, les instants r´egressifs et l’acceuil que vous m’avez apport´e. Merci `a ma m`ere pour son soutien sans faille, ses comp´etences hors normes et le refuge que sa maison `a constitu´e pour moi au fil de ces ann´ees. Un remerciement particulier ` a Julie, Xavier et Manu pour leur pr´esence et leurs conseils, sans qui je ne serais pas o` u j’en suis aujourd’hui. Je remercie aussi toute la bande: Serge, Joris, Th´eo, Jona, Guillaume et Yohan en particulier, mais aussi tous les autres, vous m’avez permis (entre autres) de me prouver que j’´etais capable d’arriver au bout de mes projets. Je souhaite remercier tout particuli`erement mon p`ere, qui `a su me pousser `a aller au bout de mes ambitions, m’` a amener `a remettre en question mes perspectives, ma fa¸con de travailler et qui m’a toujours soutenu, mˆeme bien avant le compte `a rebours de seizes semestres qu’ont ´et´e mes ´etudes `a l’Universit´e de Montpellier. Merci pour ta pr´esence et ton soutien, je te serais ´eternellement reconnaissant de l’h´eritage intellectuel et spirituel que tu m’as transmis. Je remercie mon amour, ma compagne, ma femme: Mathilde. Merci pour ta pr´esence `a mes cˆot´es pendant ces derni`eres ann´ees, ton soutien sans relˆache et nos int´eractions m’ont vraiment permis d’aller au bout de cette th`ese et de me relever apr`es les coups durs. Merci de m’avoir support´e, d’avoir continu´e `a croire en moi, d’avoir ´et´e patiente et indulgente. Tu es la personne la plus importante pour moi, et je suis persuad´e que notre avenir sera radieux : les dieux sont avec nous apr`es tout...

Montpellier, le 15 d´ecembre 2015

E. G.

Contents Abstract (English/Franc ¸ais)

iii

Acknowledgements

ix

Contents

xii

Chapter 1 Introduction 1.1 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Horiba & Horiba Medical . . . . . . . . . . . . . . . . . . . . . . 1.2 Blood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Role . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Haemopathies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Blood analyzers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Coulter effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Hydrofocalization . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Numerical simulation of red blood cells and blood analyzers . . . . . . 1.4.1 Numerical simulation of the red blood cell dynamics . . . . . . . 1.4.2 Literature review: a historical perspective . . . . . . . . . . . . . 1.4.3 Description of existing numerical approaches for deformable particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.4 Simulation of red blood cells: state of the art . . . . . . . . . . . 1.4.5 Configurations of interest in numerical studies of flows of deformable particles . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 2 2 3 3 3 6 8 8 10 11 12 12 14

Chapter 2 Numerical method 2.1 YALES2 & YALES2BIO . . . . . . . . . . . . . . . . . . . . 2.1.1 YALES2 . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 YALES2BIO . . . . . . . . . . . . . . . . . . . . . . . 2.2 Modelling framework . . . . . . . . . . . . . . . . . . . . . . 2.3 Fluid-structure coupling with the immersed boundary method 2.4 Discretization of the coupled problem . . . . . . . . . . . . . 2.4.1 Hyperelastic forces . . . . . . . . . . . . . . . . . . . 2.4.2 Curvature forces . . . . . . . . . . . . . . . . . . . . . 2.4.3 Resolution . . . . . . . . . . . . . . . . . . . . . . . .

23 24 24 24 24 26 27 28 30 32

xii

. . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

15 17 20 21

xiii

CONTENTS 2.5 2.6 2.7 2.8

Interpolation and regularization over unstructured meshes Particle volume conservation . . . . . . . . . . . . . . . . YALES2BIO timestep description . . . . . . . . . . . . . . Validation test cases . . . . . . . . . . . . . . . . . . . . . 2.8.1 Pressurized capsules . . . . . . . . . . . . . . . . . 2.8.2 Equilibrium shapes . . . . . . . . . . . . . . . . . . 2.8.3 Optical tweezers . . . . . . . . . . . . . . . . . . . 2.8.4 Spherical capsules in linear shear show . . . . . . . 2.8.5 Red blood cells in simple shear flow . . . . . . . .

. . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

33 37 38 39 39 41 44 47 61

Chapter 3 Analytical modeling of the dynamics of red blood elongational flow 3.1 Derivation of the Keller & Skalak model . . . . . . . . . . . . 3.1.1 Formulation of the problem . . . . . . . . . . . . . . . . 3.1.2 External flow and equilibrium . . . . . . . . . . . . . . . 3.1.3 Internal flow and conservation of energy . . . . . . . . . 3.1.4 Abkarian, Faivre & Viallat’s variation . . . . . . . . . . 3.1.5 Planar elongational flow derivation . . . . . . . . . . . . 3.1.6 Jeffery’s theory . . . . . . . . . . . . . . . . . . . . . . . 3.2 Models behavior in the case of a planar elongational flow . . . 3.2.1 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Effect of the capillary number . . . . . . . . . . . . . . . 3.2.3 Effect of the viscosity ratio . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

67 68 68 70 73 75 77 78 79 79 80 84

Chapter 4 Red blood cells dynamics 4.1 Numerical configuration . . . . . . . . . . 4.1.1 Red blood cell . . . . . . . . . . . 4.1.2 Configuration . . . . . . . . . . . 4.1.3 Reduced configuration . . . . . . 4.1.4 Initial state of the red blood cell . 4.2 Results on the centerline . . . . . . . . . . 4.2.1 Performed cases . . . . . . . . . . 4.2.2 General dynamics . . . . . . . . . 4.2.3 Deformation and elongation . . . . 4.2.4 Comparison with analytical models 4.2.5 Influence of the viscosity ratio . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

87 88 88 89 97 99 102 102 102 106 109 113

. . . .

117 118 118 119 119

Chapter 5 Electrical response 5.1 Red blood cells interaction with 5.1.1 AC case . . . . . . . . . 5.1.2 DC case . . . . . . . . . 5.2 Numerical method . . . . . . .

. . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . .

electrical fields . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

. . . . . . . . . . .

. . . .

cell in

. . . .

. . . .

. . . .

. . . .

xiv

CONTENTS

5.3

5.4

5.5

5.6

5.7

5.2.1 Modelling framework . . . . . . . . . . . . . . . . . . . 5.2.2 Computational treatment . . . . . . . . . . . . . . . . Validation test cases . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Numerical domain . . . . . . . . . . . . . . . . . . . . 5.3.2 Constant conductivity . . . . . . . . . . . . . . . . . . 5.3.3 Variable conductivty . . . . . . . . . . . . . . . . . . . Computation of the electrical potential in the absence of cells 5.4.1 Full configuration . . . . . . . . . . . . . . . . . . . . . 5.4.2 Reduced configuration . . . . . . . . . . . . . . . . . . Electrical pulse generated by the passage of a cell . . . . . . . 5.5.1 Effect of the passage of a cell . . . . . . . . . . . . . . 5.5.2 Procedure of generation of the electrical pulse . . . . . 5.5.3 Range of the study . . . . . . . . . . . . . . . . . . . . 5.5.4 Test case with a rigid particle . . . . . . . . . . . . . . The probability density function of volume measurements . . 5.6.1 Why a statistical approach . . . . . . . . . . . . . . . 5.6.2 Random variables: initial position and orientation . . 5.6.3 Raw results . . . . . . . . . . . . . . . . . . . . . . . 5.6.4 Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6.5 Probability density function . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Chapter 6 Conclusion 6.1 Main results . . . . . . . . . . . . 6.2 Perspectives . . . . . . . . . . . . 6.2.1 Alternative geometries . . 6.2.2 Red blood cell mechanics 6.2.3 Cell-cell interaction . . . . 6.2.4 Electrical modeling . . . . 6.3 Concluding remark . . . . . . . . Bibliography

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . .

119 122 123 124 125 126 128 128 130 133 133 133 134 135 135 135 136 142 145 146 150

. . . . . . .

151 151 153 153 153 154 154 154 155

Chapter

1

Introduction Chapter contents 1.1 1.2

1.3

1.4

1.5

Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Horiba & Horiba Medical . . . . . . . . . . . . . . . . . . . . . . Blood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Role . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Haemopathies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Blood analyzers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Coulter effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Hydrofocalization . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical simulation of red blood cells and blood analyzers . . . . . . 1.4.1 Numerical simulation of the red blood cell dynamics . . . . . . . 1.4.2 Literature review: a historical perspective . . . . . . . . . . . . . 1.4.3 Description of existing numerical approaches for deformable particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.4 Simulation of red blood cells: state of the art . . . . . . . . . . . 1.4.5 Configurations of interest in numerical studies of flows of deformable particles . . . . . . . . . . . . . . . . . . . . . . . . . . Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 2 3 3 3 6 8 8 10 11 12 12 14 15 17 20 21

The motivations behind this work and the ones of Horiba Medical, the company with which we collaborated during this thesis, are first presented. Then, the role and contents of blood as well as a few haemopathies are described. An explanation of the main principles of blood analyzers and their usage are shown. Finally, a summary of the current state of the art of numerical simulations of red blood cells is presented. 1

2

1.1

CHAPTER 1. INTRODUCTION

Motivations

The first motivation for this work originates from Horiba Medical, a company designing blood analyzers and other medical devices. Aside from the obvious technological, industrial and economical motivations from Horiba, the main objective is the understanding of the red blood cell dynamics in these devices, both in terms of pure application in the medical facilities and in terms of scientific progress in the field of deformable particles dynamics in complex flows.

1.1.1

Horiba & Horiba Medical

The Horiba Group is an international company operating worldwide with a very wide spectrum of instruments and systems for various applications (process and environmental monitoring, in-vitro medical diagnoses, semiconductor manufacturing and metrology) as well as a broad range of scientific research and development and quality control measurements. Horiba Medical is Horiba’s segment dedicated to the design, development and distribution of medical devices mainly destined to biological analysis in medical laboratories. Horiba Medical equips numerous laboratories throughout the world, producing approximately 7 500 automated analyzers per year and over 10 000 tons of reagents used to operate the analyzers. Horiba Medical is dedicated to the improvement of the quality of its devices in order to provide the most accurate, reliable and fastest products on the market. The company hires over 1,100 employees in 5 production centers (Japan, China, France, Brazil, India) and 2 R&D centers (Japan, France). Working with around 100 distributors to supply around 30 000 laboratories worldwide, the company sales reached 198 million euros in 2014 and represents 17% of the company sales for all segments. The company is also interested in multiplying the tools available for blood analysis with many innovative devices being in development and invests around 10% of the company sales in R&D. In this context, there is a need for a versatile framework of tools aimed at the design of blood analyzers of all sorts. On the one hand, Horiba Medical has dedicated teams of experimentalists focused on the production of prototypes and on the thorough testing of these units. These teams have been active for many years. On the other hand, dedicated R&D teams work on the prior design of such devices using commercial numerical tools which often show limitations in complex geometries and flows. In particular, commercial numerical tools cannot handle the complexity of the red blood cells dynamics when submitted to an outer flow. This is why a powerful numerical framework would be an important asset to Horiba Medical, as it would let the company master the functioning of their devices, allowing conscious design choices to maximize efficiency and profitability.

1.2. BLOOD

1.2

Blood

1.2.1

Role

3

Blood is a vital body fluid that aims at delivering the necessary substances (nutrients and oxygen) to the organs, and transporting waste products such as carbon dioxide out of the organism to evacuation in the kidneys, the lungs, the intestines or the spleen. It travels over the 100 000 km of the circulatory system, being pumped by the heart through the many ramifications of the vascular system vessels. Most tissues in the body have to be vascularized, otherwise they would die. Blood is a complex substance composed of a fluid, the plasma, in which various particles are suspended. These particles are white blood cells, platelets and mostly red blood cells, representing around 45% of the whole blood volume. These blood cells are created during a process called haematopoiesis, occurring in the bone marrow (and the thymus for some types of white blood cells). At the end of their life span, these blood cells are destroyed by macrophages in the kidney, the spleen or the bone marrow.

1.2.2

Contents

Plasma Plasma is the carrying fluid of blood, where the previously mentioned cells are suspended and transported inside the organism through the vascular network. It contains around 90% of water along with a wide variety of elements, such as mineral salts, dissolved ions, organic molecules (carbohydrates, lipids, proteins) and metabolic wastes like urea, uric acid or bilirubin, as well as hormones. Altogether, these components and the plasmatic proteins contribute to the acid-base homeostasis of blood (pH), its viscosity and the behavior of the blood cells membrane by maintaining the adequate osmotic pressure levels. It constitutes around 45% of the whole blood volume, with a viscosity close to the one of water.

White Blood Cells White blood cells, or leucocytes, are blood cells with a nucleus. They are dedicated to the defense of the immune system against bacterial infections and other various threats. There are many different types of white blood cells, each one being specialized in the treatment of a specific threat (granulocytes, lymphocytes and monocytes). They represent about 0.2% of the blood with 4 000 to 11 000 cells per cubic millimeter of blood. They are usually bigger than red blood cells with a diameter comprised between 7 and 20 µm depending on their type and live up to several days on average.

4

CHAPTER 1. INTRODUCTION

Platelets Platelets, or thrombocytes, are anucleate cells, involved in the coagulation mechanism (hemostatis) to repair ruptures in the vascular endothelium. They are of biconvex discoid shape and are the result of the fragmentation of bigger cells (megakaryocyte) that burst when entering the circulatory system. They are smaller than red and white blood cells, with a maximum diameter of 2 to 3 µm and live up to 10 days. They average count of platelets is 150 000 to 400 000 per cubic millimeter of blood. Red blood cells Red blood cells, or erythrocytes, are anucleate cells constituted by a membrane enclosing an internal fluid, the cytoplasm. Red blood cells travel more than 200 km in the circulatory system during their 120 days life span. The average count of red blood cells is 4 to 6 million cells per cubic millimeter of blood. At rest, red blood cells generally have the shape of a discocyte. Their average dimensions as measured by Evans and Fung [61] are presented on Fig. 1.1. Other measurements are also presented by Fung [80].

Figure 1.1: Average dimensions of a human red blood cells as reported by Evans and Fung [61]. The cytoplasm contains water and haemoglobin, a protein accounting for more than 90% of the red blood cell dry mass. Haemoglobin is the substance that carries and releases oxygen from the respiratory organs to the rest of the body, thanks to its iron core that can bind oxygen molecules. The cytoplasm viscosity is 6 to 8 times greater than the surrounding plasma viscosity. Red blood cells are subjected to high deformations when passing through the multiple vessels of the circulatory system which radius can be less than 3 µm, which is lower than the red blood cells characteristic dimension. Many obstacles are also in the way of red blood cells through the organism: atheroma plaques, clogging of the vessels with lipids and platelets, turbulence in the blood flow, vasoconstriction, spleen slits and the relative cluttering due to the presence of a large amount of cells in blood. In order to withstand these difficulties, the red blood cells are required to be able to deform and stretch, while keeping their integrity. These qualities are ensured by its geometrical nature on the one hand, with its surface being

5

1.2. BLOOD

greater than the one of a sphere with the same volume, resulting in its deflated aspect, and on the other hand with its membrane structure. The red blood cells membrane is a complex structure composed by a lipid bilayer and an underlying protein network, the cytoskeleton; these components are linked via embedded transmembrane proteins (Fig. 1.2). These two structures constituting the membrane are responsible for the complex mechanics associated to the deformation of red blood cells (see Chapter 2 of Lim et al. [119] for more details): • The lipid bilayer (A on Fig. 1.2) gives to the membrane its resistance against area dilation and bending. • The cytoskeleton (B on Fig. 1.2) gives to the membrane its resistance against shear deformations and area dilation (though with a much smaller modulus than the lipid bilayer).

Figure 1.2: Schematic representation of the red blood cells membrane as presented by Mohandas and Evans [136].

Cell and membrane geometry Red cell area (µm2 ) Red cell volume (µm3 ) Scale length (µm) Mechanical moduli

ARBC = 113.8 ± 27.6 VRBC = 89.4 ± 17.6 RRBC = 2.76 ± 0.18

Shear modulus (N.m−1 ) Area dilation modulus (N.m−1 ) Bending modulus (N.m) Surface shear viscosity (N.s/m)

Es = 5.5 ± 3.3 × 10−6 Ea = 0.39 ± 0.11 Eb = 1.14 ± 0.9 × 10−19 µS = 0.515 ± 0.19 × 10−6

Table 1.1: Typical values and dispersion of the red blood cell geometrical and mechanical parameters [119, 184]. A summary of the red blood cells geometrical and mechanical parameters are presented in Table 1.1. These parameters as well as the physiological construction of the

6

CHAPTER 1. INTRODUCTION

membrane ensure the solidity and deformability of the red blood cells all along their travel in the circulatory system. In this thesis, because of the numerical method chosen for the modeling of the red blood cell, the surface shear viscosity (the membrane viscosity) is not taken into account.

1.2.3

Haemopathies

Examples of haemopathies Many diseases, called haemopathies, are related to the blood cells population, either on the numbering itself, their mechanical and geometrical properties or the properties of the blood. Those diseases include: • Leukemia, a group of cancers resulting in an abnormal production of white blood cells. It is the most common type of cancer in children and is still the cause of more than 200, 000 deaths worldwide each year [186]. • Drepanocytosis, or sickle-cell disease, where red blood cells tend to assume a rigid and unusual shape, causes severe infections and strokes. It is identified as a haemoglobinopathy, corresponding to an anomaly in the hemoglobin molecule. It also causes more than 200, 000 deaths worldwide each year often associated with a decrease in haemoglobin. • Myelodysplastic syndrome, a medical condition resulting in a low production of blood cells causing an increased number of infections and subcutaneous hemorrhagings. Notably, signs of anemia (low red blood cells count or reduced haemoglobin) or thrombocytopenia (low platelet count) are witnessed. These diseases, along with others, are still causing an important number of deaths every year. The ability for physicians to perform a precise, fast and reliable diagnosis is thus of first importance. Anemia In several blood related diseases, signs of anemia are witnessed in patients, which can be detected using blood analyzers. Anemia comes in many forms but invariably translates in a reduced circulation of the haemoglobin in the blood, causing an inefficient transport of dioxygen in the circulatory system. A patient is considered anemic when his/her haemoglobin rate is under a threshold value, which depends on the type of population considered: adult male (13g/100ml), adult female/children (12g/100ml) or newborn (14g/100ml). The three main types of anemia are described in the following figures, where the haematocrit is the volumetric fraction of red blood cells in blood. • Microcytic anemia: the red blood cell count is normal but the mean corpuscular volume is decreased. The haemoglobin rate in a blood sample and haematocrit are thus reduced.

7

1.2. BLOOD Normal population Anemic population

Cell volume Normal (=) Cell count Normal (=)

Cell volume Decreased (↓) Cell count Normal (=)

Figure 1.3: Schematic representation of a red blood cell population under mycrocytic anemia. • Normocytic anemia: the red blood cell count is decreased but the mean corpuscular volume is normal. The haematocrit is decreased and the haemoglobin rate is decreased. Normal population Anemic population

Cell volume Normal (=) Cell count Normal (=)

Cell volume Normal (=) Cell count Decreased (↓)

Figure 1.4: Schematic representation of a red blood cell population under normocytic anemia. • Macrocytic anemia: the red blood cell count is low and the mean corpuscular volume is increased; it is often associated with a reduced haemoglobin content per cell. The haematocrit may be increased but the haemoglobin rate is decreased.

8

CHAPTER 1. INTRODUCTION Normal population Anemic population

Cell volume Normal (=) Cell count Normal (=)

Cell volume Increased (↑) Cell count Decreased (↓) + reduced Hb content

Figure 1.5: Schematic representation of a red blood cell population under macrocytic anemia.

All these conditions lead to a diminution of the haemoglobin level in blood, provoking various undesirable symptoms for the patient and being indicators of more serious diseases that would require further examination. These three types of anemia can be detected by the use of blood analyzers since they affect the red blood cell count and mean corpuscular volume values.

1.3 1.3.1

Blood analyzers Usage

An important objective of the diagnosis is to characterize the hematologic state of a patient, namely through the complete blood count (CBC). This is done thanks to blood analyzers. Blood analyzers have been used in most diagnoses for more than 40 years and their purpose is to count and size living cells and microparticles [1],[176],[199],[92], allowing the physician to obtain the CBC of a patient. A traditional complete blood count provides a detailed summary of the current state of a patient’s hematological system, allowing the determination of the following quantities: • Red blood cell count (4 to 6 million cells per mm3 ), • white blood cell count (4000 to 10000 cells per mm3 ), • platelets count (150000 to 400000 cells per mm3 ), • the total oxygen carrying capacity through the haemoglobin (120 to 175 µg/mm3 ) and haematocrit levels (31% to 53%), • the red blood cells mean corpuscular volume (MCV) (80 to 100 fL),

9

1.3. BLOOD ANALYZERS

• the mean corpuscular haemoglobin (MCH), which is the average amount of haemoglobin per red blood cell (27 to 31 picograms of haemoglobin per cell), • the mean corpuscular haemoglobin concentration (MCHC) (320 to 360 µg/mm3 ), • the red blood cell volume distribution width (RDW), which corresponds to the variation in the cells volume of the red blood cell population (11.5% to 14.5%): RDW =



Standard deviation of MCV mean MCV



× 100

(1.1)

• Mean Platelets Volume (MPV) (7.5 to 11.5 femtoliters (fL)). A variation in the number of red blood cells can be caused by various parameters, mostly pathologies. Namely, leukemias and haemopathies tend to modify the red blood cell count. The white blood cells count can be seen to vary in the event of an infection, when the organism is defending against a disease or when autoimmune disorders occur and destroy white blood cells [? ]. As for the other blood cells, the platelets population varies when their production in the bone marrow is abnormal, which can be caused by cancers or other forms of bone marrow dysfunctions. The mean corpuscular volume (MCV) and haematocrit levels which are closely related are subjected to variations in the case of various diseases that cause a decrease in the volume of the red cell population and/or a decrease of the number of red blood cells. The mean corpuscular haemoglobin (MCH) and mean corpuscular haemoglobin concentration (MCHC) values can vary because of mutations in the haemoglobin protein or diseases that cause a modification in the shape, size or integrity of red blood cells. The value of the red blood cell distribution width (RDW) expresses the variation of the corpuscular volume in a red blood cell population. High and low RDW values can be measured in various pathologies, such as anemias [? ]. When blood cells age in the circulatory system, many of their characteristics change during their life span: their volume decreases of around 20%, the haemoglobin mass decreases of 15%, the surface/volume ratio increases and the surface area decreases. Young red blood cells are known to be larger than mature cells, and the change in their volume is seen to be fast, occurring within the first week of circulation. Later volume changes are slower but always continues to decrease over a red blood cell lifetime. Red blood cells of higher (lower) volume can thus be over-represented in the volume distribution if most of the population is young (old), which can be caused by various pathologies, either way [91]. Thus, the quantities presented previously can be used in the diagnosis of several blood related diseases. In addition, the measurement of these quantities allows the detection of infections (increased white blood cell count) , drug intoxication (decreased platelets count), iron deficiency (decreased red blood cell count and decreased MCV) and other pathologies that are not directly blood related disease. A blood analyzer is able to measure some of these quantities thanks to the dielectric nature of red blood cells and

10

CHAPTER 1. INTRODUCTION

platelets, namely the MCV, RBC count, RDW, white blood cell count, platelet count and MPV. The principle of this kind of measurement is detailed in the next section. The hemoglobin content and concentration is measured through a destructive process that is not covered in this thesis, where the red blood cells are chemically forced to release their hemoglobin content, the latter being then measured with absorption spectrometry.

1.3.2

Coulter effect

A blood analyzer provides part of the CBC by means of the Coulter effect: a blood cell, considered insulating, is driven through a narrow channel filled with conductive fluid, where a constant electrical field is imposed. The perturbation created by the cell passing through this orifice is then measured and its amplitude is related to the cell volume (see Fig. 1.6) [41, 42].

a)

b)

c)

e2

C1 C2

e1

C2

Signal

C1

time

Figure 1.6: Schematic representation of the electrical measurement process inside a blood analyzer. (a) A concentrated electrical field is imposed in the orifice using electrodes e1 and e2 . b) Two red blood cells of different volumes pass through the orifice. c) Two pulses of different amplitudes are measured.

Repeating this operation for a larger sample (namely, a patient’s blood sample) gives access to the red blood cells count, the mean corpuscular volume (MCV) and red cell distribution width (RDW) defining the hematologic state of the patient. The Coulter effect is the basis of every commercialized blood analyzer to date, even though its first practical use has been conceptualized more than 60 years ago. Over the years, more and more sophisticated blood analyzers have been developed, allowing the probing of a large spectrum of quantities thanks to flow cytometry. Notably cell pigmentation, protein expression or DNA content were added as measurable quantities using optical instruments in combination with the Coulter effect. Nevertheless, a deep understanding of the dynamics of blood cells in blood analyzers and its influence on the CBC quantities is yet to be achieved.

11

1.3. BLOOD ANALYZERS

1.3.3

Hydrofocalization

The Coulter effect accurately measures the volume of red blood cells but suffers of inaccuracies depending on the cell orientation and position. These errors usually generate over-estimates of the higher volume red blood cell population. In order to reduce errors generated by the various trajectories of red blood cells within the orifice, the hydrofocalization technique is often used. Figure 1.7 shows a schematic representation of the

e2 Orifice

e1 Blood sample Sheathing Sheathing Figure 1.7: Schematic representation of a blood analyzer orifice neighborhood, with e1 and e2 the cathode and anode, respectively. hydrofocalization technique. Considering a symmetry of revolution around the blood analyzer main axis, a sheath flow is introduced all around the blood sample injector. This sheath flow has a greater velocity than the sample flow and restrains the blood sample to a defined region which depends on the respective flow rates of the blood sample and sheath flows. An example of 2-dimension hydrofocusing is shown in Fig. 1.8 with experimental results from Lee et al. [113]. It is a different setup, but it is easy to see that red blood cells will most likely stay in the darker region and be focused in the center of the orifice. This effect ensures that almost all red blood cells arrive in the measuring orifice centered and thus reduces the measurement errors. The profile of the velocity field in hydrofocalized systems suggests that red blood cells experiment an important velocity gradient in the direction of the flow. This velocity gradient, in addition to the intended centering of the cells on the micro-orifice axis, should modify the orientation and shape of the cells during their passage into the microorifice. This reorientation effect is well-known for solid particles and elongated colloids and will be studied in the case of red blood cells in Chapter 3 and 4.

12

CHAPTER 1. INTRODUCTION

Qs /Qi = 0.1

Qs /Qi = 0.5

Qs /Qi = 1.0

Qs /Qi = 2.0

Qs /Qi = 4.0

Qs /Qi = 8.0

Figure 1.8: Experimental images of the symmetric hydrodynamic focusing effect in microchannels with an aspect ratio of 1.78 and different flow rate ratios, where Qs is the sheathing flowrate and Qi the sample flowrate, from Lee et al. [113].

1.4

Numerical simulation of red blood cells and blood analyzers

As stated previously, the goal of understanding the dynamics of red blood cells in blood analyzers is tackled using numerical simulations. This section aims at presenting an overlook of the numerical simulations of particles in blood analyzers presented in the literature; Then, summarizing the numerical simulation of red blood cells and its difficulties, as well as replacing the choices presented throughout this thesis in the context of the current scientific advancement in the field. An extensive review was recently published by Freund [77] where he proposed a summary of the state-of-the-art of numerical simulations.

1.4.1

Numerical simulation of the red blood cell dynamics

The red blood cell membrane: a complex structure Red blood cells have been presented in detail above in this introduction. As a reminder, they are composed by a cytoplasm, which is an incompressible Newtonian fluid rich in hemoglobin. It is enclosed by a visco-elastic membrane, consisting of a lipid bilayer and a protein cytoskeleton. This membrane is almost incompressible as it strongly resists changes of its area. It also resists to shear and bending. These characteristics provide the red blood cells with their deformability. As long as surface and volume are conserved, they deform easily. The movement of red blood cells is usually either studied

1.4. NUMERICAL SIMULATION OF RED BLOOD CELLS AND BLOOD ANALYZERS

13

directly or by observing model particles, such as capsules or vesicles. Contrary to red blood cells and capsules, vesicles have a membrane that does not show resistance to shear, which results in vesicles experiencing slightly different movements [5]. Experiments are obviously an indispensable source of information and data. However, they have some limitations: the small length scales involved make the measurement of quantitative information difficult and the density of the red blood cells prevent a clear observation of the flow when the haematocrit is high. As a consequence, numerous analytical and numerical descriptions have been proposed over the last 30 years. Fluid-structure interaction When flowing, deformations of red blood cells can be impressive, especially in small tubes or channels, as shown by Fig. 1.9, for example.

Figure 1.9: From Abkarian et al. (2008) [2]. Time-lapse sequence of the deformation of a healthy RBC in a 5 µm channel.

This results from the fluid-structure interaction between the external fluid (plasma in blood), the internal fluid (cytoplasm) and the red blood cell membrane. The fluidstructure interaction is actually the reason why computing flows of red blood cells is a numerical challenge: • Red blood cells deform easily, and their shape may change continuously. As a consequence, they cannot be defined only with the position of the center of gravity and the orientation, as for non-deformable objects. Moreover, shape changes are large, so the small-deformation hypothesis cannot be used; • at characteristic flow scales relevant to micro-circulation, red blood cells cannot be considered as particles flowing passively with the fluid: flow deforms the red blood cells, which in turn modify the flow; • the mechanics of the red blood cell membrane is complex. ‘Macroscopic’ properties can be defined, but the membrane is inhomogeneous. Membrane modeling is thus a difficult task, where several levels of modeling, of increasing complexity, can be used; • in blood, hematocrit is of order of 45%. This high density population of red blood cells implies constant interactions between each other. At high shear rates, this

14

CHAPTER 1. INTRODUCTION interaction is limited to a dynamic interaction, but when the flow time scale is larger, red blood cells have time to aggregate, resulting in an increase of the blood viscosity [31, 80, 81, 159]; • when undergoing large deformations (often in non-physiological contexts, where shear rates can be higher than in human circulation), hemolysis can occur, where the red blood cell membrane breaks up or leaks part of the intra-cellular haemoglobin.

For all these reasons, computing the deformation and the dynamics of red blood cells in complex flows is a challenge. Such a challenge is tackled more and more efficiently, with the development of High Performance Computing over the last decades. The evolution of the analytical and numerical techniques developed over the years to predict the behavior of deformable particles in flows is reviewed in the next section.

1.4.2

Literature review: a historical perspective

Analytical prediction of the behavior of deformable particles in flows were mainly developed over the last 30 years, notably thanks to the major contribution of Barth`es-Biesel [17, 20, 21]. Barth`es-Biesel and co-workers introduced the ”capsule” model which refers to a deformable particle composed by an elastic membrane enclosing a drop of fluid. It was interestingly introduced as a model to study red blood cells in flows. Initially, Barth`es-Biesel and co-workers focused on the flow of capsules of spherical reference shape undergoing small deformations, considering an infinitely thin membrane. This can be considered as a good approximation for red blood cells, since their membrane is only a few nanometers thick. Under this hypothesis and assuming Stokes flow, they were able to determine analytically the time evolution of micro-capsules in linear flows, for different membrane mechanical laws. To tackle large deformations of capsules, the Boundary Element Method (BEM) was introduced by the group of Barth`es-Biesel [116, 161] and by Pozrikidis [152, 153, 154, 162]. The local deformation is computed over a discretized membrane and related to the local stress through the membrane mechanical law. The BEM allows the computation of the fluid flow in an infinite domain by integrating the membrane forces over the membrane surface, while considering a bounded domain increases the problem complexity and implies loss of accuracy [93, 94]. It takes advantage of the linearity of the Stokes flow describing the fluid motion. In parallel, vesicles started to receive attention as they can be considered as an alternate model for red blood cells. Vesicles refer to deformable particles consisting of a drop of fluid enclosed by a fluidic membrane: generally, the vesicle membrane is a bilayer of amphiphilic lipid molecules [167]. The vesicle membrane is different from the membrane of elastic capsules as defined by Barth`es-Biesel and is often considered as a two-dimensional perfect fluid layer: • it is mainly incompressible and the distance between lipid molecules of a layer does not change;

1.4. NUMERICAL SIMULATION OF RED BLOOD CELLS AND BLOOD ANALYZERS

15

• it does not resist shear; • however, it is not an interface: owing to its bilayer nature, it shows bending resistance, responsible for the large variety of shapes observable for vesicles [166, 167, 169, 170]. The pioneer work of Seifert is fundamental in the field, where he either used the small deformation theory and Lamb’s solution [168] (as Barth`es-Biesel earlier for capsules) or the BEM for triangulated surfaces [108]. One of his major contribution is also his work on vesicle morphology [166, 167, 169, 170]. In the 2000s, various approaches were introduced to compute the flow around and inside deformable fluidic particles undergoing large deformations. In the vast majority of the cases, they deal with three-dimensional deformable particles consisting of an infinitely thin membrane enclosing a simple Newtonian fluid. Both the internal and the external fluids are modeled and the particle can evolve in a more or less complex environment. From a global point of view, progress over the last 15 years in the simulation of deformable particles is impressive. Simulations account for: • different fluid properties inside and outside the membrane; • various mechanical behavior for the membrane: resistance to area changes, shear, bending. The membrane can be visco-elastic; • aggregation models have been proposed to model the clustering of red blood cells in low-shear flows and its consequences on blood viscosity; • numerical simulations of dense suspensions are being published in the 2010s. Simulations with several thousands of cells are now conceivable (though computationally expensive). However, simulations of red blood cells are far from being perfect: membrane mechanics can still be improved, validation of the results is still problematic and sometimes neglected to show impressive but non-validated calculations. In addition, red blood cells are often considered to flow in simple and low-Reynolds number flows. Such configurations, which correspond to physiological conditions, are not relevant to some extreme non-physiological conditions, as encountered in some medical devices. Haemolysis is also a challenge for future research on simulations of red blood cells in flows.

1.4.3

Description of existing numerical approaches for deformable particles

From the methodological point of view, the approaches encountered in the literature mainly differ from each other on the numerical treatment of the membrane. The physics of the membrane is not a criterion to differentiate the numerical methods as most papers

16

CHAPTER 1. INTRODUCTION

present results using several membrane constitutive laws. Depending on the research groups, fluid flow and membrane mechanics are treated using numerical different methods: first, the numerical methods are separated in three groups, depending on the way the membrane is treated. Then, the numerical method used to solve the fluid flow will separate the existing approaches in three other groups. Membrane approaches The continuous membrane approach This approach is the one used by the groups of Barth`es-Biesel and Pozrikidis [17, 20, 21, 108, 116, 152, 153, 154, 161, 162, 168]. The membrane is considered as a two dimensional surface of known mechanical properties: incompressible, elastic, visco-elastic, hyperelastic, resisting to bending, etc. The membrane is discretized and the local membrane forces related to its deformation state are computed. Discretization can be done through triangulation and spectral representation of the membrane [210, 219]. Forces are then calculated using finite elements, thin shell theory or others numerical methods. The majority of the numerical methods are based on the continuous membrane approach[13, 27, 51, 52, 62, 76, 93, 107, 109, 138, 145, 163, 189, 210]. The discrete membrane approach Alternatively, the membrane can be considered as a set of particles related by a spring network [64, 65, 66, 68, 103, 144, 151]. In this approach, there is no notion of discretization of a continuous membrane. Only elementary effects are defined by interaction of neighboring particles constitutive of the membrane. The physical parameters involved in the interaction between particles can be related to macroscopic properties of the cell. Gompper and co-workers use a particular membrane model, mixing the discrete membrane and the continuous membrane approach to model red blood cells as ‘elastic vesicles’ [129, 130, 140, 141, 142, 143]. A discrete membrane approach is also used by Ismail and Faure et al. [63, 97]. The implicit membrane approach In this approach, the membrane is not represented directly by a triangulated surface or a set of markers but an implicit definition of the membrane is used. A function Φ is defined over the fluid domain and an iso-value of this function marks the location of the membrane (typically, Φ = 0). The field of Φ can be transported using an Eulerian equation of transport, which is very convenient for parallel computing: only one grid is used for the computations, avoiding issues related to the exchange of information between the fluid and the membrane(s). Derivatives of Φ are conveniently used to reconstruct membrane information, as the curvature for example. A level-set approach has been developed by Cottet and co-workers [39, 40, 126, 127, 133] and more recently by Doyeux et al. [54]. A phase-field approach is also introduced by Misbah and co-workers [23, 25, 83].

1.4. NUMERICAL SIMULATION OF RED BLOOD CELLS AND BLOOD ANALYZERS

17

Fluid approaches BEM approach It has been seen that if the flow is described by the Stokes equation, the flow solution can be given using the Boundary Element Method. This is the approach notably used by Barth`es-Biesel and co-workers [58, 75, 93, 94, 110, 111, 112, 116, 161, 200, 201], Pozrikidis [155, 157, 158],Zhao and Freund [78, 139, 210, 211, 212, 213, 214, 216], Boedec and co-workers [27, 188] and others [51, 52, 145, 189]. Continuous approach More classical computational fluid dynamics (CFD) methods can also be used: finite differences are used by Popel [13, 60, 208, 209] and co-workers, Bagchi and co-workers [12, 14, 15, 16, 38, 48, 50, 207, 217] and Luo [124]. Finite differences are also used in the level-set formulation by Cottet and co-workers [39, 40, 126, 127, 133] or Kl¨opell [107]. Misbah et al. use periodical domain and Fourier transforms to discretize the space operators [23, 25, 83]. Finite element methods have also been used by several authors [54, 97, 121]. LBM approach The most famous one is the lattice-Boltzmann method (LBM) [6, 30]. Fluid flow is represented through the resolution of the motion of particles streaming and colliding. Under certain assumptions (as low Mach number), lattice Bolzmann equations are shown to tend towards the Navier-Stokes equations. LBM is used by Aidun and co-workers [6, 35, 125, 163, 204], Sui et al. [179, 180, 181] and Kr¨ uger et al. [109]. Particle approach Recently, particle methods have developed tremendously. Dissipative Particle Dynamics (DPD) is used by Karniadakis et al. [64, 65, 66, 68, 144, 151]. Gompper and co-workers use the multiparticle collision dynamics (MPC) framework [129, 130, 140, 141, 142, 143].

1.4.4

Simulation of red blood cells: state of the art

Numerical simulation is expected to provide information otherwise impossible to get analytically or by experimental means. In spite of the Stokesian character of the flow around small cells, numerical simulations are challenging due to the continuous deformation of the cells. Pozrikidis is one of the pioneers in the field of numerical simulations for red blood cells, using Stokes flow approximation and thin shell theory, his group [154, 155, 156, 157, 162] obtained a boundary integral formulation predicting the flow inside and outside a deforming membrane. The deformation of a unique cell in shear flow has been studied [157], as well as the axisymmetric movement of files of red blood cells in a cylindrical tube of small diameter [158]. This method suffered from numerical instabilities in some cases, which has led Kessler et al. [105] and Zhao et al. [210] to propose further improvements. Results shown thanks to the spectral

18

CHAPTER 1. INTRODUCTION

boundary integral method by Zhao et al. [139, 210, 215] are impressive (see Fig 1.10): this method is able to account for numerous cells, showing complex shape deformation without suffering from the numerical problems encountered in classical boundary integral methods. The group of Barth`es-Biesel has also used the boundary element method (BEM) to study the movement, buckling and interaction of capsules in shear flows [75, 110, 111, 112, 200, 201]. This method, despite its efficiency and precision, has some limitations. For example, in a situation where the zero Reynolds number assumption cannot be made (a situation encountered for blood flows in large veins and arteries, or even more in artificial devices such as blood analyzers), the resulting non-linearities of the fluid equations clearly rule out the use of the Green’s functions techniques on which the BEM method is based. Nevertheless, BEM results can be considered as reference numerical data. Few numerical errors are made and their results can be considered as ”exact” to validate other simulation techniques. This is also true for the early work from Barth`es-Biesel and co-workers, limited to small deformations [17, 20, 21]. Note that the BEM is still being developed by several groups. High-quality numerical results for vesicles were recently obtained, notably for sedimentation of vesicles, using a BEM framework [26, 27, 193, 194, 195]. In addition, sophisticated theory for small deformations around a spherical shape continues to be used by different groups to predict capsule or vesicle small deformations in flows [26, 135, 198].

Figure 1.10: From Zhao et al. (2010) [210]. Periodic three-dimensional cellular flow in a cylindrical tube of diameter 16.9 µm. Simulation with BEM. Flow is from left to right in the side view and toward the viewer in the end view.

Over the last years, several authors have focused on the use of the lattice-Boltzmann method (LBM). Sui and co-workers [179, 180, 181] have notably introduced the concept of immersed boundaries in the lattice-Boltzmann framework. They coupled this formulation with a finite element method used to obtain the force applied on the discretized surface of the cell. Aidun et al. [6, 7, 125, 204] also combine LBM and finite element analysis for the membrane, obtaining simulation results for dense suspensions (up to 50% in volume). All these simulations were able to show good qualitative behavior: Sui et al. [181] obtained the different types of motion of isolated capsule in shear flow (tumbling, tank treading and swinging) and the dependence on the shear rate applied.

1.4. NUMERICAL SIMULATION OF RED BLOOD CELLS AND BLOOD ANALYZERS

19

However, intermittent behavior observed in experiments [3] was not reproduced until recently by Cordasco et al. [37]. MacMeccan et al. [125] obtained the well-known shearthinning behavior of blood in their simulation of a suspension of capsules, the viscosity of the suspension decreasing when shear stress is increased. However, they obtained a suspension viscosity significantly under-estimating blood viscosity. More classical CFD was also used to perform simulations of red blood cells, but the deformation of the cells remains an important challenge. Lagrangian-Eulerian approaches are expected to be very precise, but reconstruction of the whole volumetric mesh at each time step is too expensive in terms of computational cost. Bi-dimensional and tri-dimensional simulations have been reported [107, 121], but they will be impossible to extend to simulations of 3D dense suspensions. Another possibility is to use full Eulerian approaches: Cottet et al. [39, 40], Maitre [126] and Maitre et al. [127] proposed one using a level set approach where the Navier-Stokes equations are modified by adding membrane forces related to membrane area changes and curvature. Though localized, these forces are smoothed over a small number of cells to enable correct resolution. They were able to obtain the biconcave discocyte shape at rest [127], as well as the tank treading and tumbling movements [126] in simple shear flow. A full Eulerian approach has also been developed by Misbah and co-workers, and applied specifically to vesicles [23, 25, 83]. Misbah also used an analytical approach to study the vacillating-breathing vesicle regime [135]. Immersed boundary methods were used by Bagchi et al. [13], Dooley and Quinlan [53] and Zhang et al. [208, 209] in 2D, using the front tracking method in a finite-difference framework. More recently, an important effort has been made in Bagchi’s group to perform 3D simulations. Doddi and Bagchi extended the method to 3D simulations [12, 14, 15, 16, 48, 50, 207, 217]. Li and Sarkar also presented front-tracking simulations [115]. The years 2010s have also seen the development of a new framework to calculate flows of red blood cells, coupling particle methods (Dissipative Particle Dynamics [64, 65, 66, 68, 144, 151], Multi-particle Collision dynamics [129, 130, 140, 141, 142, 143], Moving particle semi-implicit method [103]) and particle methods to describe the membrane behavior. The membrane is described by particles (as markers) which interact with each other locally. Particles are organized in a triangulated network. Simple interaction between particles is imposed: particles are assumed to be linked by springs to model the membrane elasticity, modified by a viscous term to recover the visco-elastic behavior of the membrane. Springs contribution also accounts for a term involving the temperature, to model thermal fluctuations of the membrane [65]. The main advantage of such a method is that macroscopic behavior of the membrane emerges from individual local interactions between particles modeling the membrane. Constants involved in the particle interaction models can be related analytically to macroscopic constants measured experimentally, for instance. To date, such models for deformable membranes are probably the most accurate to represent red blood cell membranes.

20

CHAPTER 1. INTRODUCTION

More recently, the effect of the reference shape of red blood cells has been investigated. The spheroidal stress-free configuration has been used instead of the usual biconcave red cell shape. It notably showed the importance of the reference shape in retrieving the transitions between the various regimes of the cell in shear flows [38, 107, 147, 189].

1.4.5

Configurations of interest in numerical studies of flows of deformable particles

Blood analyzers are complex devices and the study of the dynamics of cells within them is an active field of research [182]. The size of these devices, allowing no direct optical access, as well as the high velocity flows of several meters per second in the system, are important obstacles for experimentalists. An elaborate numerical framework and powerful solver are required to simulate the dynamics of deformable cells with sufficient accuracy and handle complex geometries encountered in blood analyzers. In the literature, the particles considered in numerical simulations usually are bigger and rigid. In addition, most of ex-vivo studies are performed in simplified geometries with low Reynolds number flows, which is different from the velocity fields commonly encountered in blood analyzers [24, 96, 98]. Most of the numerical studies focus on the characterization of the electrical pulse more than on the dynamics of the cells themselves [160, 178]. Many questions remain unanswered, notably on the cell dynamics in such systems and its influence on the measured electrical pulse. It is also of interest to understand the effect of the electrical field on the cell trajectory and shape with effects such as dielectrophoresis, electrorotation, etc. [47, 149, 197]. Finally, although the consequence of red blood cell deformation [11, 88] and integrity [10, 44] on its impedance has been investigated, its effect on electrical pulse measurements is still unknown. Although some of the publications discuss more complex effects such as electrokinetic forces, electrical double layer [86] or hydrofocusing [113], there is a missing link between the cells dynamics and the electrical pulse and there is no clear evidence on the nature of the dynamics of deformable particles in the high Reynolds number flows of blood analyzers. Alternatively, most numerical studies about red blood cells are focused on isolated red blood cells in simple configurations such as the simple shear flow [37, 189], red blood cells in suspension [68, 85], red blood cells in channels [67, 69, 142, 210, 212, 215] or in slits [77]. The study of the evolution of red blood cells in elongational flows is still an unexplored test case, even though it has been studied with vesicles [138, 214] or capsules [46, 51, 111]. This thesis is focused on red blood cells dynamics in such devices, as it is seen that studies about blood analyzers only consider rigid particles or no particle at all and that studies about red blood cells only consider simplified configurations or unrelated geometries.

1.5. STRATEGY

1.5

21

Strategy

The problem of the dynamics of red blood cells in a blood analyzer and their influence on the electrical measurement is tackled using a numerical method presented in chapter 2. This method is then tested against several test cases in order to assess the quality of its implementation. The dynamics of red blood cells are investigated by means of analytical developments in chapter 3: using the solid particle models of Jeffery [99], Keller & Skalak [104] and Abkarian, Faivre & Viallat [3] adapted to an elongational flow, the main physical effects behind the dynamics of red blood cells near the micro-orifice of the blood analyzer are found. In chapter 4, the dynamics of red blood cells in such a configuration are studied using 3D numerical simulations of deformable red blood cells in the blood analyzer. This study allows the understanding of the behavior of the orientation of the red blood cell when it travels inside the micro-orifice. This trend of orientation is modeled using the developments of chapter 3. Finally, the results of the simulation of the dynamics of red blood cells are used in chapter 5 to generate the electrical pulse generated by the red blood cell traveling through the micro-orifice. A statistical approach is used to understand the influence of the dynamics on this measurement by means of a Monte-Carlo algorithm, thus leading to the determination of the errors on the measurement due to the orientation, shape and trajectory of the cell.

Chapter

2

Numerical method Chapter contents 2.1

2.2 2.3 2.4

2.5 2.6 2.7 2.8

YALES2 & YALES2BIO . . . . . . . . . . . . . . . . . . . . 2.1.1 YALES2 . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 YALES2BIO . . . . . . . . . . . . . . . . . . . . . . . Modelling framework . . . . . . . . . . . . . . . . . . . . . . Fluid-structure coupling with the immersed boundary method Discretization of the coupled problem . . . . . . . . . . . . . 2.4.1 Hyperelastic forces . . . . . . . . . . . . . . . . . . . 2.4.2 Curvature forces . . . . . . . . . . . . . . . . . . . . . 2.4.3 Resolution . . . . . . . . . . . . . . . . . . . . . . . . Interpolation and regularization over unstructured meshes . Particle volume conservation . . . . . . . . . . . . . . . . . . YALES2BIO timestep description . . . . . . . . . . . . . . . . Validation test cases . . . . . . . . . . . . . . . . . . . . . . . 2.8.1 Pressurized capsules . . . . . . . . . . . . . . . . . . . 2.8.2 Equilibrium shapes . . . . . . . . . . . . . . . . . . . . 2.8.3 Optical tweezers . . . . . . . . . . . . . . . . . . . . . 2.8.4 Spherical capsules in linear shear show . . . . . . . . . 2.8.5 Red blood cells in simple shear flow . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

24 24 24 24 26 27 28 30 32 33 37 38 39 39 41 44 47 61

This chapter presents the numerical method used to perform the simulations of the dynamics of deformable particles, especially red blood cells. The method is based on the Front-Tracking Immersed Boundary Method (FT-IBM), originally developed by Peskin [148]. The whole numerical method is first presented here, as well as an additional correction algorithm used to ensure the conversation of the particle volume. Several test cases are presented in order to validate the implementation of the numerical method into the YALES2BIO solver. 23

24

2.1

CHAPTER 2. NUMERICAL METHOD

YALES2 & YALES2BIO

The numerical tool used in this thesis is YALES2BIO, an unstructured solver based on YALES2. These two numerical solvers are presented in this section.

2.1.1

YALES2

YALES2 is an unstructured solver developed at CORIA, Rouen, France (UMR CNRS 6614). It has been initially developed by Vincent Moureau in 2007. It is a finite volume unstructured solver designed for massively parallel computations on meshes that can reach a few billion elements [137]. It was originally created to study turbulent two-phase combustion, but now contains more solvers, allowing the simulation of sprays, acoustic and radiative phenomena.

2.1.2

YALES2BIO

The YALES2BIO solver is a numerical software for the simulation of blood flow [131]. It aims at helping the analysis of medical devices in contact with blood such as flow diverters, Ventricular Assist Devices, extra corporal circulation, artificial heart and valves, blood analyzers among others. YALES2BIO is based on YALES2 and inherits its massively parallel capabilities and high order finite volumes schemes for complex geometries. Both macroscale and microscale computations have been performed with YALES2BIO, such as the flow in a whole human left heart including cardiac turbulence in the work of Chanfa et al. [32, 33], inside an artificial heart, through an aortic valve or the response of a red blood cell stretched within an optical tweezers (presented as a validation test case later in this section) and the flow within an industrial blood analyzer [132], which is the subject of this thesis.

2.2

Modelling framework

The numerical method aims at modeling the transport of deformable particles, composed of an internal fluid enclosed by a flexible structure in a carrying fluid. Such a situation notably includes capsules and red blood cells. Both fluids inside and outside the cell are supposed to be incompressible and Newtonian. The fluid flow is thus governed by the continuity and the Navier-Stokes equations: ∇.~u = 0,

(2.1)

∂~u + ~u.∇~u) = −∇p + ∇.[µ(∇~u + (∇~u)T )], (2.2) ∂t where ~u is the fluid velocity, p the pressure and t denotes time. µ is the dynamic viscosity, which can be different for the internal and external fluids. Density variations are neglected so that the density ρ is a given constant. The flexible membrane is modeled as an infinitely thin massless structure, completely closed, that is transported by the ρ(

25

2.2. MODELLING FRAMEWORK

fluid. In dimension d, the membrane is a structure of dimension d − 1 (a surface in 3-D), denoted by S. Several mechanical effects can arise from the red blood cells membrane complexity: shear resistance, bending resistance and area dilation resistance. To represent the membrane mechanics, several approaches exist, which have been recently reviewed by Freund [77] and presented in Section 1.4. One can cite lattice spring models [65, 66, 68, 151, 171], multiscale mechanical approaches starting from the mechanical properties of proteins and using homogenization to obtain large-scale models [145, 146]. Several groups rely on continuous models, modeling the elastic part through an equivalent hyperelastic model and completing the approach by a model accounting for bending resistance [38, 206, 207, 210, 215]. The last framework is the one adopted here. The membrane is represented as a two-component mechanical structure, for which the in-plane resistances (to shear and area dilation) are modeled by an equivalent hyperelastic model, while the bending resistance of the lipid bilayer is treated separately. Note that the hyperelastic model accounts for the shear resistance of the cytoskeleton and the area dilation resistance of the lipid bilayer, while the bending model only accounts for the bending resistance of the lipid bilayer. Due to the value of their respective moduli [119], the bending resistance and the area dilation resistance of the cytoskeleton are neglected. The hyperelastic model used to represent the membrane of the cells was specifically developed for red blood cells by Skalak and co-workers [174]. This model is defined by the strain energy function WSK , calculated from the in-plane principal values of strain λ1 and λ2 , thanks to the membrane shear modulus Es and area dilation modulus Ea , expressed in N.m−1 . The Skalak law reads: WSK =

Es 2 Ea 2 2 [(λ1 + λ22 − 2)2 + 2(λ21 + λ22 − λ21 λ22 − 1)] + (λ λ − 1)2 . 4 4 1 2

(2.3)

It can also be written with the ratio of the area dilation modulus to the shear modulus, Ea C= : Es WSK =

Es 2 [(λ1 + λ22 − 2)2 + 2(λ21 + λ22 − λ21 λ22 − 1) + C(λ21 λ22 − 1)2 ]. 4

(2.4)

This law is considered as a good model of the mechanical properties of a red blood cell membrane. It takes into account the shear resistance through Es , representing the cytoskeleton contribution and the area dilation resistance through Ea , which represents the lipid bi-layer effects. In addition, its non-linear expression ensures a strain hardening effect [174]. This model takes into account both the membrane resistance to shear and area dilation. In addition, the bending resistance of the membrane is represented using a bending energy function Wb , proposed by Helfrich [89]: Wb =

Eb 2

Z

S

(2κ − c0 )2 dS.

(2.5)

26

CHAPTER 2. NUMERICAL METHOD

With Eb the bending modulus, κ the mean curvature respectively and c0 the spontaneous curvature. The models can be combined or used separately. A membrane only modeled as a hyperelastic shell, without bending resistance, is a good model for the thin shell of a capsule [75, 93, 94]. The capsule shell is also often modeled using another hyperelastic model, the Neo-Hookean model, which reads WN H =

Es 2 −2 (λ + λ22 + λ−2 1 λ2 − 3), 2 1

(2.6)

The Neo-Hookean model is used to model shear-softening material, while the Skalak model is used for strain-hardening membranes [19]. Lipid vesicles can also be modeled using the laws presented earlier in this section. One way to model vesicles is to use the Skalak model of Eq. (2.3) with Es = 0 and a large value of Ea [205]. In this case, the membrane has no shear resistance, but strongly resists area dilation. The Helfrich model accounts for the bending resistance.

2.3

Fluid-structure coupling with the immersed boundary method

The present framework relies on the immersed boundary method (IBM) and the fronttracking (FT) method. Let us first describe the coupling at the continuous level before discretization. The membrane movement is described in a Lagrangian manner; it is simply transported by the fluid. Adherence of the fluid over the membrane makes fluid velocity continuous at the membrane location and equal to the membrane velocity, thus, the membrane motion is described by: ~ dX(t) ~ (X(t), ~ =U t) = dt

Z

Ω

~ ~u(~x, t) δ(~x − X(t)) d~x.

(2.7)

~ denotes the membrane coordinates, U ~ the membrane velocity and ~u the fluid velocity. X The IBM is a one-fluid formalism, thus the fluid domain Ω is the union of the internal and external fluid domains. The velocity of the membrane is calculated from the fluid velocity field, using the Dirac function δ, ensuring the continuity of velocities at the interface. The effect of the particle deformation on the fluid motion is twofold: • Fluid properties are variable in space (e.g. viscosity may differ inside and outside the membrane, whose position varies over time) • As a reaction against its deformation, the membrane exerts a force f~ on the fluid. Thus, the fluid momentum conservation equation reads: ρ(

∂~u + ~u.∇~u) = −∇p + ∇.[µ(∇~u + (∇~u)T )] + f~. ∂t

(2.8)

2.4. DISCRETIZATION OF THE COUPLED PROBLEM

27

The force per unit volume f~ is obtained from the membrane forces per unit surface F~ , through: Z ~ ~ ~ f (~x, t) = F~ (X(t), t) δ(~x − X(t)) dS. (2.9) S

The addition of this source term to the fluid momentum equation ensures the right stress discontinuity across the membrane.

2.4

Discretization of the coupled problem

Figure 2.1: Schematic representation of the discretized problem in 2D. The membrane is discretized using a set of M Lagrangian markers (Fig. 2.1). The −−→ discrete unknowns are the coordinates of each marker m : Xm = (Xm , Ym , Zm )t . The t fluid domain Ω is discretized with N grid vertices located at − x→ n = (xn , yn , zn ) . The markers transport equation, Eq. (2.7) is advanced using an explicit Euler scheme. The discretization of the fluid equations over unstructured grids makes the transport equation for the membrane markers involve discrete Dirac functions wm , to interpolate the fluid velocity at the markers location from the fluid velocity at the t neighboring grid vertices, located at − x→ n = (xn , yn , zn ) : N − → −−→ X −→ − − →, t)w ( ||xn − Xm (t)|| )dV , Um (t) = u→ ( x n n m n h n=1

(2.10)

dVn is the volume of the element associated with the fluid grid node n. The definition of the discrete Dirac function wm is the object of section 2.5. The force per unit volume applied on the fluid by the membrane is computed from the discrete forces per unit surface on the membrane markers, through the process of force regularization, which also involves wm : −−→ M X − → −→ −−→ ||− x→ n − Xm (t)|| )dSm . fn (t) = Fm (Xm (t), t)wm ( h m=1

(2.11)

28

CHAPTER 2. NUMERICAL METHOD

−→ dSm is the surface element associated with marker m. Fm is the discrete force per unit − → surface applied by the membrane on the fluid grid around marker m. fn is the force per unit volume at grid vertex n. In two dimension, the membrane is discretized as a set of markers connected by edges as shown in Fig 2.1. In 3 dimension, the membrane surface is triangulated as presented on Fig. 2.2.

Figure 2.2: Triangulated red blood cell membrane.

2.4.1

Hyperelastic forces

To compute the elastic membrane forces, the method introduced by Charrier et al. [29] is used. It is a simple first-order finite element method. In order to present the whole algorithm, the method is presented for a flat membrane, deformed only in its plane. The presentation can thus be made in 2D. This method is easily extended to a 3-dimensional membrane by means of solid body transformations on the deformed elements, as will be explained at the end of the section. The method was already fully described by Charrier et al. [29], but it is detailed here for the sake of completeness and clarity.

(X i , Y i)

(x j, y j)

(X j , Y j)

(x i, y i)

y, Y

(x k, y k)

(X k , Y k ) vk uk

x, X Figure 2.3: Schematic of a membrane element deformation.

29

2.4. DISCRETIZATION OF THE COUPLED PROBLEM

In this section, we will use the following notations: a flat membrane is deformed in its plane. A membrane marker is described by its reference coordinates x and y, in the undeformed configuration and the current coordinates X and Y after deformation. The notation from Charrieret al. [29] is used, column vectors are denoted by { }. The displacements are noted u and v in the x and y directions, respectively (see Fig. 2.3). 1. The relation between the deformed and undeformed coordinates can be written as: X = x + u, (2.12) Y = y + v.

(2.13)

2. The membrane being modeled with a 2-dimensional hyperelastic law W , function of λ1 and λ2 the principal values of strain, one can find the force exerted by the elements with: a) The first order development of this function, which reads: δW

∂W ∂λ1 ∂W ∂λ2 + = {δu} ∂λ ∂u ∂λ2 ∂u  1    ∂W ∂W ∂λ ∂λ2 1 + {δv}T + , ∂λ1 ∂v ∂λ2 ∂v T











(2.14)

b) and the principle of virtual work: δWe = {δu}T {Fx } + {δv}T {Fy },

(2.15)

The element is assumed to undergo homogeneous deformation, thus: δWe = Ve δW,

(2.16)

with Ve the original volume of the element. Hence, one has: ∂W ∂λ1



∂λ1 ∂u



+ Ve

∂W ∂λ2



∂λ2 ∂u



,

(2.17)

∂W {Fy } = Ve ∂λ1



∂λ1 ∂v



∂W + Ve ∂λ2



∂λ2 ∂v



,

(2.18)

{Fx } = Ve

with Ve the original volume of the element. 3. In order to compute Eqs. (2.17) and (2.18), an expression for the principal values of stretch λ1 and λ2 is needed. Their value is computed for each triangle of vertices i, j, k. It is possible to find an explicit relation between the principal values of strain and the displacements using a linearity hypothesis on the displacements {u} and {v} in the triangle: {u} = Ni ui + Nj uj + Nk uk = {uT }{N },

(2.19)

30

CHAPTER 2. NUMERICAL METHOD {v} = Ni vi + Nj vj + Nk vk = {v T }{N },

(2.20)

with Ni , Nj and Nk being shape functions such that Ni = 1 on marker i and Ni = 0 on markers j and k, Nj and Nk being defined by cyclic changes of indices. To relate the coordinates of the undeformed (x, y) and deformed (X, Y ) elements markers, one writes: {dX} = [F ]{dx}, (2.21) with [F ] the transformation gradient which reads, [F ] =



∂X ∂x



= [I] +



∂u . ∂x 

(2.22)

4. The right Cauchy-Green deformation tensor, [G] = [F ]T [F ], can then be computed and principal values of strain λ1 and λ2 are related to the eigenvalues of [G], so that: (2.23) [G] = [F ]T [F ]. q 1 (2.24) λ21 = [G11 + G22 + (G11 − G22 )2 + 4G212 ], 2 q 1 λ22 = [G11 + G22 − (G11 − G22 )2 + 4G212 ]. (2.25) 2 G11 , G22 and G12 have explicit expressions as a function of the displacements and the shape functions used in the linearity hypothesis (see Charrier et al. [29] for details).

5. Once λ1 and λ2 are computed, it is possible to use Eqs. (2.17) and (2.18) and compute the contribution of one neighbouring element on the marker. In order to obtain the resulting force on a marker, the procedure has to be repeated for each neighbouring element and summed on the considered marker.

The method is easily extended to general 3-dimensional membranes as the calculations are made element by element. Thus, to compute the hyperelastic forces for a 3-dimensional membrane, solid body transformations are performed on the deformed elements, in order to put both the undeformed and deformed elements in the same plane. Since solid body transformations do not change the stretching of the elements, the problem remains the same.

2.4.2

Curvature forces

When the curvature energy of the membrane is described by the Helfrich equation of formula 2.5, the curvature force on the membrane, as found by Zhong-can et al. [218] can be written as: F~b = −Eb [(2κ − c0 )(2κ2 − 2κg + κc0 ) + 2∆s κ]~n.

(2.26)

2.4. DISCRETIZATION OF THE COUPLED PROBLEM

31

With Eb the bending modulus, κ and κg the mean and Gaussian curvatures, c0 the spontaneous curvature and ∆s the Laplace-Beltrami operator, which is a surface Laplacian operator. The method described in Farutin et al. [62] is used to compute κ, κg and ∆s κ from the coordinates of the membrane markers. A local coordinate system associated to each marker is introduced. It is composed of the approximate normal to the surface n~O for node a, the average of the neighboring faces normals, and two unit vectors ξ~O and η~O orthogonal to n~O and to each other. The origin of this local coordinate system is node ~ has the following O. In this coordinates system, a neighbor of node O, located at X, tangential coordinates: ~ − X~O ).ξ~O , sη = (X ~ − X~O ).η~O . sξ = (X

(2.27)

The surface is then locally characterized as a function of the tangential coordinates sξ and sη . In order to obtain the curvatures, a quadratic approximation of the global coordinates around node O is made (Xi , i = 1, 2, 3 denote the coordinates in the global system): Xi (sξ , sη ) = XiO + (∂ξ Xi )O sξ + (∂η Xi )O sη 1 [(∂ξξ Xi )O s2ξ + (∂ηη Xi )O s2η + 2(∂ξη Xi )O sξ sη ]. + 2

(2.28)

The five unknown coefficients of Eq. (2.28) are computed using least squares fitting of the approximation by comparison with the actual coordinates of the markers around marker O (see details at the end of the subsection). The mean and Gaussian curvatures at node O are then defined as: 1 κO = T r[cO (g O )−1 ], 2

O O −1 κO g = det[c (g ) ],

(2.29)

with, O = (∂ X )O (∂ X )O , gαβ α i β i

O O cO αβ = ni (∂αβ Xi ) ,

(2.30)

(α, β ∈ {ξ, η}, i = 1, 2, 3). O and cO are the metric and curvature tensors, respectively. Once κ is calculated for gαβ αβ each marker, the same procedure is applied to calculate ∆s κ. The mean curvature is approximated using a quadratic approximation:

κ(sξ , sη ) = κO + (∂ξ κ)O sξ + (∂η κ)O sη 1 + [(∂ξξ κ)O s2ξ + (∂ηη κ)O s2η + 2(∂ξη κ)O sξ sη ]. 2

(2.31)

The Laplace-Beltrami operator then writes: ∆s κ = p

q 1 −1 ∂β κ), ∂α ( |detg|gαβ |detg|

(2.32)

32

CHAPTER 2. NUMERICAL METHOD

which also reads: O −1 O O −1 O O ∆s κO = (∂αβ κ)O (g O )−1 αβ − [(g )αβ (∂αβ Xi ) ][(g )γδ (∂γ κ) (∂δ Xi ) ].

(2.33)

The fitting procedure consists in finding a paraboloid surface patch approximating the membrane locally, around each node. It is based on the procedure described in Garimella et al. [82]. A quadric patch f = f (sη , sξ ) defined as, 1 2 1 asη + bsη sξ + cs2ξ + dsη + esξ = f − f O , 2 2

(2.34)

is searched to fit the membrane around the node. In the case of Eq. (2.28) for instance, a, b, c, d and e correspond to (∂αα Xi )O , (∂αβ Xi )O , (∂ββ Xi )O , (∂α Xi )O , (∂β Xi )O , respectively. To find the values of the coefficients best fitting the variables of interest, Eq. (2.34) is written for a set of m markers neighboring O. The equations can be gathered in a unique system of the form S~a = f~

1

′2 2 sη1

  

s′η1 s′ξ1 .. .

.. . 1 ′2 ′ ′ s 2 ηm sηm sξn

1 ′2 2 sξ1

s′η1 .. .

.. . 1 ′2 ′ s 2 ξm sηm

    a  s′ξ1  f1 b    ..   c =  .      fm s′ξm d

e

− fO  .. . .  

−

(2.35)

fO

This system is over-constrained as soon as more than five neighbours of O are used. To find the values of the coefficients (gathered in vector ~a) that provide the best fit, a least-square fitting approach is used. The least-square problem can be written as a linear system to solve [59]: ~a is then obtained explicitly, provided that the inverse of (ST S) is calculated: ~a = (ST S)−1 ST f. For each node, this fitting procedure is applied four times for each marker: for each of the three global coordinates and for the mean curvature.

2.4.3

Resolution

The incompressible Navier-Stokes equations are solved using a projection method [34]. The fluid velocity is first advanced using a 4th-order central scheme in space, and a 4th-order Runge-Kutta like scheme in time [32]. In YALES2BIO, the forces exerted by the membrane on the fluid are included in the prediction step. Divergence-free velocity is obtained at the end of the time-step by solving a Poisson equation for pressure, to correct the predicted velocity. A Deflated Preconditioned Conjugate Gradient algorithm is used to solve this Poisson equation [128]. The correction step is not modified by the IBM. The fluid velocity at the end of each time step is interpolated to the markers, to convect the membrane.

2.5. INTERPOLATION AND REGULARIZATION OVER UNSTRUCTURED MESHES

2.5

33

Interpolation and regularization over unstructured meshes

In the original IBM formulation, as in a majority of studies using IBM, finite differences are used to solve the fluid equations, using a regular Cartesian grid of constant grid size h. In this case, the discrete Dirac function wm involved in the velocity interpolation (Eq. (2.10)) and force regularization (Eq. (2.11)) can be easily defined. wm is taken as a product of one-dimensional delta functions: → − − ||→ x n − X m || xn − Xm yn − Ym zn − Zm wm ( ) = D( )D( )D( ). (2.36) h h h h The cosine representation is often used in each direction:

D(r) =

    πr 1   if |r| < 2   4h 1 + cos 2     0

(2.37)

if |r| ≥ 2

The volume of the grid cells is also defined in terms of grid size: dVn = h3 . Note that more sophisticated delta functions were also proposed: see for example the discussions by Roma et al. [164] and Peskin [148] about the discrete Dirac functions and their properties. Eventually, note that over regular Cartesian grids, the discrete Dirac functions can be defined independently of the position of the markers. Over unstructured grids, Cartesian versions of wm cannot be used. This issue was dealt with in the context of finite element by Liu and coworkers [202] and their ideas were adapted later to the finite difference/finite volume context by Pinelli et al. [150]. Adaptation of the immersed boundary formalism to unstructured meshes relies on the Reproducing Kernel Particle Method (RKPM) [120]. The principle of the interpolation and regularization procedure with the Reproducing Kernel Particle Method of second order (RKPM2) is briefly explained in 1D before detailing the implementation in the 3D case. The effect of the RKPM and its order will be studied in section 2.8.4. ~ a force is Let us consider the regularization problem: at one marker located at X, known and has to be regularized over the fluid grid. For a fluid grid vertex of coordinate ~x, one has to determine the regularization weight, the value of the regularized discrete Delta function (also called window function here). Starting with an initial window function φX~ as the cosine one (Eq. (2.37)), one can introduce the following modified window function φX~ : K ~ X ~ ~ ~x − X ~x − X ~ ~x − X )k . ) = φX~ ( ) βk (X)( φX~ ( h h h k=0

(2.38)

~ would be the marker coordinate and ~x the coordinate of the grid vertex where the X weight φX~ is to be calculated. βk are the coefficients of the polynomial correction of

34

CHAPTER 2. NUMERICAL METHOD

the original window function. They are calculated by imposing conditions on the first moments of φX~ . The p-th moment of the function φX~ and φX~ is defined as: ~ ~ ~x − X ~x − X )p φX~ ( ) dx, h h Ω Z ~ ~ ~x − X ~x − X ~ mp (X) = ( )p φX~ ( ) dx. h h Ω

~ = mp (X)

Z

(

(2.39)

By plugging the definition of φ into Eq. (2.39), one easily obtains that:

~ = mp (X)

n X

~ m(k+p) (X), ~ βk (X)

(2.40)

k=0

~ the (k + p)th moment of φ at the location of the marker, X. ~ For a with m(k+p) (X) ~ moments calculated at X ~ are known: the first unit point force applied at coordinate X, moment is 1 and the following ones are all 0. These conditions, which are naturally true for the force to regularize, are imposed on the regularized force. Verification of the moments conditions ensures the quality of the regularization process. In 1D, imposing ~ = 1 and mp (X) ~ = 0 for 1 ≤ p ≤ P , one can determine the corresponding that m0 (X) βk for 0 ≤ k ≤ P , by inverting the (P + 1) × (P + 1) symmetric matrix [150]. When ~ =1 discrete Dirac functions are used for regularization of a marker force, having m0 (X) ~ = 0 guarantees that the regularized force has the same intensity and the and m1 (X) same application point as the marker force.

4h

marker m

node j y

x

dVj

Figure 2.4: Schematic representation of the procedure to compute the window function for regularization and interpolation over unstructured grids.

Extension to three dimensions is shown in Pinelli et al. [150]. In practice, the following algorithm is used to calculate the weights of the window function associated

2.5. INTERPOLATION AND REGULARIZATION OVER UNSTRUCTURED MESHES

35

with each marker. The a priori window function is: ~ m , ~xj ) w(X = + + +

~m ~xj − X xj − Xm yj − Ym zj − Zm )[(β0 )m + (β1 )m + (β2 )m + (β3 )m h h  h     h  xj − Xm yj − Ym zj − Zm yj − Ym (β4 )m + (β5 )m h h h h   2   xj − Xm xj − Xm zj − Zm + (β7 )m (β6 )m h h h   2 2 yj − Ym zj − Zm (β8 )m + (β9 )m ]. (2.41) h h

~ m, w(X

For one marker m, the algorithm is:

1. A local grid size h is defined. Its value is unique, typical of the fluid grid size at the particle location. Different values for each marker could be used but have not been tested.

2. For all the J fluid nodes located in the sphere of diameter 4×h centered on marker m, the intermediate window function weights are calculated (j = 1, ..., J): ~ ~ m , ~xj − Xm ) = 1 D w(X h 2d−1

! → − − || X m − → x j || h

(2.42)

3. Moments of the intermediate window function are calculated as:

~ m) = mp,q,r (X

     J  X xj − Xm p yj − Ym q zj − Zm r

j=1

h

h

h

~ m, w(X

~m ~xj − X )dVj , h (2.43)

where dVj is the volume of the element associated with vertex j.

4. Coefficients of the modified window function are calculated: for the modified window function to be a good regularization of a point force, the values of moments of order 0, 1 and 2 can be imposed: m0,0,0 = 1 while the others moments are set to 0. This necessitates, for each marker, the use of 10 coefficients (βk )m , (0 ≤ k ≤ 9 in 3D). Using Eq. (2.41), one can obtain the equivalent of Eq. (2.40) to relate the (βk )m

36

CHAPTER 2. NUMERICAL METHOD coefficients to the imposed moments: 











m0,0,0 (β0 )m (β0 )m       (β1 )m  m1,0,0  (β1 )m        (β )  m  (β )   0,1,0      2 m 2 m       m0,0,0 m1,0,0 · · · m0,0,2  m0,0,1  (β3 )m  (β3 )m         m  (β )   m1,0,0 m2,0,0 · · · m1,0,2   (β )  1,1,0       4 m 4 m  = .  ≡ M3D   ..  .. ..  . m0,1,1       (β ) (β ) . .   5 m .    .  5 m m  (β )  (β )   1,0,1   6 m 6 m m0,0,2 m1,0,2 · · · m0,0,4        m2,0,0  (β7 )m  (β7 )m        m (β )   (β )   0,2,0   8 m  8 m m0,0,2 (β9 )m (β9 )m

(2.44)

In order to calculate the (βk )m coefficients, the 3-D moments matrix M3D is ~ m is omitted) from the moments of the built (the dependence of moments on X intermediate window function calculated using Eq. (2.43). M3D is then inverted for each marker. The coefficients (βk )m are then calculated using: 







 

1 m0,0,0 (β0 )m       0 (β1 )m  m1,0,0         0 (β )  m    0,1,0   2 m       0 m0,0,1  (β3 )m           m (β )   4 m −1 0 −1  1,1,0   = M3D   .   = M3D  0 m0,1,1  (β5 )m        0  (β )  m    1,0,1   6 m       0 m2,0,0  (β7 )m         0 m (β )     0,2,0   8 m m0,0,2 0 (β9 )m

(2.45)

5. The weights of the modified window function can thus be calculated using Eq. (2.41). These weights are used to regularize the membrane forces on the fluid grids as well as for the interpolation of the fluid velocity to the membrane markers, using: f~n (~xn , t) =

M X

~ ~ m , ~xn − Xm ), ~ m (t), t)dSm w(X F~m (X h m=1

~ m (X ~ m , t) = U

N X

n=1

~ m, ~un (~xn , t)w(X

~m ~xn − X )dVn . h

(2.46)

(2.47)

Note that the window function w depends both on the fluid grid and on the markers. To avoid computing the weights twice per iteration (regularization and interpolation, their size being N × M ), (βk )m (data of size 10 × M ) are stored after the regularization step of Eq. (2.46) to be used for the interpolation step Eq. (2.47).

2.6. PARTICLE VOLUME CONSERVATION

37

Note that this algorithm has been modified by Pinelli et al. to avoid ill-conditioned M3D matrices. Such an issue has not been encountered in the present study, presumably because of the larger stencil used here. The discrete Dirac function is also used to update the viscosity field. As in the fronttracking method presented by Unverdi and Tryggvason [190], an indicator function In = I(~xn , t) is calculated to determine if ~xn is located inside (In = 0) or outside (In = 1) the particle. The dynamic viscosity thus reads: µ(x~n , t) = µint + (µext − µint )In ,

(2.48)

with µint the internal viscosity and µext the viscosity outside the particle. I(~xn , t) is determined as the solution of a Poisson equation: ~ ∆I = ∇.G,

~ x~n ) = with G(

M X

~m ~xn − X ~nm dSm w(X~m , ). h m=1

(2.49)

~ the right hand side of the Poisson equation, can be seen as the regularization field of G, the membrane unit normal vectors (~nm , m = 1, ..., M ). This technique ensures that the viscosity has the right value inside and outside the particle, with a smooth transition in the membrane region, where the forces are applied.

2.6

Particle volume conservation

The IBM does not intrinsically conserve the volume enclosed by the flexible membrane as the interpolation process does not conserve the divergence-free character of the carrying fluid flow. A Lagrange Multiplier method is used to find the smallest correction to the markers location ensuring the conservation of the particle volume. In three dimensions, the cell membrane is a mesh of F triangular elements. Each face f has three nodes i, j and k. The volume enclosed in the membrane can be calculated as: F − → − → − → − → − → − → − → − → − → 1X V (~x) = [xfi .(xfj × xfk ) − xfj .(xfk × xfi ) + xfk .(xfi × xfj )]. (2.50) 6 f =1 At the beginning of the computation, the initial volume of the particle V0 is computed. At the end of each time step, the current volume of the particle is also computed. Due to the numerical integration, they can be different. In order to obtain a particle of volume ~ which ensure that V (~x + δx) ~ = V0 . V0 , we search the smallest markers displacements δx ~ minimize Introducing a Lagrange multiplier Λ, the sought correction displacements δx the following cost function: ~ = JΛ (δx)

M X

m=1

~ − V0 ). ((δxm )2 + (δym )2 + (δzm )2 ) + Λ(V (~x + δx)

(2.51)

Partial derivatives of JΛ are calculated against δxm , δym and δzm . As these derivatives are zero at the minimum, one finds: δxm = −

Λ → Λ → Λ −→ [αm ]x ; δym = − [− αm ]y ; δzm = − [− αm ]z , 12 12 12

(2.52)

38

CHAPTER 2. NUMERICAL METHOD

− with → αi defined for marker i as: F X → − − − αi = (→ x fj × → x fk ).

(2.53)

f =1

Λ When JΛ is minimum, ∂J ∂Λ = 0, which means that the current volume is equal to the initial volume. Using Eq. (2.50) and Eqs. (2.52), the following third-order polynomial equation in Λ is obtained: AΛ3 + BΛ2 + CΛ + D = 0, with:

M −→ −→ −→ −→ −→ −→ −→ −→ −→ 1 X [αim .(αjm × αkm ) − αjm .(αkm × αim ) + αkm .(αim × αjm )], 3 12 m=1

(2.54)

M −→ −→ −→ −→ −→ −→ −→ −→ 1 X m m −→ m m m m m m [(xm i + xj + xk ).Vaa + (αi + αj + αk ).(Vax + Vxa )], 2 12 m=1

(2.55)

A=

B=

M −→ −→ −→ −→ −→ −→ −→ −→ −→ 1 X m m m m m [(αm + αjm + αkm ).Vaa + (xm C= i + xj + xk ).(Vxa + Vax )], 12 m=1 i

D = V (~x) − V0 . With,

(2.56) (2.57)

−→ −→ −→ −→ −→ m = (αkm − αjm ) × (αim − αjm ), Vaa −→ −→ −→ −→ −→ m m m m = (xm Vxx k − xj ) × (xi − xj ),

−→ −→ −→ −→ −→ m m Vax = (αkm − αjm ) × (xm i − xj ),

−→ −→ −→ −→ −→ m m m m Vxa = (xm k − xj ) × (αi − αj ).

This third-order polynomial equation is then solved using Cardano’s method. Once Λ is found, the membrane markers positions are updated to ensure volume conservation, as follows: − −→ Λ − α→ (2.58) x→ m. m = xm − 12

2.7

YALES2BIO timestep description

In this section, a regular timestep of YALES2BIO is detailed. Considering an initial configuration of a red blood cell inside a fluid domain, it goes: 1. The viscosity field is calculated by computing the indicator function (see section 2.5). 2. Hyperelastic and bending forces are computed on each marker of the membrane (see section 2.4). 3. These forces are regularized over the fluid grid (see section 2.5).

2.8. VALIDATION TEST CASES

39

4. The Navier-Stokes equations are explicitly advanced using the regularized forces as a source term using Chorin’s projection method (see section 2.4). 5. The new velocity field is interpolated on the membrane markers (see section 2.5). 6. The membrane markers are convected using an Euler time-advancement scheme (see section 2.3). 7. The markers position is corrected by the volume conservation algorithm (see section 2.6).

2.8

Validation test cases

The two dimension validation of the implementation of the presented numerical method has been presented in a previous publication, Mendez et al. [132]. It shows multiple test cases that prove the quality of the numerical results, the coupling and the algorithm in two dimension. In this section, the accuracy of the three dimension extension of the method and of the three dimension membrane mechanics are tested. First, a few static test cases are presented: pressurized capsules, equilibrium shapes of vesicles and red blood cells and a red blood cell undergoing the optical tweezers test case. Finally, an extensive test case of the dynamics is presented with spherical capsules in linear shear flow.

2.8.1

Pressurized capsules

The pressurization test case is often performed in the literature [201]. A spherical capsule is placed in a cubic fluid domain of side length 2L and is inflated from a radius a to radius ap = (1 + α)a by an internal pressure, the surrounding fluid is at rest (see Fig. 2.5). The capsule undergoes an isotropic tension associated to the stretch ratio λ = 1 + α. One can analytically relate the tension and stretch ratio [201]. In the case of the Neo-Hookean law, this relation reads: 1 (2.59) T = Es (1 − 6 ), λ and in the case of the Skalak law, T = Es (λ2 − 1 + Cλ2 (λ4 − 1)).

(2.60)

The isotropic tension of the membrane is directly related to the pressure difference across the membrane using Laplace’s law: 2T 2T p= = . (2.61) ap (1 + α)a Simulations are performed with a capsule which surface mesh contains 4580 elements and a cubic fluid grid of side 2L = 4a with solid walls and a uniform mesh size of ap . Without changing ap , different values of α are used according to table 2.1. dx = 12.5

40

CHAPTER 2. NUMERICAL METHOD

ap

2L

a

y x z Figure 2.5: Schematic representation of the pressurization test case.

Simulations parameters ap Es (N.m−1 ) C ap dx 12.5 1.0 1.0 0.5

L ap

2

α 0−0.5

Table 2.1: Summary of the simulations parameters for the pressurized capsule test case. ap , with p the pressure The comparison is done with the non-dimensional pressure E s difference between the inner and outer fluids. Using Eq. (2.59), Eq. (2.60) and Eq. (2.61), the analytical non-dimensional pressure difference is computed and compared with the numerical results on Fig. 2.6. A good agreement is found between the results computed with YALES2BIO and the analytical ones, the maximum error is found to be less than 0.5%.

2.8. VALIDATION TEST CASES

41

Figure 2.6: Non-dimensional pressure difference across the membrane for different values of the inflation factor α in the cases of the Skalak law (dashed) and the Neo-Hookean law (solid).

2.8.2

Equilibrium shapes

Equilibrium shapes of vesicles Equilibrium shapes of vesicles are computed. Vesicle dynamics are governed by bending forces as well as a surface incompressibility constraint, which implies that the vesicle with an initial ellipsoidal shape will deflate to a discocyte shape. Thus, the Skalak hyperelastic model is used with a shear modulus Es = 0 (providing no shear resistance), and high area dilatation modulus Ea (giving the surface incompressibility constraint). The computed equilibrium shapes are compared with semi-analytical results extracted in Boedec et al. [26]. The equilibrium shapes should depend only on the reduced volume v, defined as: V , (2.62) v= 4 3/2 3 π (A/4π) where V , A are the volume and area of the vesicle. Starting from prolate initial shapes of large radius a, the vesicles gradually deform to reach their equilibrium shapes. The time needed for the vesicle to reach the equilibrium shape depends on its size, the membrane moduli and the viscosity of the fluid. However, only the equilibrium shape, which depends on v, is compared with the reference data. The vesicles are meshed with 5028 and 7930 elements for v = 0.783 and v = 0.902, respectively. The fluid domain is a cube of size 2L = 2.4a, closed by solid calls, in a is used. Figures 2.7 and 2.8 show comparisons which a uniform grid of size dx = 12.5 of the computed equilibrium shapes with semi-analytical results. The shapes are found

42

CHAPTER 2. NUMERICAL METHOD

y/a

to be in excellent agreement.

x/a

y/a

Figure 2.7: Comparison of the final shape from YALES2BIO calculation with the semianalytical shape for v = 0.783.

x/a Figure 2.8: Comparison of the final shape from YALES2BIO calculation with the semianalytical shape for v = 0.902.

43

2.8. VALIDATION TEST CASES Equilibrium shapes of red blood cells

The reference shape of the cytoskeleton of red blood cells is not known and is still debated [38, 55, 73, 119, 147]. Although the regular discocyte shape of the red blood cell could be the obvious choice, one of the candidate for the resting shape of red blood cells is an ellipsoid with a reduced volume close to 1.0, implying pre-stress in the red blood cells (Fig. 2.9). This reference shape has notably allowed to reproduce the wellknown Stomatocyte-Discocyte-Echinocyte (SDE) sequence using shape optimization algorithms. The SDE sequence cannot be reproduced with the discocyte as a reference shape. In addition, using a spheroidal stress free shape has enabled recent simulations to obtain results in shear flow much closer to experimental data [38, 147, 189]. The test case consists in showing that the red blood cell shape can be obtained by deflating an ellipsoid and letting it evolve to a stable state [106, 118, 119]. In this test case, an initially ellipsoidal red blood cell is suspended into a resting fluid. Its reference shape is imposed to be an ellipsoid with a reduced volume of 0.95 and an area equal to the one of a red blood cell, ARBC = 133.4 µm2 . At the beginning of the calculation, the volume of the capsule is instantaneously imposed to be Vref = VRBC = 93.5 µm3 . This corresponds to a violent deflation of the capsule. The initial volume is thus smaller than the reference one. The capsule progressively relaxes to recover its reference volume, but at a fixed area.

Figure 2.9: Initial and final shape of the red blood cell.

It deforms into a discocyte in order to minimize the energy associated to the elastic and bending forces (Fig. 2.9). In this test case, the red blood cell is modelled using the Skalak law, with Es = 3.6 × 10−6 N.m−1 , C = 100 and Eb = 2 × 10−19 N.m. A parametric expression of the discocyte shape of red blood cells has been calculated from experimental measurement by Evans and Fung [61]. It reads:

Z = ±0.5R0

"

X2 + Y 2 1− R02

#1  2 X2 + Y 2 C0 + C1 + C2 2

R0

X2 + Y 2 R02

!2  ,

(2.63)

44

CHAPTER 2. NUMERICAL METHOD

with R0 = 3.91 µm, C0 = 0.207161, C1 = 2.002558 and C2 = −1.122762.

(2.64)

The profile of the red blood cell membrane obtained in YALES2BIO is compared to the parametric profile from Evans and Fung [61] on Fig. 2.10. In order to quantify the mesh refinement, a characteristic length is introduced: RRBC =



3VRBC 4π

1/3

,

(2.65)

With VRBC the volume of the red blood cell. The fluid mesh used is a Cartesian mesh with RRBC dx = 13, the membrane is meshed with 5120 surface elements.

Figure 2.10: Comparison of the parametric profile of a red blood cell and the numerical profile obtained. A good agreement is found between the two profiles. A shape very close to the parametric shape of Evans and Fung [61] can thus be recovered by minimizing the membrane energy (elastic and curvature energy) with an ellipsoidal reference shape for the cytoskeleton. As shown by Cordasco et al. [38] and Peng et al. [147], one can also impose a non-zero spontaneous curvature to obtain an even better match with the parametric shape.

2.8.3

Optical tweezers

The optical tweezers test case is an experimental setup where a red blood cell is stretched. Two silica microbeads, each 4.12 µm in diameter, are attached to the cell at diametrically opposite points.

45

2.8. VALIDATION TEST CASES

← ← ← ←

Laser beam

Figure 2.11: Schematic representation of the optical tweezers test case.

The left bead is anchored to the surface of the glass slide, and the right bead is trapped by a laser beam. As the trapped bead remains stationary, moving the slide and the attached left bead stretches the cell until the trapped bead escapes the laser, as the force exerted by the cell on the bead is higher than the one exerted by the laser. The cell is stretched at different laser power to record cell deformation over a range of forces. For each stretch test, the axial diameter (in the direction of stretch), and the transverse diameter (orthogonal to stretch direction) of the cell are measured when the trapped bead just escapes the trap. Such an experimental device is used to produce strain-stress curves for the whole cell to characterize membrane mechanics in healthy and pathological cases [43, 90, 114, 134, 183]. Computational domain For the cell geometry, Eqs. (2.63) and (2.64) are used. The fluid geometry is a parallelepipedic domain centered on the cell. Assuming that x is the direction of the stretching and y the axis of revolution of the red blood cell, the dimensions of the parallelepiped in x, y, and z directions are respectively Lx , Ly and Lz . The red blood cell surface is meshed with 3310 triangular faces and the fluid grid is discretized such as Lx = 6RRBC , Ly = 2RRBC and Lz = 4RRBC . Boundary conditions Specific conditions have to be defined on the cell membrane in order to stretch it. The assumption made is that the contact area between the beads and the cell can be considered as the intersection between the cell surface, and a plane perpendicular to the direction of stretch. The plane location is chosen such as the contact area length is equal to 2 µm in the z direction (shown on Fig. 2.12). To stretch the cell, the chosen method consists in applying an external force on the edge nodes of the contact area.

46

CHAPTER 2. NUMERICAL METHOD

Figure 2.12: Contact area between the beads and the red blood cell.

This method has been designed by J. Sig¨ uenza [172, 173] and has shown to eliminate the problem of applying a constant stress on the zone of the cell in contact with the bead. Such a technique generates unphysical pointwise oscillations, as previously shown [62, 107]. Numerical results The red blood cell is modeled with the Skalak law: Simulations parameters R0 Es (N.m−1 ) C dx −6 12.5 2.5−5.7×10 1.0

Eb 0.0

Law Skalak

Table 2.2: Summary of the simulations parameters for the optical tweezers test case. Transverse and axial diameters are measured and compared to experimental results from Mills et al. [134] over a defined range of stretching force values.

47

2.8. VALIDATION TEST CASES 20 Mills et al. (2004)

17,5

E =2.5E-6 2.5 µN.m−1 E_s s = E_s E =3.65E-6 3.65 µN.m−1 s = E_s E =5.7E-6AAAAAAA 5.7 µN.m−1 s =

Diameter YY(µm)

15 12,5 10 7,5 5 2,5 0

0

25

50

75

100

125

150

175

200

AppliedXXforce (pN) Figure 2.13: Axial (top) and transverse (bottom) diameters of the red blood cell for different applied forces in the optical tweezers test case and different shear modulii Es compared to the experimental results of Millset al. [134].

A good agreement is found between the experimental results from Mills et al. [134] and those produced by YALES2BIO.

2.8.4

Spherical capsules in linear shear show

In this test case, which setup is represented in figure 2.14, an initially spherical capsule of radius a is freely suspended in a linear shear flow, defined by ~u = (ky, 0, 0), with k the shear rate. The membrane is an infinitely thin hyperelastic structure.

48

CHAPTER 2. NUMERICAL METHOD

y z

x

Figure 2.14: Schematic representation of the test case setup.

In this prescribed flow, the capsule deforms in the shear plane (in this case, xOy) during an initial transient phase.

Figure 2.15: Successive states of a capsule during the transient phase.

t0

t1 > t0

t2 > t1

t3 > t2

t4 > t3

Figure 2.16: Schematic representation of the tank-treading motion. Successive positions of a specific marker (black square) on the capsule membrane along a run. After this transient phase (Fig. 2.15 and 2.17), it stabilizes at a steady orientation and shape and starts a tank-treading motion (Fig. 2.16 and 2.17). This configuration

49

2.8. VALIDATION TEST CASES

is very important in the context of computation of deformable particles under flow. It has been extensively studied in the zero Reynolds limit [17, 18, 20, 75, 111, 162] and now continues to serve as a validation case for numerical methods [49, 117, 124].

D

θ

kt

kt

Figure 2.17: Typical evolution of the deformation parameter D and orientation θ of a capsule versus time in simple shear flow. In addition to the tank-treading movement, Lac et al. [111] have shown another interesting feature: under some conditions, the membrane is subjected to compressive stress and buckles. This buckling may lead to instabilities in some numerical methods. The modeling of finite-thickness membrane capsules (Dupont et al. [58]) shows that the global capsule shape is quasi-identical whether the thickness of the membrane is modeled or not. The membrane thickness essentially impacts the appearance and the characteristics of wrinkles present in the buckling situation. To define a particular test case, one needs to introduce these several parameters:

• The viscosity ratio λ:

λ=

µint . µext

(2.66)

• The Reynolds number, in Eq. (2.67), is defined with the capsule radius a, the shear rate k associated to the prescribed flow and the external dynamic viscosity µext . It compares inertial effects to viscous effects at the scale of the capsule: Re =

ρka2 . µext

(2.67)

In all the test cases performed in this study, the Reynolds number is considered low, as in the literature.

• The capillary number is also introduced, as in Eq. (2.68) with Es , the shear modulus of the capsule membrane. It compares viscous forces to the membrane’s

50

CHAPTER 2. NUMERICAL METHOD elastic forces: Ca =

µext ka . Es

(2.68)

• In the case of a membrane with its mechanics defined with Skalak law, an important parameter is C, it is the ratio of the area dilation modulus Ea and the shear modulus Es : Ea C= . (2.69) Es Numerical configuration and parameters In numerous studies from the literature, the numerical domain is infinite, which ensures the absence of boundary effects on the dynamics of the capsule. In our numerical framework the domain is finite and meshed, thus the larger this domain, the less these effects affect the results, while increasing computational cost. The domain size ratio is defined as: L , (2.70) a with the cubic domain having a size of 2L. The mesh size ratio, defining the number of fluid mesh nodes on a characteristic length, is defined as (with a the radius of the capsule and dx the fluid mesh size): a . (2.71) dx All the results will be displayed in terms of a deformation parameter D and orientation of the capsule. D defined as: A−B , (2.72) A+B where A and B are the small and large axes of the ellipsoid that has the same matrix of moments of inertia as the one computed for the deformed capsule. Its asymptotic value, reached when the simulation has converged (i.e, when the shape is stabilized) will be particularly studied. In cases were buckling occurs, notably for small Ca [18], our simulations are unstable in the long run. Values are extracted from the simulation after the stabilization and before the instability, typically for kt ∼ 10. The orientation of the membrane is defined as the angle between the direction of the flow and the large axis of the membrane’s inertia ellipsoid in the shear plane. D=

Domain size convergence In this section, the effect of the domain size on the results is investigated. The initially spherical capsule is modelled with the Neo-Hookean Law. It is suspended in a domain for which the size is varied. The capsule deformation and orientation after stabilization are then reported. The variation of the value for the deformation or orientation of the capsule is represented by the following quantity: xcurrent − xmax × 100. (2.73) var(x) = xmax

51

2.8. VALIDATION TEST CASES Simulations parameters L a Ca λ dx a 6.25 0.075 1 4 − 20

Law Neo-Hookean

Table 2.3: Summary of the red blood cell parameters for the domain size convergence.

These

var(D) (%)

a = 25). With xmax being value of x in the extreme case (i.e La = 20 or dx Figure 2.18 and 2.19 present var(D) and var(θ) versus the domain size ratio.

L a

Figure 2.18: Percentage of relative variation of the deformation versus domain size ratio.

results show that the boundary effects are almost negligible when the domain size ratio is greater than 8.

52

var(θ) (%)

CHAPTER 2. NUMERICAL METHOD

L a

Figure 2.19: Percentage of relative variation of the orientation versus domain size ratio.

Mesh size convergence Simulations parameters L a Ca λ dx a 6 − 25 0.5 1 4

Law Neo-Hookean

Table 2.4: Summary of the red blood cell parameters for the mesh size convergence. In this section, the effect of the mesh size on the results is investigated. The initially spherical capsule is modelled with the Neo-Hookean law. Both mesh sizes for the fluid and membrane are the same, decreasing with each successive case. The deformation and orientation of the capsule after stabilization is then reported in figures 2.20 and 2.21.

53

var(D) (%)

2.8. VALIDATION TEST CASES

a dx

var(θ) (%)

Figure 2.20: Percentage of relative variation of the deformation versus mesh size.

a dx

Figure 2.21: Percentage of relative variation of the orientation versus mesh size. Figure 2.20 and 2.21 shows that when the fluid and membrane meshes are not refined enough, the associated variation in the results can be of the order of 5%. This indicates a = 10 to ensure acceptable that the membrane and fluid needs to be refined at least to dx accuracy of the results. A study of the over-refinement of the membrane with respect to the fluid refinement has been performed and a maximum variation of 0.5% has been

54

CHAPTER 2. NUMERICAL METHOD

found. In the next sections, all test cases will be performed in a fixed domain (Fig. 2.22 and a Fig. 2.23) with dx = 10 and La = 16.

Figure 2.22: Numerical domain chosen for the linear shear flow test cases.

Figure 2.23: Zoom on the chosen numerical domain. The membrane used in the following test cases is meshed with 2906 triangular elements as presented on Fig. 2.24.

55

2.8. VALIDATION TEST CASES

Figure 2.24: Membrane mesh in the linear shear flow test cases.

Neo-Hookean capsules in linear shear flow In this section, the quality of the results is investigated in the case of the Neo-Hookean Law for different values of the capillary number (Fig. 2.25). Our results are compared to Walter et al. [201], Li et al. [115], Lac et al. [111] and Luo et al. [124]. The simulations are performed in the limit of zero Reynolds number; as presented in Mendez et al. [132] the inertial terms are not computed in order to reduce the computation time. This approximation can create overshoots on the deformation and orientation of the cell during the transient part of the simulation. It is not modifying the static results that are reported here after stabilization of the cell. Simulations parameters a Ca λ dx 10 0.15 − 0.9 1

L a

16

Law NH

Table 2.5: Summary of the red blood cell parameters for the Neo-Hookean cases. A good agreement with the literature is found for the asymptotic deformation parameter of the cell. As expected, the tank-treading behaviour of the cell (Fig. 2.16) is observed in all theses cases. Some cases encountered did not reach convergence because of the absence of bending resistance, allowing the creation of unwanted folds on the extremities of the deformed capsules, the capillary number associated to these cases are in agreement with the range of stability (0.45 < Ca < 0.63) found in Walter et al. [201] .

56

Dasymp.

CHAPTER 2. NUMERICAL METHOD

Ca Figure 2.25: Asymptotic deformation versus capillary number in the case of the NeoHookean law. Skalak law capsules in linear shear flow In this section, the quality of the results is investigated in the case of the Skalak Law for different values of the capillary number (Fig. 2.26). Our results are compared to Walter et al. [201], Li et al. [115] and Lac et al [111] over the range of parameters presented in Table 2.6. Simulations parameters a Ca λ dx 10 0.15 − 1.8 1

L a

16

Law SK

Table 2.6: Summary of the red blood cell parameters for the Skalak law cases. Figure 2.26 shows a good agreement with the literature. The tank-treading behaviour of the cell is also observed in all the cases. A slight under estimation is found for the a a case where Ca = 1.8. In this one, the ratio dx had to be increased to dx = 16 in order to capture the large deformation of the cell. Further increase of this ratio did not reduce the under estimation.

57

2.8. VALIDATION TEST CASES

Figure 2.26: Asymptotic deformation versus capillary number.

Effect of the calculation of viscosity Using the Skalak law in a large domain, two ways of calculating the viscosity are compared: eight calculations are run, at viscosity ratios 0.2 and 5.0 for capillary numbers 0.6 and 1.5 as presented in Table 2.7. Simulations parameters a Ca λ dx 10 0.6 and 1.5 0.2 and 5.0

L a

16

Law SK

Table 2.7: Summary of the red blood cell parameters for the study of the effect of the calculation of the viscosity. The way of calculating the variable viscosity differs: we either use the classical approach, in which the discrete Dirac function is also used to update the field of viscosity. As in the front-tracking method presented in section 2.5, ensuring a smooth transition of the viscosity field in the membrane region. It will be referred as Smooth in the following tables. The second method is to determine if xn (a fluid grid node coordinates) is located inside or outside the capsule. The viscosity is thus imposed without smooth transition, with a jump through the membrane. In other words, the indicator function is here a Heaviside function. This method will be called the Jump method and referred as Jump in the following tables.

58

CHAPTER 2. NUMERICAL METHOD

Simulations parameters and results D Ca = 0.6, λ = Ca = 1.5, λ = 0.2 0.2 Smooth 0.4527 0.5755 Jump 0.4521 0.5735

Ca = 0.6, λ = 5.0 0.267 0.268

Ca = 1.5, λ = 5.0 0.282 0.281

Table 2.8: Results of the asymptotic deformation D for the study of the calculation of the viscosity.

Simulations parameters and results θ/π Ca = 0.6, λ = Ca = 1.5, λ = 0.2 0.2 Smooth 0.1588 0.1307 Jump 0.1587 0.1306

Ca = 0.6, λ = 5.0 0.0543 0.0548

Ca = 1.5, λ = 5.0 0.032 0.033

Table 2.9: Results of the asymptotic orientation θ/π for the study of the calculation of the viscosity. Results displayed in Tables 2.8 and 2.9 show that the method to compute the viscosity field has a very limited impact on the results in terms of deformation and inclination angles of the capsules. In the following, the jump method is used since it is computationally cheaper. Comparison with an equivalent numerical method and effect of the regularization procedure In order to estimate the effect of the regularization procedure on the results, the results of simulations for the case also used by Luo et al. [123] are reported. In that case, C = 1 and λ = 1. The size of the domain Lx × Ly × Lz is Lx = 8a, Ly = 8a, Lz = 5a in the paper of Luo et al. [123], it is kept the same here. A Cartesian mesh is used to compare the original regularization by Peskin (STRUCTURED) to the unstructured versions at the first order (RKPM1) and the second order (RKPM2). The type of regularization performed is varied, and results in terms of deformation and orientation are displayed in Figures 2.27 and 2.28. Simulations parameters a Ca λ dx 10 0.15 − 1.8 1

Law SK

Table 2.10: Summary of the red blood cell parameters for the regularization procedure comparison cases.

59

Dasymp

2.8. VALIDATION TEST CASES

Ca

θ π

Figure 2.27: Asymptotic deformation versus capillary number. Comparisons of results by Foessel et al. [75], Luo et al. [123] and YALES2BIO. YALES2BIO results are generated using three different regularizations procedures: the classical cosine function for STRUCTURED meshes, and the unstructured method at the first order (RKPM1) and second order (RKPM2).

Ca Figure 2.28: Asymptotic orientation versus capillary number. Comparisons of results by Foessel et al. [75], Luo et al. [123] and YALES2BIO.

60

CHAPTER 2. NUMERICAL METHOD

YALES2BIO results are generated with a rather coarse mesh: the number of points per radius is 10, while Luo et al. [123] present results for 8, 12 and 16. Using a classical cosine window, results are close to the ones obtained by Luo et al. using a/dx = 12. However, using the unstructured regularization procedure at the second order, the deformation parameter D at high capillary number gets much closer to the reference data by Foessel et al. [75]. The inclination angles are well predicted by the different methods. Note also that at low capillary numbers, the structured procedure and the RKPM1 procedure overestimate D, while using RKPM2 improves the results. Comparisons with the Foessel data at various viscosity ratio Simulations parameters a Ca λ dx 10 0.15 − 1.8 0.2 − 5

L a

16

Law SK

Table 2.11: Summary of the capsule parameters for the comparison with Foessel data cases.

Dasymp

Here, the ‘best’ setup is used as the viscosity ratio is varied: the domain is large, with L/a = 16 to ensure that the size of the domain has no influence on the results.

λ = 0.2 λ = 1.0 λ = 5.0 λ = 0.2 λ = 1.0 λ = 5.0

Ca Figure 2.29: Asymptotic deformation versus capillary number. Comparisons of YALES2BIO results with data by Foessel et al. [75], at three different viscosity ratios. In addition, an unstructured grid is used with a regularization at the second order.

61

2.8. VALIDATION TEST CASES

The Jump method is used to calculate the viscosity (sudden jump in the viscosity of the fluid is imposed at the membrane location between the internal and the external values).

θ π

λ = 0.2 λ = 1.0 λ = 5.0 λ = 0.2 λ = 1.0 λ = 5.0

Ca Figure 2.30: Asymptotic orientation versus capillary number. Comparisons of YALES2BIO results with data by Foessel et al. [75], at three different viscosity ratios. Results are in good agreement with the reference data generated by Foessel et al. [75]. The largest differences are obtained at λ = 5.0: the deformation parameter is underestimated compared to the reference data. However, inclination angles are well reproduced. Except for D at λ = 5.0 and for Ca higher than 0.6, where errors are of the order of 3%, the agreement is excellent.

2.8.5

Red blood cells in simple shear flow

As a canonical configuration to understand the behavior of red blood cells under flow, the case of an isolated red blood cell flowing in an unconfined pure shear flow has been extensively studied over the last decades. However, in spite of numerous theoretical, experimental and numerical studies on the topic, this simple configuration continues to reveal new aspects of the red blood cell dynamics. Behavior of a red blood cell in a simple shear flow As it is relatively easy to perform experimentally, this test case has been investigated extensively over the years. Since the 1970s, two types of motion have been known to occur: the tumbling motion was observed for instance by Goldsmith and Marlow [84],

62

CHAPTER 2. NUMERICAL METHOD

where the red blood cell flips in the flow like a rigid body. Another regime was also identified, where the red blood cells ‘oriented themselves at a constant angle to the flow and their membrane appeared to rotate about the interior’, as stated by Goldsmith and Marlow [84], analogously of a droplet. This regime, referred to as tank-treading to describe the movement of the membrane about the cell interior, was especially visible in viscous Dextran solutions. This movement was then observed in details by Fischer and collaborators [70, 71, 74], in particular. Their detailed visualizations of tank-treading red blood cells were performed in a viscous suspending medium, typically twenty to forty times the viscosity of plasma [74]. The biconcave shape of the red blood cell was observed to be conserved for low and moderate shear rates. For high shear rates, the red blood cell strongly deforms to an ellipsoid. The dynamics of red blood cells was notably shown to depend on the viscosity ratio between the internal medium and the suspending fluid. In 1982, Keller & Skalak [104] published an analytical model, able to predict either a steady orientation or a tumbling regime depending on the viscosity ratio, notably, but failed to predict the tumbling to tank-treading transition as a function of the shear rate [3]. The case of high viscosity of the suspending medium and moderate shear rates has been investigated in the 2000s, in particular to study the tumbling to tank-treading transition, for shear rates of 1.0 s−1 approximately. By doing so, and thanks to sideview microscopic imaging, a new motion has been unraveled: the swinging, where the orientation of the cell in the tank-treading motion is seen to oscillate around a mean value [3]. It has been described by analytical models by Skotheim & Secomb [175] and Abkarian, Faivre & Viallat [3] by adding the contribution of the elasticity of the cell membrane to the Keller & Skalak model (more details in Chapter 3). These models were able to predict the amplitude of the oscillation of the orientation of the cell and the transition between the tank-treading and tumbling motions as a function of the shear rate. Numerically, the challenge has been to take all the mechanical effects into account. When the viscosity of the cytoplasm (inner fluid), bending and shear elasticity were accounted for, the swinging motion and transitions between motions were retrieved. However, until recently, the shear rates for the transition were overestimated. Recent studies showed the importance of the stress-free shape of the cell on these results [147]: the biconcave stress-free shape implied large deformations at the transition between tank-treading and flipping, as well as inaccuracies in retrieving the critical values of the shear rate of tumbling to tank-treading transition. The recent success of numerous groups in retrieving experimental results has led to community to adopt the following strategy for the a priori unknown cytoskeleton stress-free shape: in numerical simulation a spheroidal stress-free shape is desired. It cannot be a sphere, as there would be no shape memory [71], but it should be close to a sphere, typically an oblate ellipsoid of volume ratio higher than 0.95 [38, 147, 189]. This choice enables to better recover the tumbling to tank-treading transition, predict

2.8. VALIDATION TEST CASES

63

the swinging motion and also to maintain the biconcave shape under moderate shear stresses. It also improves predictions of analytical models derived from the Keller & Skalak model [56]. More recently, additional experiments for physiological viscosity of the suspending medium revealed that tumbling, swinging and tank-treading were not the only regimes present for a red blood cell. In particular the off-shear plane dynamics was studied, and Dupire, Socol and Viallat [57] showed that when increasing the shear rate, the tumbling motion is unstable and the movement changes to rolling, where the red cells align their small axis with the vorticity direction. In the rolling motion, the cell is almost undeformed. Additional movements of so-called kayaking/hovering/frisbee-like movements were also reported at higher shear rates [36, 38, 57] The difficulty of imaging flowing cells at high velocity explains that the regime at high shear stresses and physiological viscosity ratio has not been studied until 2014-2015, by Lanotte and Abkarian [4]. New visualizations show that the dynamics revealed by Dupire et al. does not fully describe the movements of red blood cells. Lanotte and Abkarian [4] notably show that when red blood cells are suspended in plasma (or a fluid of similar viscosity), at shear rates higher than a few hundreds, polylobed complex shapes are observed, with a typical sequence of stomatocytes, trilobe, tetrahedra, startetralobes, etc. These shapes have also been observed in numerical simulations by Mauer, Fedosov and Gompper and by Mendez (using YALES2BIO). Publications are in preparation. Numerical simulations The dynamics of red blood cells in a simple unconfined shear flow is thus extremely complex and a full description of this flow is not the objective of this thesis. Only a few cases are presented to show the ability of YALES2BIO to reproduce typical dynamics, in the case where the viscosity of the suspending medium is higher than that of the cytoplasm. The computational domain is shown in Fig. 2.31a). It consists of a cube of edge L, with the origin at its center, closed by two sliding walls at y = ± L/2 and by periodic boundary conditions in the other directions. x is the direction of the main flow, y is the direction of the velocity gradient and z the vorticity direction. In the absence of red blood cell, the flow is linear ~v = ky e~x . At the beginning of the simulations, a cell is deposited in its resting shape, with its center of mass at the origin and its small axis aligned with the x direction. The Reynolds number is 0.1 in all the simulations, and the viscosity ratio (internal viscosity µint divided by the external viscosity µext ) is always 0.2, the viscosity inside the cell being 5 times smaller than the external viscosity. Results will be displayed in terms of the capillary number, defined as Ca = µext ka/Es , where a is an equivalent radius of the particle (a = 2.82 µm) and Gs the shear modulus. Here, the in-plane elasticity is described by the Skalak model, with C = 30, which is sufficient to ensure a good

64

CHAPTER 2. NUMERICAL METHOD

conservation of the total area of the red blood cell. The bending model by Helfrich is used, with Eb = 6.0×10−19 J. A positive spontaneous curvature is imposed, of C0 = 1.41 µm−1 . The stress-free shape is an ellipsoid with volume ratio 0.99. Values are identical to the ones chosen by Cordasco et al. [38] in their analysis of the effect of stress-free shapes on RBC in shear flow. The fluid domain is discretized by a Cartesian grid of 603 elements (L/a is of the order of 7.5, which yields a resolution of 0.35 µm). The membrane is discretized with a similar resolution (approximately triangular 2000 faces). b)

a)

y

x x

L

lx

z

z

lz Figure 2.31: Computational configuration of interest. Computational domain definition (a) and definition of the extensions of the cell in the tank-treading regime (b). Typical dynamics is displayed for three cases in Fig. 2.32, at capillary numbers a) Ca = 0.02 , b) Ca = 0.04 and c) Ca = 0.28, for which tumbling, swinging and tank-treading (with very small swinging oscillations) are observed. The regimes are consistent with the observations of Abkarian, Faivre & Viallat [3] (in a suspending medium of viscosity µext = 0.031 Pa.s for instance, Ca = 0.02, 0.04 and 0.28 correspond to k = 0.6 s−1 , 1.2 s−1 and 8.0 s−1 , respectively. At higher shear rates, the red blood cell tank-treads. Figure 2.31b shows how simulations are compared to experimental data: lengths in the flow (lx ) and vorticity (lz ) directions are measured and a projected deformation rate D = (lx − lz )/(lx + lz ) is calculated. This quantity is compared to experimental data from direct visualizations and ektacytometry [70, 87, 187] in Fig. 2.33. The experimental capillary numbers are calculated by assuming the same value for the shear modulus as in the computations, Gs = 2.5 µN.m−1 . As for swinging, small oscillations of D are observed over time, so only the average value is displayed. Good agreement with experimental data is shown,

65

2.8. VALIDATION TEST CASES a) Re = 0.1, λ = 0.2, Ca = 0.02. Images every kt = 1.25.

b) Re = 0.1, λ = 0.2, Ca = 0.04. Images every kt = 2.5.

c) Re = 0.1, λ = 0.2, Ca = 0.28. Images every kt = 2.5.

Figure 2.32: a) Tumbling, b) swinging and c) tank-treading of red blood cells, shown by a sequence of shapes in calculations. 0.6

0.5

0.4

D

0.3

0.2

Fischer 1980 Tran-Son Tay et al. 1984 young cells Tran-Son Tay et al. 1984 old cells

0.1

Hardeman et al. 1994 YALES2BIO

0.0 7

8

9

2

3

4

5

1

Ca

Figure 2.33: Projected deformation index D of tank-treading red blood cells as a function of the capillary number Ca . Numerical results are compared with data from Fischer (1980) [70], Tran-Son-Tay et al. (1984) [187] and Hardeman et al. (1994) [87].

albeit with over-estimations of the deformation at high capillary numbers. In conclusion, the numerical method introduced in this chapter is able to simulate the dynamics of deformable particles, particularly red blood cells in flow at arbitrary Reynolds number, which is an important feature since it is often overlooked in favor of simulations in the physiological context. The solver has been validated against several test cases, assessing of the quality of the results in dynamical and stationary setups.

Chapter

3

Analytical modeling of the dynamics of red blood cell in elongational flow Chapter contents 3.1

3.2

Derivation of the Keller & Skalak model . . . . . . . . . 3.1.1 Formulation of the problem . . . . . . . . . . . . . 3.1.2 External flow and equilibrium . . . . . . . . . . . . 3.1.3 Internal flow and conservation of energy . . . . . . 3.1.4 Abkarian, Faivre & Viallat’s variation . . . . . . . 3.1.5 Planar elongational flow derivation . . . . . . . . . 3.1.6 Jeffery’s theory . . . . . . . . . . . . . . . . . . . . Models behavior in the case of a planar elongational flow 3.2.1 Parameters . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Effect of the capillary number . . . . . . . . . . . . 3.2.3 Effect of the viscosity ratio . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

68 68 70 73 75 77 78 79 79 80 84

Since the early twentieth century, a number of studies focused on the modeling of the dynamics of the orientation of solid ellipsoidal particles in various configurations. Notably, Jeffery [99] opened the way with an insightful paper published in 1922 describing the orientation of a solid ellipsoidal particle in general flows. Roscoe [165] later extended the scope of this paper to viscoelastic particles suspensions with surface motion, which is one of the field of application of Jeffery’s theory nowadays. Keller & Skalak [104] introduced a model based upon Jeffery’s and Roscoe’s works for modeling the dynamics of red blood cells in a linear shear flow, using the model of a fixed-shape ellipsoidal particle with an imposed membrane velocity, reproducing the well-known tank-treading movement witnessed in experiments [70, 72, 74, 187]. While the Keller & Skalak model is able to reproduce the tank-treading and the tumbling movements and their transition with the internal-to-external viscosity ratio, the model results do no 67

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 68 CELL IN ELONGATIONAL FLOW depend on the shear rate, while a shear rate dependency is obtained experimentally. In 2007, two groups proposed a modification of the Keller & Skalak model, introducing the elasticity of the membrane [3, 175]. This improvement enabled to recover the tumblingto-tank-treading transition with the shear rate. The framework of Jeffery’s theory is general and several papers proposed a new analytical derivation of it, such as the work of Jeˇzek [100] or Junk [102]. However, models relative to red blood cells were only derived for the linear shear flow. Other models, such as the one of Vlahovska et al. [198], consider non-fixed shaped particles without an imposed membrane velocity and could have been used to model the dynamics of red blood cells in the blood analyzer. However, since these models consider that the cell is quasi-spherical, which is a strong approximation on the shape of the red blood cell, they have not been used in this present work. This section aims at adaptating the Keller & Skalak model and the extension of Abkarian, Faivre & Viallat in the case of a hyperbolic flow which is typical to the velocity field red blood cells are undergoing in the blood analyzer. This work relies on the paper of Keller & Skalak [104] and Abkarian, Faivre & Viallat [3], but using a hyperbolic flow as the carrying flow. The derivation of the Keller & Skalak model and the extension of Abkarian, Faivre & Viallat are first presented in the case of an axi-symmetric elongational flow and thoroughly explained, as it has been found that some analytical tricks and expressions were not detailed in the original paper of Keller & Skalak. Then, the models are investigated for an ellipsoidal particle in a planar elongational flow.

3.1

Derivation of the Keller & Skalak model

In order to facilitate the comparison with the Keller & Skalak model derived for pure shear flow, the same structure is used for presenting the analysis, and notations are conserve as much as possible.

3.1.1

Formulation of the problem

An ellipsoidal membrane containing a viscous fluid is subjected to an elongational flow. As in the seminal work of Keller & Skalak, the shape of this ellipsoidal membrane is considered fixed and axes are defined as in Fig. 3.2, where xˆi denote coordinates in a fixed Cartesian coordinate system, and xi denote coordinates in a second Cartesian coordinate system of same origin, whose axes correspond to the principal axes ai of the ellipsoidal shape and referred to as the body frame. xˆ3 and x3 are assumed to coincide, when the x1 and x2 axes are rotated by an angle θ with respect to the xˆ1 and xˆ2 axes. Angle θ is positive for a counter-clockwise rotation and depends on time, thus, the angular velocity θ˙ is defined as: dθ θ˙ = , (3.1) dt

69

3.1. DERIVATION OF THE KELLER & SKALAK MODEL xˆ2

x2

µ x1

n

θ a2

a1 xˆ1

E

µ′

∂E Figure 3.1: Schematic of an ellipsoidal membrane filled with fluid of viscosity µ′ suspended in a fluid of viscosity µ, with the definition of the coordinate frames used in the analytical developments.

1

The quantity a0 = (a1 a2 a3 ) 3 is introduced for the ellipsoidal membrane of semi-axes ai , allowing the definition of αi as: αi =

ai , a0

(3.2)

The velocity of the fluid external to the particle relative to the fixed frame and expressed in the body frame is denoted by vˆi . Alternatively, the velocity relative and expressed in the body frame is denoted by vi . The components of the velocity of the internal fluid are denoted by vˆi ′ and vi′ . The flow is supposed to be governed by Stokes equations and the continuity equation, thus: µvi,jj = p,i

vi,i = 0,

(3.3)

′ µ′ vi,jj = p′,i

′ vi,i = 0,

(3.4)

where p and p′ are the pressure fields and µ and µ′ are the viscosity associated with the external fluid and the fluid inside the membrane, respectively. The notation of the Keller & Skalak paper is used, where a comma followed by subscript i denotes the derivative with respect to xi . The undisturbed velocity field (far from the particle) relative to the fixed frame is an axi-symmetric elongational flow, thus:   uˆ1 0 = κxˆ1             

uˆ2 0 = − κ2 xˆ2 ,

uˆ3 0 = − κ2 xˆ3

(3.5)

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 70 CELL IN ELONGATIONAL FLOW with κ the extensional rate. Referred to the body frame this undisturbed velocity field writes:  i  h   0 3 2 θ − 1 sin2 θ + x  sin 2θ − v ˆ = κ x cos 2 1 1  2 4      

h







vˆ2 0 = κ x1 − 34 sin 2θ + x2 sin2 θ − 12 cos2 θ        

i

,

(3.6)

vˆ3 0 = − κ2 x3

Keller & Skalak [104] define the membrane surface velocity vim relative to and referred to the body frame as, v1m

a1 = −ω˙ x2 , a2 



v2m

a2 = ω˙ x1 , a1 



v3m = 0,

(3.7)

with ω˙ the tank-treading frequency. The components of the unit outer normal to ∂E (the surface of the particle) in the body frame write ni =

xi a2i

!

x2j a4j

!− 1 2

(sum on j, no sum on i),

(3.8)

which implies vim ni = 0,

(3.9)

corresponding to the fact that the surface velocity is everywhere tangent to ∂E. As mentioned in [104], since v3m = 0, material surface points move along elliptical paths in planes parallel to the (x1 ,x2 )-plane. For positive values of ω, ˙ the membrane motion is in the counter-clockwise direction when viewed from the positive x3 axis. The boundary conditions imposed on the internal and external fluids are:   vi = vim      

3.1.2

      

on

∂E

vi → vi0

as

|x| → ∞

vi′ = vim

on

∂E

(3.10)

External flow and equilibrium

The components of the surface-stress vector referred to the body frame are introduced as   (3.11) Ti = −p′′ ni + µ A∗ij + 2em ij nj , where p′′ is an arbitrary constant pressure and the strain rate em ij is defined as, em ij =

 1 m m . vj,i + vi,j 2

(3.12)

The first term of Eq. (3.11) represents the surface stress created on the particle by the surrounding pressure field. The second term gathers the effect of the velocity field on

3.1. DERIVATION OF THE KELLER & SKALAK MODEL

71

the particle through the strain rate em ij and a tensor Aij ∗ that depends on the shape of the particle. Introducing the following geometrical integrals only depending on the shape of the membrane, Z ∞  ds   g = 1  2  (α1 + s)∆  0       Z ∞  ds

g′ =

1  (α22 + s)(α32 + s)∆ 0       Z ∞   sds    g1′′ = 2 + s)(α2 + s)∆ (α 0 2 3

with,

∆2 = (α12 + s)(α22 + s)(α32 + s),

(3.13)

(3.14)

the tensor element A∗11 is defined as A∗11 =

4 2g1′′ e∗11 − g2′′ e∗22 − g3′′ e∗33 . 3 g2′′ g3′′ + g3′′ g1′′ + g1′′ g2′′

(3.15)

Note that in the original paper of Keller & Skalak [104], the expression for A∗11 is erroneous, with gi in the numerator instead of the present gi′′ as found in Jeffery [99], Roscoe [165] and Jeˇzek [100]. The expressions for A∗22 and A∗33 can be obtained by cyclic permutations. The tensor elements A∗12 and A∗21 read A∗12 = 4

∗ −ǫ g1 e∗12 − α22 g3′ (ζ12 12k ωk ) , ′ 2 2 g3 (α1 g1 + α2 g2 )

(3.16)

A∗21 = 4

∗ −ǫ g2 e∗21 − α12 g3′ (ζ21 21k ωk ) . ′ 2 g3 (α1 g1 + α22 g2 )

(3.17)

with ǫijk the alternating tensor (Levi-Civita symbol). For i 6= j, the other elements of A∗ij can be obtained by permutations. The quantities in all elements of A∗ij are defined as  1 0 0 0 (3.18) v ˆ + v ˆ e∗ij = e0ij − em , e = i,j , ij ij 2 j,i and, m 0 ∗ , (3.19) − ζij = ζij ζij with,

  1 m 1 0 0 m 0 . ζij = , (3.20) vˆj,i − vˆi,j vˆj,i − vˆi,j 2 2 which will be detailed later in the case of the hyperbolic flow. Finally, ωk is the angular velocity of the ellipsoid with respect to the fixed frame. For the case considered, the following values for ωk are straight-forward: m ζij =

ω1 = 0,

ω2 = 0,

˙ ω3 = θ.

(3.21)

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 72 CELL IN ELONGATIONAL FLOW The expression for Ti in Eq. (3.11) is valid for any undisturbed velocity flow and membrane velocity vim that are both linear in xi . Any angular velocity ωi can also be considered. The stress vector Ti may be written as 

cij = −p′′ δij + µ A∗ij + 2em ij

Ti = cij nj , with



(3.22)

The ellipsoidal particle membrane undergoes a force exerted by the external fluid, which components in the body frame are Fi =

Z

∂E

cij nj dA.

(3.23)

Fi = 0 as the contour ∂E is closed and cij are constant with respect to xi . The circulation of the normal nj is obviously always null. A moment about the origin is exerted by the external liquid on the ellipsoidal membrane. At any instant, the components of this resultant moment in the body frame are Mi = Also, note that

Z

∂E

Z

∂E

ǫijk xj ckl nl dA.

xi ni dA = V

where V is the volume of the ellipsoid. Since

(no sum), Z

∂E

state that

(3.24)

(3.25)

xi nj dA = 0 if i 6= j, it is possible to

M3 = V (c21 − c12 ).

(3.26)

The considered hyperbolic flow, Eq. (3.5), implies the following results for e∗ and ζ ∗ :  

κ(cos2 θ − 12 sin2 θ)

3 e∗ =  − 4 κ sin 2θ − 0



ω˙ 2



 

ω˙ ζ∗ =  2

a22 −a21 a1 a2



0 a22 +a21 a1 a2

0



− 34 κ sin 2θ − κ(sin2 θ −

2 2 ω˙ a2 −a1 2 a1 a2 1 2 2 cos θ)





a22 +a21 a1 a2

0



0

0

0

 

 

0  

− κ2

0

− ω2˙



(3.27)

0 

(3.28)

,

(3.30)

0

Using the fact that αi2 gj′ + gj′′ = gk with i 6= j 6= k and substituting the expressions of e0 , em , e∗ , ζ 0 , ζ m and ζ ∗ for the carrying hyperbolic flow in c12 and c21 , it may be shown that, M3 = M3H + M3F + M3T , (3.29) with,   M3H = −C 34 κ sin 2θ(a21 − a22 )             

˙ 2 + a2 ) M3F = −C θ(a 1 2 M3T = −C2a1 a2 ω˙

3.1. DERIVATION OF THE KELLER & SKALAK MODEL and, C=

4µV . + a22 g2

73

(3.31)

a21 g1

M3H corresponds to the effect of the external hyperbolic flow on the stationary particle at angle θ. M3F is the moment acting on the rigid ellipsoid flipping around the x3 axis. M3T is the moment on the ellipsoidal particle undergoing a possible tank-treading motion. The membrane is considered neutrally buoyant, so that no forces or couples are exerted and that inertial effects can be neglected. Each fluid element inside the membrane is thus in equilibrium, which implies that the whole membrane is also in equilibrium at each instant of time. Thus, the resultant force and moment exerted by the external fluid must vanish at each instant, which leads to M3 = M3H + M3F + M3T = 0. Substituting the values for M3H , M3F and M3T in the previous equation yields: ˜ sin 2θ, θ˙ = A˜ + B with

2a1 a2 A˜ = − 2 ω˙ a1 + a22

(3.32) 2

2

˜ = − 3 a1 − a2 κ. B 4 a21 + a22

(3.33)

It is interesting to note that only the expression for M3H changes when using a hyperbolic carrying flow instead of the linear shear flow, M3S in Keller & Skalak’s analysis. It is justified by the fact that this term is associated with the moment due to the shear flow in the reference paper. It results in the presence of a sine term instead of a cosine term ˜ in Eq. (3.32). and a modification of factors A˜ and B

3.1.3

Internal flow and conservation of energy

The surface velocity from Eq. (3.7) is extended as the solution for the internal velocity field vi′ :     a1 a2 v1′ = −ω˙ x2 , v2′ = ω˙ x1 , v3′ = 0. (3.34) a2 a1 The internal velocity field from Eq. (3.34) satisfies Eq. (3.4) and the last line of Eq. (3.10). The dissipation function corresponding to vi′ writes: Φ′ = µ′ f1 ω˙ 2 , where, f1 = (r2 − r2−1 )2 ,

r2 =

(3.35) a2 , a1

(3.36)

Since Φ′ is considered homogeneous in space, the rate D′ at which energy is dissipated in the internal fluid is defined as: D′ = V µ′ f1 ω˙ 2 .

(3.37)

The reader is referred to the paper of Keller & Skalak [104] for a discussion about the energy dissipation of the membrane which may be even greater than the internal fluid

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 74 CELL IN ELONGATIONAL FLOW dissipation. For the sake of simplicity and clarity, the membrane energy dissipation is considered to be zero here, even though it is possible to add the contribution of the membrane by increasing the inner viscosity µ′ to a greater value µapp . The total rate at which the work exerted by the external flow over the whole particle is done can be expressed as Z Wp =

∂E

cij nj vˆim dA.

(3.38)

Thanks to the equilibrium condition (M3 = 0), it is possible to write Eq. (3.38) as: Wp =

Z

∂E

cij nj vim dA.

(3.39)

Starting from Eq. (3.39), we compute: Wp =

Z

c11 n1 v1m + c12 n2 v1m + c13 n3 v3m +

(

∂E

(3.40)

c21 n1 v2m + c22 n2 v2m + c23 n3 v3m + c31 n1 v3m + c32 n2 v3m + c33 n3 v3m ) dA.

Since v3m = 0, the last line of Eq. (3.40) is null. Also, from the expressions of e0 , em , e∗ , ζ 0 , ζ m and ζ ∗ , it is easy to find that c23 = c13 = 0, so that: Wp =

Z

∂E

(c11 n1 v1m + c12 n2 v1m + c21 n1 v2m + c22 n2 v2m ) dA,

(3.41)

From the expression of v1m and v2m , it yields that the terms in c11 and c22 are zero, since Z ∂E

ni xj dA = 0. Thus,

Wp = Using the fact that g3′ = Wp = µV ω˙

2

a21 − a22 a1 a2

g2 −g1 , α21 −α22

!2

Z

∂E

(c12 n2 v1m + c21 n1 v2m ) dA,

(3.42)

the following expression is found:

2µV ω˙ 2 − ′ 2 g3 (α1 + α22 )

a21 − a22 a1 a2

!2

3µV ωκ ˙ sin 2θ + ′ 2 g3 (α1 + α22 )

a21 − a22 a1 a2

!

, (3.43)

which can be rewritten as: Wp = µV



3 f2 ω˙ − f3 κω˙ sin 2θ 2 2



(3.44)

with 1 z1 = (r2−1 − r2 ) and z2 = g3′ (α12 + α22 ) 2 2 4z1 f2 = 4z12 (1 − ) and f3 = − . z2 z2

(3.45) (3.46) (3.47)

75

3.1. DERIVATION OF THE KELLER & SKALAK MODEL

Thus, the work of the external force must equal the internal dissipation at equilibrium (Wp = D′ ), we obtain a quadratic equation in ω, ˙ which is satisfied by either ω˙ = 0 or ω˙ =

3f3 κ sin 2θ , 2(f2 − λf1 )

(3.48)

Using this value of ω, ˙ Eq. (3.32) finally writes: θ˙ = A˜ sin 2θ,

(3.49)

with, !

a1 a2 f3 3(a21 − a22 ) κ A˜ = − 2 3 . + (f2 − λf1 ) 4 (a1 + a22 )

(3.50)

This differential equation can be solved analytically, resulting in: ˜

θ(t) = arctan(tan(θ0 )e2At )

(3.51)

It is obvious that θ = 0◦ is a stable solution for θ˙ = 0, the case of θ = 90◦ is an unstable equilibrium state since θ˙ is negative for θ < 90◦ and positive for θ > 90◦ thus bringing the capsule to θ = 180◦ over time, which is qualitatively equal to θ = 0◦ due to symmetries. Also, the tank-treading frequency ω˙ can be written as, ˜

3f3 κ sin(2 arctan(tan(θ0 )e2At )) . ω˙ = 2(f2 − λf1 )

3.1.4

(3.52)

Abkarian, Faivre & Viallat’s variation

The extension of Abkarian, Faivre & Viallat introduced in their recent work [3] is used to complete the model. It adds the elastic contribution of the membrane to energy equation. This extension enabled to recover the swinging motion witnessed in experiment, together with other complex features as an intermittent swinging-tumbling movements. In concrete terms, the balance of energy is modified by a term representing the elastic storage of energy in the membrane: Wp = D′ + Pel = V µ′ f1 ω˙ 2 + Pel

(3.53)

with, Pel =

Z

Ω

T r(σ :

em ij )dΩ

1 = ω˙ 2

a22 − a21 a1 a2

!2

Gs sin 2ω Ω,

(3.54)

with Ω, the volume of the cell membrane, and Gs the shear modulus. The shear stress tensor σ is defined as for a simple linear 3D thin visco-elastic with the Kelvin-Voigt model. The same idea was also pursued by Skotheim and Secomb [175] at the same time. There, the elastic energy term was estimated globally, without a calculation from a specific law. More advanced mechanical behaviors would add complexity to the

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 76 CELL IN ELONGATIONAL FLOW calculation, whereas the simple Kelvin-voigt model already yields very good results [3, 56]. Hence, σ = 2µm em (3.55) ij + 2Gs E. with E the Euler-Almansi strain tensor obtained from the displacements induced by the membrane velocity and µm the membrane viscosity. The shear stress tensor σ can be written as the sum of a fluid and an elastic part: σ = σel + σf luid , corresponding to the contributions of the lipid bilayer and the cytoskeleton in the case of a red blood cell. The fluid contribution can be written as, σf luid = −pI + 2µem ij .

(3.56)

where p is the pressure in the membrane. Since the membrane is considered as an isotropic, incompressible, homogeneous and linear elastic material, the elastic contribution to the strain tensor can be expressed using the Piola-Kirchhoff stress tensor: πel = π 0 I + 2Gs e, where π 0 is the initial stress of the membrane. π 0 is supposed to be 0, as in the initial work of Abkarian et al. [3], thus assuming the ellipsoidal shape is the stress-free shape. e is the strain tensor, which can be computed from the trajectory of a point on the membrane xk : a1 0 x sin(ω) a2 2 a2 x2 (t) = x02 cos(ω) + x01 sin(ω) a1 0 x3 (t) = x3 x1 (t) = x01 cos(ω) −

(3.57)

where ω corresponds to the phase angle of a membrane element and x0k are the initial vectorial positions. Thus, ! 1 ∂xk ∂xk eij = (3.58) − δij . 2 ∂x0i ∂x0j Therefore, πel is simply expressed using an infinitesimal tangential transformation F ≡ ∂xi : ∂x0 j

σel =

T

F−1 : 2Gs e : F−1 .

(3.59)

with, cos ω − aa12 sin ω 0   F =  aa21 sin ω cos ω 0 0 0 1 



(3.60)

Thus, the shear stress tensor σ writes:

T −1 σ = −pI + 2µem : 2Gs e : F−1 ij + F

(3.61)

From Eq. (3.61), it is easy to compute the integrand T r(σ : em ij ) of Eq. (3.54): m m T −1 T r(σij : em : e : F−1 : em ij ) = 2µT r(eij : eij ) + 2Gs T r( F ij )

(3.62)

77

3.1. DERIVATION OF THE KELLER & SKALAK MODEL

which when integrated, results in the right hand-side of Eq. (3.54). Using Eq. (3.53) and the expression for Pel of Eq. (3.54), the model turns into two differential equations: −2a1 a2 ω˙ 3(a21 − a22 ) θ˙ = 2 − κ sin 2θ, a1 + a22 4(a21 + a22 )

(3.63)

and, f3 κ 3 f1 Es Σ ω˙ = sin 2θ + sin 2ω , (f2 − λf1 ) 2 2f3 µκV 



(3.64)

where Σ is the surface of the membrane. It is easy to see that if the surface shear modulus is null (Es = 0), Eq. (3.48) is retrieved.

3.1.5

Planar elongational flow derivation

In this section, the results of the analysis in a planar elongational flow are presented for the sake of completeness. The carrying flow is considered two-dimensional for the sake of simplicity. The previous models are thus derived again following the previously presented method, only applied to a velocity field defined as:   u ˆ0 = κˆ x             

Keller & Skalak

vˆ0 = −κˆ y

(3.65)

w ˆ0 = 0

The Keller & Skalak model is derived once again, the development in 2D is straightforward and similar to the 3D case, as only the terms e∗ and ζ ∗ change. It finally results in, θ˙2D = A2D + B2D sin 2θ2D , (3.66) with, A2D = −

2a1 a2 ω˙ 2D a21 + a22

B2D = −

a21 − a22 κ. a21 + a22

(3.67)

where, f3 κ sin 2θ2D , (f2 − λf1 )

(3.68)

θ˙2D = A˜2D sin 2θ2D ,

(3.69)

ω˙ 2D =

The equation for θ˙2D can thus be written as:

with, A˜2D =

2a1 a2 f3 κ − − (a21 − a22 ) . 2 2 (f2 − λf1 ) (a1 + a2 ) 



(3.70)

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 78 CELL IN ELONGATIONAL FLOW Note that the difference in the evolution equation for θ between the 3-dimensional case, in Eq. (3.49), and the 2-dimensional case, in Eq. (3.69), is small. Again, this differential equation can be solved analytically, resulting in: ˜

θ2D (t) = arctan(tan(θ0 )e2A2D t )

(3.71)

Abkarian, Faivre & Viallat Once again, the Abkarian, Faivre & Viallat variation can be introduced in the Keller & Skalak model for a 2D hyperbolic flow, without taking the membrane viscosity into account, i.e µm = 0. (a21 − a22 ) −2a1 a2 ω˙ − κ sin 2θ, (3.72) θ˙2D = 2 a1 + a22 (a21 + a22 ) and, ω˙ 2D =

f3 κ f1 Es Σ sin 2ω , sin 2θ + (f2 − λf1 ) 2f3 µκV 



(3.73)

with Σ the surface of the cell membrane and µ the external viscosity. Again, if the surface shear modulus is null (Es = 0), Eq. (3.69) is retrieved.

3.1.6

Jeffery’s theory

This model has been introduced by Jeffery [99]. It is here presented in a simplified two-dimensional version. It is used to predict the orientation of a rigid ellipsoid in a Stokes flow. It differs from the previous models as it does not take a membrane velocity or an elastic contribution into account. The direction of the cell is represented by the orientation vector p~, which is identical to the x1 vector from Fig. 3.2, in the xˆ1 − xˆ2 plane: ! cos θ p~ = . (3.74) sin θ with θ the angle between p~ and the xˆ1 -axis. From the fluid flow, the variation of this vector can be calculated using: ∂pi = Rij pj + D(Sij pj − Skl pk pl pi ), ∂t

(3.75)

with, a2 − a22 , D = 21 a1 + a22

R=

!

0 0 , 0 0

S=

!

κ 0 . 0 −κ

(3.76)

with R and S, the antisymmetric and symmetric parts of the velocity gradient tensor, respectively. D is the aspect ratio of the particle in the plane. If Eq. (3.75) is expanded with the values for S in the case of a 2D hyperbolic flow, the following expression is easily found: a2 − a22 sin 2θ, (3.77) θ˙ = −κ 21 a1 + a22

3.2. MODELS BEHAVIOR IN THE CASE OF A PLANAR ELONGATIONAL FLOW

79

which corresponds to the Keller & Skalak and Abkarian, Faivre & Viallat models with particular hypotheses on the viscosity ratio or capillary number, as explained in the next section. Jeffery’s model is easily obtained from the K&S model and the AFV model by imposing ω˙ = 0 in Eq. (3.66).

3.2

Models behavior in the case of a planar elongational flow

The behavior of the previously derived models are presented in the case of a particle deposited in a 2-dimensional hyperbolic flow with a constant extensional rate κ. From its initial orientation θ = 45◦ , it is expected to orient so that its longer axis is aligned with the direction of extension, defined as θ = 0◦ . The influence of two parameters are investigated: the capillary number Ca , which influences the Abkarian, Faivre & Viallat model through the surface shear modulus and the viscosity ratio λ, which is accounted for in both models. The Jeffery model is also used.

3.2.1

Parameters

In order to define the case considered, several parameters are introduced: • The capillary number associated to the particle Ca =

µκa1 , Es

(3.78)

with Es the surface shear modulus of the capsule membrane. Ca is the ratio of the viscous effects over the mechanical elastic effects. This parameter will vary between the different cases presented, as it is a parameter of the Abkarian, Faivre & Viallat model through the surface shear modulus Es . • The viscosity ratio:

λ=

µint , µext

(3.79)

which is the ratio of the dynamic viscosity between the internal and external fluids. It will also be varied between the different cases since it is a parameter of all the models derived previously. The values and ranges of all the parameters is summarized in table 3.1.

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 80 CELL IN ELONGATIONAL FLOW Particle parameters Axes ratio a1 /a2 Axes ratio a2 /a3 Deformation parameter D = Membrane surface Capillary number Viscosity ratio Initial angle

a21 −a22 a21 +a22

a1 /a2 = 2 a2 /a3 = 0.5 D = 0.6 Σ = 5.38 × 10−13 m2 Ca ∈ [0.001; 5] λ ∈ [1; 100] θ0 = 45◦

Table 3.1: Models parameters: Particle and flow. . xˆ2

x2

µext x1 a2

µint

θinit a1

xˆ1

Figure 3.2: Schematic of an ellipsoidal membrane suspended in a fluid of viscosity µext , with aa12 = 2 as used for the following behavior study. These parameters will be used for the Keller & Skalak model as well as the Abkarian, Faivre & Viallat model and the resulting behavior will be displayed in the next section. In all cases, the models are seen to predict the expected equilibrium state at θ = 0◦ . The aim of this study is to understand the trends of both models in an elongational flow for varying capillary number and viscosity ratios. In all the predictions, whatever the model used, the ellipsoid angle decreases from its initial value and tends to 0. The different models and the value of the parameters changes the dynamics of the inclination angle as detailed in the next section. Given the typical behavior of the angle over time (Eq. 3.71), results will be presented in terms of tan(θ).

3.2.2

Effect of the capillary number

In this section, the viscosity ratio λ is kept constant at λ = 1 for all cases. The capillary number is then varied between 0.001 and 0.1. Figure 3.3 shows that the Abkarian, Faivre & Viallat model converges to the solid regime of Jeffery’s theory for decreasing capillary number. In the case of low capillary number, the surface shear modulus Es value is high. The elastic contribution to the Abkarian, Faivre & Viallat model, the

3.2. MODELS BEHAVIOR IN THE CASE OF A PLANAR ELONGATIONAL FLOW

81

second term of Eq. (3.73), dominates the term in sin 2θ. For low capillary number limit (denoted by lc), ω˙lc has the following limit: ω˙ lc = Alc sin 2ωlc , with, Alc =

Es Σ f1 , 2(f2 − λf1 ) µV

(3.80)

(3.81)

which can be solved analytically, as previously presented for θ: ωlc = arctan(tan(ω0 )e2Alc t ).

(3.82)

Since Alc is negative, ω which corresponds to the phase angle of membrane elements will converge towards 0. The dependence of Alc on the surface shear modulus Es shows that for decreasing capillary number, ω will vanish faster and faster. As ω = 0 is an equilibrium state of Eq. (3.80), ω, ˙ the tank-treading rate, also drops to 0. Thus, the initial dynamics of the particle will be described by the complete AFV model but as soon as ω˙ = 0, the resulting dynamics is perfectly described by Jeffery’s theory since Eq. (3.72) becomes Eq. (3.77) in this case. The time for which the dynamics of the particle are described by the complete AFV model, before matching the rigid case, can be estimated using the associated characteristic time scale: τlc ∼

1 2µV (f2 − λf1 ) = , Alc f1 Es Σ

(3.83)

which interestingly does not depend on the extensional rate κ and only compares viscous and elastic effects. This time is a viscoelastic time characterizing the delay between the initial situation and the time where the rigid behavior is obtained.

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 82 CELL IN ELONGATIONAL FLOW 1 K&S AFV Ca = 0.001 AFV Ca = 0.01 AFV Ca = 0.05 AFV Ca = 0.1 AFV Ca = 1.0 Jeffery

tan YY θ

0,1

0,01

0,001

0,0001

0

1

2

3

4 XX κt

5

6

7

8

Figure 3.3: Comparison of the tangent of the orientation of the Keller & Skalak model (K&S) and the Abkarian, Faivre & Viallat model (AFV) for different capillary numbers, for λ = 1 versus the non-dimensional time. The Keller & Skalak model is plotted only for λ = 1 as it does not depend on the capillary number. In the high capillary number limit, the elastic power stored by the membrane in the AFV model is small, especially at the beginning of the simulation. Thus, the results of the AFV model matches the ones from the K&S model when Ca ≫ 1.

3.2. MODELS BEHAVIOR IN THE CASE OF A PLANAR ELONGATIONAL FLOW

83

K&S AFV Ca = 0.001 AFV Ca = 0.01 AFV Ca = 0.05 AFV Ca = 0.1 AFV Ca = 1.0

0,6 0,5

ω˙ YY κ

0,4 0,3 0,2 0,1 0 -0,1

0

1

2

3

4 XX κt

5

6

7

8

Figure 3.4: Comparison of the tank-treading frequency ω˙ of the Keller & Skalak model (K&S) and the Abkarian, Faivre & Viallat model (AFV) for different capillary numbers, for λ = 1 versus the non-dimensional time.

Fig. 3.4 shows the instantaneous tank-treading frequency for the Keller & Skalak and the Abkarian, Faivre & Viallat models. In the case of the Keller & Skalak model, the value of ω˙ will not be affected by the capillary number. It is interesting to note that the tank-treading frequency is seen to be only positive for the K&S model, and vanishes for longer times. This implies that the tank-treading motion will be faster at the beginning of the reorientation of the cell and will slow down rapidly during the motion, keeping the same direction at all times until a full stop when the orientation is completed. For the Abkarian, Faivre & Viallat model, ω˙ does depend on the capillary number. Since the factor in front of Eq. (3.73) is negative, higher values of Es (lower capillary numbers) will result in a quicker decrease of the value of ω˙ over time. This conveniently confirms the behavior witnessed in Fig. 3.3 where the AFV model converges to Jeffery’s theory, as if the value of ω˙ is 0, the resulting model is found to be the same as Jeffery’s (Eq. (3.77)). The value of ω˙ for κt = 0 is the same for all cases, it does not depend on the capillary number. The value of ω˙ converges towards 0 for the AFV model because the elastic energy stored by the membrane is released to the flow during the reorientation, which is why negative values of ω˙ are present. The reorientation process is slowed down by negative values of ω. ˙ The integration of ω˙ shows the elastic energy stored (for positive values of ω) ˙ and released (for negative values of ω). ˙ These two contributions should be equal. The less elastic energy is stored in the membrane, the less time it takes to release it with lower intensity.

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 84 CELL IN ELONGATIONAL FLOW

3.2.3

Effect of the viscosity ratio

In this section, the viscosity ratio λ is varied in the range presented in Table 3.1 in the Keller & Skalak and Abkarian, Faivre & Viallat models while the capillary number is kept constant at Ca = 5. Figure 3.5 shows the tangent of the orientation of the particle with non-dimensional time. It is first seen that with increasing viscosity ratio, the models will converge towards Jeffery’s theory, as expected. The second obvious observation is that the K&S model and the AFV model agree at small times, before the elasticity gets important in the AFV model. 1

0,1

tan YY θ

0,01

K&S L λ ==1 1 λ ==3 3 K&S L λ ==1010 K&S L λ ==100 100 K&S L λ==1 1 AFV L λ==3 3 AFV L AFV L λ==1010 AFV L aa| λ==100100 Jeffery

0,001

0,0001

1e-05

0

1

2

3

4 XX κt

5

6

7

8

Figure 3.5: Comparison of the tangent of the orientation of the Keller & Skalak model (K&S) and the Abkarian, Faivre & Viallat model (AFV) for different viscosity ratios, for Ca = 5 versus the non-dimensional time. Both models depend on λ, as seen in Eq. (3.68) and (3.73). These equations show that for increasing λ, the resulting equations for θ˙ of both models will become the same, namely Eq. (3.77) of Jeffery’s theory. A striking results is the effect of λ on the speed of the reorientation. Due to increasing dissipation with λ, the particle reorients less rapidly when λ increases. The dependence on λ is less significant on the results of the AFV model for lower capillary numbers, as the elastic contribution dominates the ω and ω˙ terms. Figure 3.6 shows the value of the tank-treading frequency ω˙ for varying viscosity ratio λ. It is observed that once again, the Keller & Skalak model only shows positive values of ω˙ as opposed to the positive and negative ones of the Abkarian, Faivre & Viallat model. This corresponds to the membrane rotating only in one direction for the K&S model, and changing direction over time for the Abkarian, Faivre & Viallat

3.2. MODELS BEHAVIOR IN THE CASE OF A PLANAR ELONGATIONAL FLOW

85

model. It is seen that the greater the viscosity ratio λ, the quicker the vanishing of ω. ˙ As explained previously, if the tank-treading frequency ω˙ is 0, the regime of Jeffery’s theory is retrieved, since the contribution of the velocity of the membrane in the K&S and AFV models is negligible in front of viscous effects. This behavior for ω˙ is easily found analytically in Eq. (3.68) and (3.73), where it is seen to vanish in the case of the AFV model and resulting in Jeffery’s theory in both cases. It is also seen that ω˙ for κt = 0 strongly depends on the value of λ. For lower values of λ, the membrane is seen to reorient and spin rapidly. In the AFV model, an important amount of elastic energy is thus stored which strongly slows down θ˙ upon releasing. This is why the curve for λ = 1 in Fig. 3.5 rapidly shows an almost flat trend.

ω˙ κ

λ=1 λ=3 λ = 10 λ = 100 λ=1 λ=3 λ = 10 λ = 100

κt Figure 3.6: Comparison of the tank-treading frequency ω˙ of the Keller & Skalak model (K&S) and the Abkarian, Faivre & Viallat model (AFV) for different viscosity ratios, for Ca = 5 versus the non-dimensional time. The derivation of the Keller & Skalak model and the addition of Abkarian, Faivre & Viallat in the context of a 2D and 3D hyperbolic flow shows that the resulting model is different both from the one derived by Keller et al. [104] and from the one in Abkarian et al. [3]. The case of the elongational flow is simpler, the particles reorienting systematically in the direction of the elongation. The models show the expected behavior, correctly predicting the equilibrium state at θ = 0◦ (the particle is aligned with the direction of extension). It is seen that they also converge towards Jeffery’s theory for high values of the viscosity ratio λ. For increasing capillary number, the AFV model is seen to converge towards Jeffery’s theory again, where the elastic contribution of the membrane will dominate the membrane velocity contribution. In conclusion, the AFV

CHAPTER 3. ANALYTICAL MODELING OF THE DYNAMICS OF RED BLOOD 86 CELL IN ELONGATIONAL FLOW model is found to converge towards Jeffery’s theory for vanishing capillary number or high viscosity ratio. The K&S model is found to have the same behavior in the case of high viscosity ratio. The AFV model is found to converge towards the K&S model for high values of the capillary number. Also, the K&S model is found to reorient faster than the AFV model, as the membrane stores energy as it rotates and releases it by slowing down the reorientation. Finally, the viscosity ratio is seen to have an important influence on the results, θ˙ will be lower for higher values of the viscosity ratio, which corresponds to a slower reorientation. In the next section, these analytical models will be used and compared against the orientation of a red blood cell during its passage inside the micro-orifice.

Chapter

4

Red blood cells dynamics Chapter contents

4.1

4.2

Numerical configuration . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

4.1.1

Red blood cell . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88

4.1.2

Configuration

89

4.1.3

Reduced configuration

. . . . . . . . . . . . . . . . . . . . . . .

97

4.1.4

Initial state of the red blood cell . . . . . . . . . . . . . . . . . .

99

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Results on the centerline . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.2.1

Performed cases . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

4.2.2

General dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 102

4.2.3

Deformation and elongation . . . . . . . . . . . . . . . . . . . . . 106

4.2.4

Comparison with analytical models

4.2.5

Influence of the viscosity ratio . . . . . . . . . . . . . . . . . . . . 113

. . . . . . . . . . . . . . . . 109

This section investigates the dynamics of red blood cells within a blood analyzer geometry using numerical simulations. Results about the deformation, trajectory and orientation of cells while they travel through the device are shown and compared to some of the analytical models presented in Chapter 3. First, the numerical configurations are presented: the red blood cell geometry and the fluid domain, the blood analyzer geometry. Then, the simulation setup along with the associated results about the orientation of centered red blood cells are detailed. Finally, the behavior of the cells is closely analyzed. 87

88

CHAPTER 4. RED BLOOD CELLS DYNAMICS

4.1 4.1.1

Numerical configuration Red blood cell

A single red blood cell is used for all simulations; its shape is a discocyte as presented on Fig. 4.1.

Figure 4.1: Shape of the red blood cell used for the simulations. Whole red blood cell (left), half red blood cell (right). This red blood cell is defined by various parameters as presented in Table 4.1. These parameters are close to physiological parameters experimentally measured and remain identical for all cases performed. Cell and membrane geometry Red cell area Red cell volume Length scale Mechanical parameters

ARBC = 133.4 µm2 VRBC = 93.5 µm3 RRBC = 2.81 µm

Model mechanical law Shear modulus Area dilation modulus Bending modulus Others

Skalak Es = 5.7 × 10−6 N.m−1 Ea = 5.7 × 10−1 N.m−1 Eb = 1.7 × 10−19 N.m

Viscosity ratio Interior viscosity Exterior viscosity Membrane viscosity Conductivity Membrane surface elements

λ=6 νint = 6 × 10−6 m2 s−1 νint = 10−6 m2 s−1 νmemb = 0 σm = 2.27 × 10−12 S.m−1 Nmemb = 5120 faces

Table 4.1: Summary of the chosen red blood cell parameters . RRBC is the characteristic length scale of the red blood cell (radius of the sphere of

89

4.1. NUMERICAL CONFIGURATION same volume) and is defined as RRBC =



3VRBC 4π

1/3

.

(4.1)

It is of interest to note that this discocyte shape for the red blood cell is the result of the deflation simulation presented in Section 2.8.2, which means that it is a pre-stressed equilibrium state of the cell.

4.1.2

Configuration

The simulation of the flow of red blood cells inside the blood analyzer is done using the numerical method described in Chapter 2. The blood analyzer geometry has been provided by Horiba Medical, it is shown in Fig. 4.2. Y

Z

X

Figure 4.2: Picture of the blood analyzer whole geometry. (A) is the sample injector, (B) and (C) are inlets for the sheath fluid and (D) is the outlet of the system. The micro-orifice where the counting/sizing takes place, is not visible in this view. In figure 4.3, a slice of the whole geometry is shown. It is about 2.5 cm long, and 5 mm in diameter. The blood samples are injected with a flowrate of Qsample = 0.51 µL.s−1 in (A) of Fig. 4.3. The sheath flow imposed on (B) is used to hydrofocalize the blood sample, as presented in Section 1.3.3. The sheath flow imposed at (C) is used to evacuate the sample rapidly after the measurement is performed in the microorifice. These sheath flows are injected with flow rates of Qsh1 = 8.38 µL.s−1 and Qsh2 = 75.4 µL.s−1 on (B) and (C) respectively.

90

CHAPTER 4. RED BLOOD CELLS DYNAMICS Anode

Cathode

1 mm

Figure 4.3: Slice of the full geometry. (A) is the sample injector, (B) and (C) are inlets for the sheath fluid and (D) is the outlet of the system. The red circle is centered on the micro-orifice. Bulk velocities are calculated from these flow rates. We find Usample = 1.8 × m.s−1 , Ush1 = 10−2 m.s−1 and Ush2 = 1.5 × 10−1 m.s−1 . The micro-orifice where the measurement is performed is a cylindrical region of diameter 2R = 50 µm and length l = 65 µm. The curvature radius of the enter of the entrance of the orifice is Rc = 15 µm.

10−3

Rc

R Centerline

zf = 2.6R

O

l = 2.6R

Figure 4.4: Slice of the micro-orifice. point O is the entrance of the orifice on the centerline, chosen as the origin of the z coordinate. zf is the exit of the orifice on the centerline. The resulting maximum velocity inside the orifice is about U = 7 m.s−1 . The bulk Reynolds number in the micro-orifice is then: Re =

2U R ≈ 350. νext

(4.2)

4.1. NUMERICAL CONFIGURATION

91

In the same manner, the Reynolds number inside the blood sample injection tube (Res ) and in the region where the blood sample and sheath flow meet (Rem ) can be defined. The resulting Reynolds number are Res ∼ 1 and Rem ∼ 5. The Reynolds number near the sheath fluid injection (B) is ReB ∼ 10. The geometry is discretized using an unstructured mesh of approximately 14 million elements, which is more refined around the micro-orifice. Inlets conditions are imposed on (A), (B) and (C), and an outlet condition is imposed on (D). The physical time simulated in the full geometry presented here is of a few microseconds. These values of the Reynolds numbers in the sample injection tube and the meeting region cause the flow to establish rapidly; even the sheath flow is seen to quickly converge to its axi-symmetric profile with a local injection in (B) (not shown). The value of the Reynolds number in the micro-orifice explains the formation of a jet flow downstream of the micro-orifice, as observed in Fig. 4.5.

Figure 4.5: Fluid flow close to the blood analyzer micro-orifice. Figure 4.6 shows the pressure field inside the micro-orifice. Since the value of the pressure far from this region is constant, the overall pressure drop is around ∆P = 27 000 Pa in the system.

92

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Figure 4.6: Pressure field (Pa) inside the micro-orifice. Pressure is defned up to an arbitray additive constant. z Figure 4.7 shows the particle acceleration ∂u ∂z along the main axis of the blood analyzer. It illustrates the large acceleration that the red blood cells undergo while traveling through the micro-orifice. The position along the main axis is made nondimensional using the radius of the micro-orifice, z = 0 (point O in Fig. 4.4) being the entrance of the orifice. It is seen that the red blood cells are subjected to a strong acceleration and are expected to have an interesting behavior in the neighborhood of the micro-orifice, in terms of deformation (notably elongation) and trajectory.

93

4.1. NUMERICAL CONFIGURATION

∂uz ∂z

(s−1 )

← Orifice →

Longitudinal position (z/R) Figure 4.7: Axial gradient of the axial component of the velocity before and inside the orifice on the centerline.

Figure 4.8 shows the streamlines inside the blood analyzer, when computed from the exit of the sample injection tube. It gives information about the type of trajectories the red blood cells will most likely have, assuming that the red blood cells follow the streamlines in this region. Figure 4.9 shows a close-up of the streamlines inside the micro-orifice and reveals that the sheath flows does focus the red blood cells close to the micro-orifice center axis.

Figure 4.8: Streamlines from the blood sample injection tube exit to the micro-orifice. A local capillary number associated to the red blood cell area dialation modulus can be computed. It is based on the local velocity gradient along the blood analyzer main axis, denoted by k, and a red blood cell area dilation modulus Ea : Caa =

µext kRRBC , Ea

(4.3)

94

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Figure 4.9: Streamlines in the micro-orifice.

with RRBC the characteristic length of the red blood cell, defined as, RRBC =



3VRBC 4π

1/3

.

(4.4)

The field of the capillary number in the micro-orifice is shown on figure 4.10. In the region where red blood cells are present, Caa reaches 0.5 which is a relatively high value, given that it is based on Ea . Thus, small area changes may be present. The value of the capillary number based on the shear modulus would be C = 100000 times larger; significant deformations are thus be expected.

Figure 4.10: Field of local capillary number inside the blood analyzer, in the microorifice area. The extreme streamlines of Fig. 4.8 are reported to visualize the area where red blood cells are most probably found.

4.1. NUMERICAL CONFIGURATION

95

The whole configuration is computationally costly and will not be used for the numerical simulations of the red blood cells in the blood analyzer as it would increase the time of computation dramatically. A reduced configuration with smaller computational domain is thus desirable. The axial gradient of the axial component of the velocity over the centerline of the whole configuration is shown in Fig. 4.11 and 4.12. The former, which presents the axial gradient of the axial component of the velocity on the centerline between the exit of the injection tube and the coordinate z/R = −25, shows that the magnitude of the velocity gradient in the orifice is about 10000 greater than the one experienced by red blood cells after the exit of the injection tube. It is seen that the velocity gradient starts by being negative and increases. It is caused by the fact that the flow of the injection tube meets the sheath flow in this region. Fig. 4.12 shows the axial gradient of the axial component of the velocity on the centerline between the coordinate z/R = −25 and the coordinate of the deposition (z/R = 1) as it will be presented later in this chapter. It is seen that the velocity gradient does not increase significantly for most of the length traveled by the red blood cells. The influence of the velocity gradient could have more important effects from z/R = −10. Even if the dynamics of the red blood cells from z/R = −10 to z/R = −1 would be interesting to investigate, the time of transit in this section can be approximated to be about 3 ms, which is 100 times greater than the time spent in the region considered in this chapter (from z/R = −1 to z/R = 2.6). It was thus decided to reduce the computational domain even further, as detailed in section 4.1.3.

96

CHAPTER 4. RED BLOOD CELLS DYNAMICS 15

∂uz −1 ∂zYY(s )

10

5

0

-5 -125

-100

-75

-50

-25

XX Longitudinal position (z/R)

Figure 4.11: Axial gradient of the axial component of the velocity far from the entrance the orifice on the centerline. z/R ∈ [−125; −25]

50000

∂uz −1 ∂z YY(s )

40000

30000

20000

10000

0 -25

-22,5

-20

-17,5

-15

-12,5

-10

-7,5

-5

-2,5

LongitudinalXX position (z/R) Figure 4.12: Axial gradient of the axial component of the velocity near deposition of the red blood cell, before the orifice on the centerline. z/R ∈ [−25; −1]

4.1. NUMERICAL CONFIGURATION

4.1.3

97

Reduced configuration

As previously stated, a reduced configuration is used in order to reduce computational cost and because the velocity field is virtually homogeneous in most of the blood analyzer (see Fig. 4.3 and above). It is extracted from the whole system. The flow inside the whole blood analyzer is first computed and converged in the absence of cells, as shown in the previous section. The resulting velocity field is then interpolated on a reduced configuration.

Figure 4.13: Crossed-section of the reduced configuration shown over the full configuration. The geometry and flow are axi-symmetric. The indicated boundaries are set as flow inlets. A velocity profile is imposed on these boundaries at all times during the simulations, the boundary on the right is set to be an oulet boundary condition. The velocity field in the reduced configuration is interpolated from the velocity field computed in the full configuration. The mesh contains around 8 million tetrahedral elements, the region where the red blood cells are expected to travel through the system being more refined (see Fig. 4.9). The mesh inside this refined region is fixed to RRBC = 6. The reduced configuration dx and its velocity field, centered on the aperture, are presented on Fig. 4.14.

98

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Figure 4.14: (a) Reduced configuration extracted from the whole blood analyzer. (b) Fluid flow inside the reduced configuration.

uz (m/s)

The velocity field from a reduced domain simulation in the absence of a particle is compared to the one found in the full configuration in Fig 4.15. The maximum difference between the two profiles is found to be of 5%. As we are more interested by the qualitative behavior than by precise numbers, such a difference is considered to be acceptable. This reduced configuration will be used in all the following simulations.

Longitudinal position (z/R) Figure 4.15: Comparison of the z-component of the velocity on the main axis of the full and reduced configurations.

4.1. NUMERICAL CONFIGURATION

4.1.4

99

Initial state of the red blood cell

The initial state of the red blood cells arriving near the orifice is actually unknown. In order to study the dynamics of the red blood cells in the orifice, a parametric study is performed, as a function of the initial position and orientation of the cell. The cell will be supposed to be undeformed where it is deposited in the flow. The initial position and orientation of the red blood cell in the numerical domain is usually defined by six variables: three coordinates (x0 , y0 , z0 ) for the initial position and three angles for the orientation, (θ0 , ϕ0 , ψ0 ). As the initial shape of the red blood cell is axisymmetric (see Fig. 4.1) only two orientation angles are necessary, namely θ0 and ϕ0 . As the geometry and the flow are q axisymmetric, variables x0 and y0 can be gathered using the distance to the axis r = x20 + y02 . Without loss of generality, the red blood cells will be deposited in the y = 0 plane. The initial position of the red blood cell in the numerical domain is first defined by its position z0 along the blood analyzer axis (z-axis), which is fixed for all the simulations presented: z0 = −R, (4.5) with R = 25 µm, the radius of the micro-orifice. This position was chosen so that most of the acceleration phase is computed at a relatively reduced computational cost. This initial coordinate being fixed, three independent parameters are still necessary to describe the red blood cell position and orientation: • r, the initial distance from the z-axis, taken along the x-axis. • θ0 and ϕ0 , two angles used to define the red blood cell orientation. Due to the hydrofocalization, the red blood cells travel will be restricted to the domain bounded by the streamlines, as seen in Fig. 4.9. The position r is found to take all the possible values inside the domain bounded by the streamlines in Fig. 4.16, since the other half would generate redundant results. The resulting range for r is: • r ∈ [0; 15] µm. In order to define the angles θ and φ, a local frame with respect to the streamlines has to be introduced. These axes are (see Fig 4.18): • ~zstream , which is aligned with the streamline for the current position r of the cell. • ~ystream = ~y , which is defined as the vector normal to the plane presented in Fig. 4.17. It matches the vector of the canonical basis. • ~xstream , the axis completing the direct orthonormal basis formed with ~ystream .

100

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Figure 4.16: Reduced configuration around the aperture with streamlines and position marks.

~ystream ~zstream ~ystream ~xstream

~zstream

~xstream

Figure 4.17: Definition of axes ~xstream , ~ystream , ~zstream in the reduced configuration. A local frame for the red blood cell is also introduced. Figure 4.18a shows the ~xRBC and ~zRBC axes which are centered on the center of mass of the red blood cell and aligned with the long and small radii of the cell, respectively. These two axes are chosen to be in the (~xstream , ~zstream ) plane, so that ~yRBC = ~ystream and completes the direct (~xRBC , ~yRBC , ~zRBC ) orthonormal basis.

4.1. NUMERICAL CONFIGURATION

101

Figure 4.18: Representation of the red blood cell local frame and definition of θ. (a) RBC local frame, (b) θ > 0, (c) θ < 0.

As presented on Fig. 4.18b and 4.18c, the angle θ is defined as, cos θ = ~zRBC . ~zstream .

(4.6)

The second angle, ϕ, is the amount of rotation around ~zstream . To define this angle, it proves useful to introduce the projection ~xp of ~xRBC on the plane formed by ~xstream and ~ystream = ~y , which writes, ~xp = (~xRBC . ~xstream )~xstream + (~xRBC . ~ystream )~ystream .

(4.7)

Thus, as seen on Fig. 4.19, the angle ϕ is defined as, cos ϕ = ~xp . ~xstream .

(4.8)

Figure 4.19: Representation of the red blood cell local frame and definition of θ and ϕ.

Both θ0 and ϕ0 should vary inside their own ranges defined by considering the geometrical symmetries of the blood analyzer itself, as well as the ones of the red blood cell. The resulting ranges are: • θ0 ∈ [−90◦ ; 90◦ ]. • ϕ0 ∈ [0◦ ; 90◦ ].

102

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Note that when the red blood cell is located on the z-axis (r = 0 µm), φ0 is no longer useful and the range of θ0 can be reduced to [0◦ ; 90◦ ] for symmetry reasons. This particular configuration will be studied in details in the remaining of this chapter. The cases where r 6= 0 will be considered in Chapter 5.

4.2 4.2.1

Results on the centerline Performed cases

As a first step, the range of explored initial states is restricted in this chapter, decreasing the number of needed variables to uniquely define a single case. The initial position of the red blood cell is fixed on the centerline (r = 0 µm). The performed simulations consider a centered red blood cell, with the longitudinal position still imposed as z0 = −R: (x0 = 0, y0 = 0, z0 = −R), (4.9) In this case (r = 0 µm), ~zstream is aligned with ~z, the main axis of the blood analyzer. The initial orientation of the red blood cell is thus defined only by θ0 , as ϕ0 only describes redundant cases due to the cylindrical symmetry of the geometry. It is thus fixed as ϕ0 = 0◦ so that ~yRBC = ~ystream for any value of θ0 . The angle θ0 will thus vary inside a range defined by considering the geometrical symmetries of the blood analyzer and the red blood cell. The resulting range is: • θ0 ∈ [0◦ ; 90◦ ]. In this section, several cases are presented. The value of θ0 will change from θ0 = 0◦ to θ0 = 90◦ while the initial position (x0 = 0, y0 = 0, z0 = −R) remains unchanged.

4.2.2

General dynamics

The focus is set on the time evolutions of the longitudinal position and the orientation of the cells. General dynamics and behaviors are extracted from the results in an attempt to better understand the dynamics of the cells. Time evolution of the longitudinal position of the cells The results are extracted from two simulations where the initial orientation of the cell θ0 is set at θ0 = 0◦ (long axis of the cell perpendicular to the direction of the flow) and θ0 = 90◦ (long axis aligned with the direction of the flow). The results are shown in Fig. 4.20. The time evolution of the longitudinal position of the cell is seen to be almost perfectly identical for the two cells of very different orientation. More interestingly, the time evolution of the longitudinal position of a tracer convected by the velocity field

103

Longitudinal position (z/R)

4.2. RESULTS ON THE CENTERLINE

θ0 = 0◦ θ0 = 90◦

Time (µs) Figure 4.20: Time evolution of the longitudinal position (z/R) for red blood cell with an initial orientation of θ0 = 0◦ and θ0 = 90◦ . The time evolution of the longitudinal of a tracer deposited in the velocity field is also plotted. The dashed lines represent the coordinate of the exit of the orifice and the corresponding time at which the red blood cell center of mass reaches this coordinate.

is reported on Fig. 4.20 and also shows to be identical. This shows that the initial orientation of the cell, the deformation and rotations it undergoes or the presence of the cell itself during the travel in the orifice has little to no effect on the trajectory of its center of mass. In addition, it allows the measurement of the average time of transit in the micro-orifice τt which is found to be: τt = 25.9 µs. Time evolution of the orientation of the cells The deformation states of the red blood cells along their pathlines are presented on Fig. 4.21 for four different initial orientations (θ0 = 0◦ , 30◦ , 60◦ , 90◦ ). Figure 4.21b and 4.21c clearly show that the velocity field in the micro-orifice is reorienting the red blood cell towards the inclination θ0 = 90◦ (see Fig. 4.21d). The first case, with θ0 = 0◦ is not reoriented because of the symmetry of the flow and of the red blood cell. In such a case, the cell slightly elongates in the streamwise direction and contracts in the other directions due to the action of the flow. Note that such a symmetric case can be viewed as a pathological case, probably never observed in reality. The simulations have been repeated for more values of θ0 . The orientation of the cell with respect to the longitudinal axis is plotted on Fig. 4.22. This figure shows

104

CHAPTER 4. RED BLOOD CELLS DYNAMICS

a)

b)

c)

d)

Figure 4.21: Successive positions of the red blood cell with initial orientations: a) θ0 = 0◦ , b) θ0 = 30◦ , c) θ0 = 60◦ , d) θ0 = 90◦ .

that for all values of θ0 , a reorientation of the cell is observed. Depending on its initial inclination, the final orientation of the cell will be more or less close to 90◦ . The performed cases have been labeled as RX on Fig. 4.22 and 4.23, with X the value for θ0 . This reorientation phenomenon is the result of the interaction between the fluid (hydrofocalization and resulting hydrodynamic forces) and the cell mechanics: membrane resistance to deformation and viscosity difference. It is investigated in details in Section 4.2.4 using some of analytical models developed in Chapter 3.

105

Orientation (◦ )

4.2. RESULTS ON THE CENTERLINE

Longitudinal position (z/R)

Orientation (◦ )

Figure 4.22: Instantaneous orientation versus position along the axis of the blood analyzer, for different initial orientations. The cases are labeled as RX, with X the value for θ0 .

Time (µs) Figure 4.23: Instantaneous orientation versus time, for different initial orientations. The cases are labeled as RX, with X the value for θ0 .

106

4.2.3

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Deformation and elongation

The deformation and elongation state of the red blood cell for θ0 = 90◦ is investigated in this section. In this case, the axis of the red blood cell coincide with the canonical axes of the geometry so that: ~xRBC = ~zstream , ~zRBC = ~xstream and ~yRBC = ~ystream . First, the classical Taylor deformation parameter used in the linear shear flow test case of Section 2.8.4 is introduced again: D=

A−B , A+B

(4.10)

where A and B are the semi-axes of the ellipsoid that has the same tensor of the moments of inertia as the red blood cell. In this chapter, three Taylor parameters are defined, each of them corresponding to the deformation parameter in one plane. The first one, Dxy , using the semi-axes of the ellipsoid in the x − y plane, which normal is the blood analyzer main axis, the second one, Dxz , using the semi-axes of the ellipsoid in the x − z plane, which is used to present the results throughout this thesis and the third one, Dzy , the plane perpendicular to the latter plane, containing the main axis of the blood analyzer, ~z. These planes are represented in Fig. 4.24, along with a red blood cell orientated at θ0 = 90◦ and the associated deformation parameters. ~zstream

Dxz ~xstream Dxy

~ystream

Dzy

Figure 4.24: Representation of the planes associated with Dxy , Dzy and Dxz , with a red blood cell with θ0 = 90◦ . The results for the evolution of Dxy , Dzy and Dxz along the trajectory of a red blood cell with θ0 = 90◦ are presented in Fig. 4.25.

107

4.2. RESULTS ON THE CENTERLINE

Taylor parameters

Dxy Dzy Dxz

Longitudinal position (z/R) Figure 4.25: Taylor parameters versus longitudinal position in the micro-orifice, for θ0 = 90◦ .

Fig. 4.25 shows that the red blood cell deforms in all planes. It can be related to Fig 4.21d, where the red blood cells shows an elongation along the z-axis while undergoing compression in all directions in the x − y plane, correctly corresponding to an increase in Dxz and Dyz and a decrease in Dxy . Dzy shows a initial value close to 0 because of the red blood cell symmetry. All of the deformation parameters show a slow decrease to an asymptotic value after the first part of the orifice (z/R > 1.5), which can be seen as a relaxation of the red blood cell to an asymptotic deformation state, after being subjected to the velocity gradient at the enter of the orifice (see Fig. 4.15). This suggests that when a red blood cell is deposited upstream from the current location (z0 = −R), it could converge to this asymptotic deformation state as well. Dzy and Dxz reach close values in the perforation while Dxy gets small: the red blood cell reaches an elongated quasi-axisymmetric shape. In Fig. 4.26, the Taylor parameters of a red blood cell with an initial orientation of θ0 = 90◦ deposited at z/R = −1.8 is plotted. It shows that the Taylor parameters keep the same trend as in the case of a later deposition as in Fig. 4.25. Some differences in the values of the parameters are to be noted: Dxz reaches a lower value (the same as Dzy ) than in the case of a later deposition. This indicates that the red blood cells reaches a shape that is almost perfectly axi-symmetric, even more than in the case of a late deposition, probably because the red blood cell is subjected to a velocity gradient in the direction of the flow for a longer time. Dxy reaches a smaller value than in the earlier case, and is seen to hit the value of Dxy = 0, indicating that the cell is of circular

108

CHAPTER 4. RED BLOOD CELLS DYNAMICS

shape when projected in the xy-plane. In this case again, the Taylor parameters reach an asymptotic value after they crossed the entrance of the orifice.

Taylor parameters YY

0,6

D D11 xy D22 D zy D D33 xz

0,5

0,4

0,3

0,2

0,1

0 -2

-1,5

-1

-0,5

0

0,5

1

1,5

2

2,5

LongitudinalXXposition (z/R) Figure 4.26: Taylor parameters versus longitudinal position in the micro-orifice, for θ0 = 90◦ and a deposition coordinate of z/R = −1.8.

Elongation

lx /lx0 ly /ly0 lz /lz0

Longitudinal position (z/R) Figure 4.27: Elongation versus longitudinal position in the micro-orifice, for θ0 = 90◦ . Figure 4.27 shows the elongation in each direction. It is easy to see that the elon-

109

4.2. RESULTS ON THE CENTERLINE

gation along the z-axis is correctly observed here, with a maximum elongation of less than 2. Along the y-axis (perpendicular to the plane used in all figures), a decrease of the elongation is clearly seen, with this axis shrinking to about 40% of its initial length. Along the x-axis, no strong variation is observed. This can be justified by the fact that along this axis, the red blood cell already has its smaller length (the dimples of the red blood cell), and due to the incompressibility of the cell and its membrane mechanics, most of the compression will act along the y-axis. 2

1,75

llzz x /lx0 llyy y /ly0 lz /lz0 lxxaaaa

Elongation YY

1,5

1,25

1

0,75

0,5

0,25 -2

-1,5

-1

-0,5

0

0,5

1

1,5

2

2,5

LongitudinalXXposition (z/R) Figure 4.28: Elongation versus longitudinal position in the micro-orifice, for θ0 = 90◦ . Again, a red blood cell with an initial orientation of θ0 = 90◦ is deposited at z/R = −1.8. The results in terms of elongation are plotted in Fig. 4.28. The trend of the elongation is again conserved, even though the dynamics are different, with more important variations in the elongation over time. Nevertheless, the value of the elongation at the end of the orifice is found to be identical with a red blood cell deposited at z/R = −1. This supports the fact that the red blood cells reaches an asymptotic deformation state in the orifice which does not depend on the initial position of the cell for z0 ≪ −1, thus justifying the size of the reduced computational domain used in this study (Fig. 4.13).

4.2.4

Comparison with analytical models

In this section, the analytical models presented in Chapter 3 will be used to model the orientation of the red blood cell along its travel through the micro-orifice. The three models will be used: Jeffery’s theory [99], the Keller & Skalak model [104] and the

110

CHAPTER 4. RED BLOOD CELLS DYNAMICS

Abkarian, Faivre & Viallat model [3] in the context of a 3D elongational flow. Each of these models predicts the orientation of fixed-shape ellipsoids with different membrane properties. They take as input parameters the shape of the particle and the velocity gradient of the flow. Along their trajectories, red blood cells are subjected to a varying elongation rate. In order to reproduce this effect, the trajectory of the particles is first calculated from the carrying flow and the local elongation rate is obtained. The extensional rate κ is thus updated depending on the position of the cell. The value for a1 , a2 and a3 that correspond to the semi-axes of the ellipsoid in the analytical models are set according to the initial value of the semi-axes of the ellipsoid of inertia that has the same moments of inertia as the red blood cell. All the other parameters for the cell (volume, membrane surface, etc.) are set equal to the simulations parameters. Note that some attempts were made to update the shape parameters with simulation data, without convincing improvements. The shape is thus fixed in the models. The goal of this section is to be able to model the trend of the reorientation of the red blood cells in the blood analyzer, and to determine if the previous models are accurate enough to predict the reorientation phenomenon in detail. Note that the computations do no verify the hypothesis of zero Reynolds number of all analytical models, as the maximum particle Reynolds number is of the order of 1.0. Results The three models have been used to compute θ(z) for the cases with θ0 = 10◦ , θ0 = 40◦ and θ0 = 70◦ . Results are presented in Fig. 4.29, 4.30 and 4.31, respectively. It is seen that the general trend of the reorientation is predicted in all cases. However, the details of the phenomenon are not described: notably, for θ0 = 10◦ , the models do not capture the stabilization of the orientation seen in the simulation around z/R = 1.25−1.5, where the models keep producing an increasing orientation. This stabilization is most likely due to mechanical effects related to the shape deformation which are not accounted for in the analytical models.

111

4.2. RESULTS ON THE CENTERLINE 70

Orientation (◦ ) YY

60

K&S AFV Jeffery YALES2BIO

50

40

30

20

10 -1

-0,5

0

1

0,5

1,5

2

2,5

XX

Longitudinal position (z/R) Figure 4.29: Comparison of the instantaneous orientation of the cell from the YALES2BIO simulation with the analytical models for θ0 = 10◦ .

Figures 4.29, 4.30 and 4.31 show that the initial slope of the orientation is closer to the numerical results as the initial orientation θ0 increases. This can be justified by the fact that for higher values of θ0 , the deformation is less significant at the beginning of the simulation and thus, does not play an important role in the reorientation phenomenon at small time scales. It is also of interest to note that Jeffery’s theory is always found to underestimate the reorientation more than the Keller & Skalak and Abkarian, Faivre & Viallat models which is consistent with the parametric comparison of the models in Chapter 3. In addition, the Keller & Skalak and Abkarian, Faivre & Viallat models show an almost exactly similar behavior, which leads to think that the elastic contribution of the membrane is not important for the reorientation in this case. This can be understood by the fact that the red blood cell is subjected to a flow with high capillary numbers, thus minimizing the elastic contribution to the model. The differences in reorientation between the Keller & Skalak or Abkarian, Faivre & Viallat models and Jeffery’s theory are thus related to the consideration of the internal fluid.

112

CHAPTER 4. RED BLOOD CELLS DYNAMICS 90

Orientation (◦ ) YY

80

70

60

K&S AFV Jeffery YALES2BIO

50

40 -1

-0,5

0

1

0,5

2

1,5

2,5

XX

Longitudinal position (z/R) Figure 4.30: Comparison of the instantaneous orientation of the cell from the YALES2BIO simulation with the analytical models for θ0 = 40◦ .

90

Orientation (◦ ) YY

85

80

K&S AFV Jeffery YALES2BIO

75

70 -1

-0,5

0

1

0,5

1,5

2

2,5

XX

Longitudinal position (z/R) Figure 4.31: Comparison of the instantaneous orientation of the cell from the YALES2BIO simulation with the analytical models for θ0 = 70◦ . These results suggest that the reorientation phenomenon depends on the viscosity ratio between the inner and outer fluids, since the K&S and AFV models are found

113

4.2. RESULTS ON THE CENTERLINE

to predict more correctly the orientation. Such a statement will be examined in more details in the next section. However, the elastic energy storage of the AFV model is not found to play a role in the reorientation. In addition, some physical effect is obviously absent from the analytical models, as all of them miss some features of the reorientation. This effect is most likely the deformation of the cell, that will have a direct impact on its orientation. Thus, unless an analytical model with deformable particles is developed, performing numerical simulations of the kind presented here is crucial to get details on the orientation of red blood cells in such a geometry and a flow.

4.2.5

Influence of the viscosity ratio

As presented in Table 4.1, the viscosity ratio of a healthy red blood cell is around 6, which is the value that has been used in all simulations so far. In this section, the influence of the viscosity ratio is studied in more details, independently of its industrial relevance in this particular problem. The aim is to further understand the dynamics of red blood cells in elongational flows. Jeffery’s theory is thus not considered anymore, as it does not account for on the value of the viscosity ratio. a)

b)

c)

d)

Figure 4.32: Successive deformation states of the cell for different viscosity ratios, for θ0 = 10◦ . a) λ = 1, b) λ = 2.5, c) λ = 6 and d) λ = 10. Additional simulations have been performed with different values for the viscosity ratio. They show great impact of the value of this ratio on the dynamics of the orientation of the red blood cell as well as on its deformations, as seen in Fig. 4.32 and 4.33. It is seen that the higher the viscosity ratio, the smaller the reorientation. This effect is more important for lower angles. This is not surprising since lower angles cases, red blood cells will undergo an important compression in the x − y plane, in all directions (perpendicular to the blood analyzer axis), resulting in greater deformation and thus

114

CHAPTER 4. RED BLOOD CELLS DYNAMICS

contributing more importantly to the reorientation phenomenon as the orientation is calculated from the tensor of the moments of inertia. Thus, if the viscosity ratio increases, the deformation and reorientation of the red blood cells will be less pronounced. The 3D Abkarian, Faivre & Viallat extension of Keller & Skalak’s model derived for a hyperbolic flow in Chapter 3 is now used to model the orientation of the red blood cell during its passage in the micro-orifice. The original Keller & Skalak model is not used as it is close to the Abkarian, Faivre & Viallat model for all values of the viscosity ratio λ considered. The results are presented for 3 values of the viscosity ratio: λ = 1, 2.5, 6 and 10 with two initial orientations: θ0 = 10◦ , 40◦ in Fig. 4.33, 4.34 and 4.35, respectively. 90 80

Orientation (◦ ) YY

70 60 50 YALES2BIO Lλ== 11 YALES2BIO Lλ==2.52.5 YALES2BIO Lλ== 66 YALES2BIO Lλ== 1010 AA AFV Lλ== 11 AFV Lλ==2.52.5 AFV Lλ== 66 AFV Lλ== 1010

40 30 20 10 -1

-0,5

0

0,5

1

1,5

2

2,5

XXposition (z/R) Longitudinal

Figure 4.33: Instantaneous orientation of the cell along its trajectory for different viscosity ratios versus the longitudinal position in the blood analyzer, for θ0 = 10◦ . Figures 4.33, 4.34 and 4.35 show that the AFV model reproduces the dependence of the results on the viscosity ratio: the higher the viscosity ratio, the higher the reorientation phenomenon. However, the case of θ0 = 10◦ is found not to show great agreement for viscosity ratios different of λ = 6. This is due to the fact that the deformation gets more important for low values of the viscosity ratio.

115

4.2. RESULTS ON THE CENTERLINE 90

◦ Orientation Orientation ( )

80

70

YALES2BIO Lλ== 11 YALES2BIO Lλ==2.52.5 YALES2BIO Lλ== 66 YALES2BIO Lλ== 1010 AA AFV Lλ==1 1 AFV Lλ==2.52.5 AFV Lλ==6 AA 6 AFV Lλ==1010

60

50

40 -1

-0,5

0

1

0,5

1,5

2

2,5

Positionposition (R) Longitudinal (z/R)

Figure 4.34: Instantaneous orientation of the cell along its trajectory for different viscosity ratios versus the longitudinal position in the blood analyzer, for θ0 = 40◦ . 90 88

◦ Orientation Orientation ( )

86 84 82 YALES2BIO Lλ== 1 1 YALES2BIO Lλ== 2.52.5 YALES2BIO Lλ== 6 6 YALES2BIO Lλ== 1010 AA AFV Lλ== 1 1 AFV Lλ== 2.52.5 AFV Lλ== 6 AA 6 AFV Lλ== 1010

80 78 76 2

2,1

2,2

2,3

2,4

2,5

2,6

Positionposition (R) Longitudinal (z/R)

Figure 4.35: Instantaneous orientation of the cell along its trajectory for different viscosity ratios versus the longitudinal position in the blood analyzer, for θ0 = 40◦ for z/R ∈ [2; 2.5]. The deformations undergone by the red blood cell have an impact on the reorientation process itself, but also on the measure of this orientation. As it was previously

116

CHAPTER 4. RED BLOOD CELLS DYNAMICS

z y x

Before z/R = −0.2

z y x

After z/R = 0.1

Figure 4.36: Pictures of the red blood cells extracted from numerical simulations, before and after the instant of the bump on the orientation, for θ0 = 40◦ and λ = 1.

mentioned, the orientation is measured by means of the matrix of inertia and the bumps on Fig 4.34, for λ = 1 at z ≃ 0.15R (and less significantly for λ = 2.5) show the dependence of the orientation on the shape of the cell. Indeed, at the instant of these bumps, the red blood cell dimples are seen to buckle as seen on Fig. 4.36. Aside from this bump effect on orientation, the Abkarian & Viallat variation model is seen to be in good agreement with the results for θ0 = 40◦ with different viscosity ratios. This is due to the fact that for higher values of θ0 , the deformation will most likely have less impact on the orientation and will thus allow the model to capture its trend. This shows that once this model has been adapted to strain flow, it is able to predict the trend of the orientation of red blood cells in this type of flows, as long as deformation does not play an important role in this process.

Chapter

5

Electrical response Chapter contents 5.1

5.2

5.3

5.4

5.5

5.6

5.7

Red blood cells interaction with electrical fields . . . . . . . . 5.1.1 AC case . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 DC case . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical method . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Modelling framework . . . . . . . . . . . . . . . . . . . 5.2.2 Computational treatment . . . . . . . . . . . . . . . . Validation test cases . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Numerical domain . . . . . . . . . . . . . . . . . . . . 5.3.2 Constant conductivity . . . . . . . . . . . . . . . . . . 5.3.3 Variable conductivty . . . . . . . . . . . . . . . . . . . Computation of the electrical potential in the absence of cells 5.4.1 Full configuration . . . . . . . . . . . . . . . . . . . . . 5.4.2 Reduced configuration . . . . . . . . . . . . . . . . . . Electrical pulse generated by the passage of a cell . . . . . . . 5.5.1 Effect of the passage of a cell . . . . . . . . . . . . . . 5.5.2 Procedure of generation of the electrical pulse . . . . . 5.5.3 Range of the study . . . . . . . . . . . . . . . . . . . . 5.5.4 Test case with a rigid particle . . . . . . . . . . . . . . The probability density function of volume measurements . . 5.6.1 Why a statistical approach . . . . . . . . . . . . . . . 5.6.2 Random variables: initial position and orientation . . 5.6.3 Raw results . . . . . . . . . . . . . . . . . . . . . . . 5.6.4 Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6.5 Probability density function . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

117

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

118 118 119 119 119 122 123 124 125 126 128 128 130 133 133 133 134 135 135 135 136 142 145 146 150

118

CHAPTER 5. ELECTRICAL RESPONSE

This chapter presents the modeling of the electrical measurement of the volume of red blood cells passing through a blood analyzer. In a blood analyzer, the Coulter effect allows the measurement of the volume of cells. The insulating nature of the red blood cells membrane creates a perturbation in the resistance of the system that varies as they travel through the device. This perturbation, in the form of an electrical pulse, is measured and its amplitude is directly related to the cell volume. The measure of this amplitude allows the characterization of possible errors on the measurement of the volume of red blood cell. In this chapter, a brief summary of the interaction of red blood cells with electrical fields is first presented. Then, the numerical method is described. Several test cases are used to assess the quality of the implementation of this method, notably for the more extreme case of discontinuous conductivity. Finally, the method is applied to the blood analyzer problem: the numerical framework and setup used to perform the numerical simulation is first described, followed by a thorough description of the post-processing procedure (namely fittings and a statistical approach) allowing the extraction of meaningful data on the quality of the electrical measurement of the volume of red blood cells in the blood analyzer.

5.1

Red blood cells interaction with electrical fields

Red blood cells are complex deformable particles that have been thoroughly described from the mechanical point of view in Chapter 4. The interaction between red blood cells and an electrical field is now considered. It can take different forms depending on the nature of the current.

5.1.1

AC case

AC stands for alternating current, where the flow of electrical charges periodically changes its direction. It is defined by its intensity and its frequency. This case may be related to particle analyzers, discriminators and manipulators, where it is possible to witness various phenomena such as electrorotation, dielectrophoresis and electroporation [8, 9, 44, 45, 47, 197]. This is due to the dielectric nature of red blood cells in AC electrical field, where the cell is almost an insulator while being able to be polarized. Electrorotation Electrorotation is the rotational movement induced by a polarizing electrical field. It may be caused by the interaction between the alternative electrical field phase and the polarization of the cell [101]. It allows the measurement of various properties of the cell, such as the conductivity and permittivity of the ’inside’ or membrane of cells. Dielectrophoresis Dielectrophoresis is a phenomenon that happens when a cell is subjected to a non-uniform electrical field. It will cause a force to be applied to this cell:

5.2. NUMERICAL METHOD

119

it is polarized and the electric field is different on each pole thus creating a movement whose direction and intensity are determined by the cell electrical characteristics [101]. Electroporation Electroporation is an interaction of an AC generated electrical field with a cell membrane, particularly the lipid bi-layer. It consists of pores opening in the lipid bi-layer. It can be used to introduce DNA inside cells for medical purposes, but depending on the intensity of the electrical field applied, it can cause an irreversible electroporation, compromising the cell viability [203].

5.1.2

DC case

DC stands for direct current, where the flow of electrical charges keeps the same direction at all times. It is only defined by its intensity. In this case, red blood cells are considered as pure insulators that cannot be polarized, thus no effect of the electrical field on the cell has to be taken into account. This chapter is focused on the DC case and the red blood cell will be considered as a pure insulating particle, as the DC context is used in the actual configuration of the blood analyzer designed by Horiba Medical.

5.2

Numerical method

The numerical method chosen for the simulation of the electrical response of red blood cells flowing in a blood analyzer is based on the electrostatic approximation. The construction of the model equation is first shown, and its numerical resolution is presented. This part of the work has been done in collaboration with M. Martins Afonso (postdoctoral fellow at IMAG 2011-2013).

5.2.1

Modelling framework

Model equation In order to model the electrical response of red blood cells flowing in a blood analyzer, the electrostatic approximation is made. Indeed, the physical velocities encountered in the system are several orders of magnitude smaller than the speed of light. In this case, the electric field adjusts instantaneously to any change in the spatial configuration. Also, any Doppler effect or two-way coupling between the electric field and the hydrodynamic or magnetic ones are neglected, following [96]. This approximation allows the computation of the dynamics of the cells to be performed separately from the electrostatic computations, since the former is entirely governed by the fluid-structure interaction treated in Chapter 4. There are two relevant electric parameters in this framework: • the static conductivity σ, • the static relative dielectric constant ǫR .

120

CHAPTER 5. ELECTRICAL RESPONSE

These two parameters vary in space and time, taking one of three distinct values depending on the position considered: either inside the cell (cytoplasm), outside of the cell (carrying fluid) or inside the membrane. These values can be smoothed using a numerical interpolation in order to avoid discontinuities. This smoothing operation is presented in the following part of this section. The fourth of Maxwell’s equations (Amp`ere’s law) in a conductive medium writes: ~ = ~j + ∂t D, ~ ∇×H

(5.1)

~ the magnetizing, ~j the free electric current density and D ~ the electric displacewith H ment. The divergence of Eq. (5.1) reads: ~ = 0. ∇ · (~j + ∂t D)

(5.2)

Two constitutive relations can be used in the case of a conductive medium. First, Ohm’s law, which writes: ~ ~j = σ E, (5.3) ~ the electrical field. The second relation holds for linear and isotropic media, with E which is a reasonable assumption since it is the case for both the interior and exterior fluids. See [149] for a discussion about the possible electric anisotropy of the cell membrane. This relation reads: ~ ~ = ǫ0 ǫR E, (5.4) D where ǫ0 ≃ 8.85 × 10−12 F.m−1 is the electric permittivity of vacuum. Using Eq. (5.3) and (5.4) in Eq. (5.2), the following expression is obtained: ~ = 0. ~ + ∂t (ǫ0 ǫR E)] ∇ · [σ E

(5.5)

Note that Eq. (5.5) does not account for the magnetic field anymore. A consequence of the electrostatic approximation is the irrotationality of the electrical field (because of the Maxwell-Faraday equation), leading to the introduction of the electrical potential V as, ~ = −∇V. E (5.6) Therefore, the following equation is easily obtained using Eq. (5.6) and Eq. (5.5): ∇ · [σ∇V + ǫ0 ∂t (ǫR ∇V )] = 0.

(5.7)

A summary of the dielectric parameters for all the different elements of the system is presented in Table 5.1.

121

5.2. NUMERICAL METHOD Dielectric parameters Component External fluid Cytoplasm Membrane

Conductivity (S.m−1 ) 2.27 0.31 1 × 10−6

Permittivity (F.m−1 ) 80 59 4.44

Table 5.1: Dielectric parameters in the system as provided by Horiba Medical and found in Valero et al. [191].

Both liquids can thus be seen as excellent conductors and the membrane as an insulator, as the conductivity and permittivity values in Table 5.1 suggest. The insulating nature of the red cell membrane allows a simplifying assumption that will be used in the electrostatic solver: the electrical field in the system will not interact with the interior fluid (cytoplasm) of the red cell, since it is wrapped in the insulating membrane. Considering the non-dimensional counterpart of Eq. (5.7) with respect to the time (t = τf t∗ ), it can be written: ∇ · [∇V +

ǫ0 ǫR ∂ ∇V ] = 0. τf σext ∂t∗

(5.8)

with τf the characteristic time scale associated to the flow. From Chapter 4, one can extract the mean velocity inside the orifice (Umean = 3.5 m.s−1 ). The time τf needed for a red blood cell to travel for a distance corresponding to its own size is thus of the order of: τf ∼ 10−6 s.

(5.9)

R will In Eq. (5.8), both ∇V and ∂t∂∗ ∇V are of the same magnitude. Thus, the term τfǫ0σǫext dictate which terms dominates the other one. It can be seen as a ratio of characteristic time scales: τelec ǫ0 ǫR , (5.10) = τf σext τf

where τelec is the electrical time scale of the external medium, since it is the only medium considered as previously justified by the insulating nature of the membrane. Using the values provided by Horiba Medical, the external fluid electric time scale is found: τelec = τexternal =

ǫ0 ǫext R ∼ 10−10 s. σext

(5.11)

Eq. (5.11) shows that the electrical time scale is several orders of magnitude lower than the time scale associated to the flow. It means that the electrical effects will happen significantly faster than the displacement of the particle and thus allows the simplification of Eq. (5.7) into: ∇ · [σ∇V ] = 0,

(5.12)

122

CHAPTER 5. ELECTRICAL RESPONSE

which is the equation that will be solved in the electrostatic simulations. The boundary conditions at the interfaces between the different materials impose the continuity of the tangential component of the electric field, resulting in, Etint = Etext ,

(5.13)

and a possible jump condition for the normal component of the electric-displacement field due to the surface charge density Σ [79]: Σ = Dnint − Dnext = ǫRint Enint − ǫRext Enext .

(5.14)

Another important condition is the conservation of charge at the interfaces. In general, this gives rise to a transport equation for Σ, here interpreted as a space-time field defined on the interface (see [185], [192], [122]). However, in the present study, following [196] and [22], the normal component of the electric current is considered continuous: jnint = jnext ⇒ σint Enint = σext Enext .

5.2.2

(5.15)

Computational treatment

In reality, red blood cells show a conductivity field as presented in Fig. 5.1a, where the conductivity inside the membrane is low and the cytoplasm is conducting with σcyto = 0.31 S.m−1 as recapitulated in Table 5.1. In the numerical simulations performed with YALES2BIO, following the previously stated assumption from the insulating nature of the membrane, the conductivity field σint inside the membrane is set to σint ∼ 10−12 S.m−1 (arbitrarily small value to mimic the insulating nature of the membrane) and is smoothed from the exterior value of σext = 2.27 S.m−1 through the membrane, as seen in Fig. 5.1b and presented in the next section. This approximation corresponds to σext Enext ≃ 0 in Eq. (5.15) because the conductivity in the membrane is low (σint ≪ 1 on the membrane).

123

5.3. VALIDATION TEST CASES a)

σ(r)

b)

σ(r)

σext

σext σcyto σint σmemb

0

Ext.

Membrane

r 0

Int.

Ext.

r

Membrane

Int.

Figure 5.1: Schematic representation of the conductivity field σ. a) In an actual red blood cell. b) For the YALES2BIO simulations.

In Fig. 5.2, the conductivity field σ in the blood analyzer is presented. The model equation, ∇.(σ∇V ) = 0, (5.16) is solved everywhere over the volume of the configuration. A discretization method of second order based on the DGA [177] is used, along with the DPCG procedure [128] to solve the linear system. Boundary conditions are imposed on the walls of the configuration: Dirichlet conditions for the electrodes with V = Velectrode , and homogeneous Neumann conditions on the insulating walls with ∇V = ~0.

σint σext

Insulating walls Anode

Cathode

Figure 5.2: Schematic representation of the conductivity field σ in the blood analyzer.

5.3

Validation test cases

This section aims at providing a set of test cases demonstrating the ability of the solver to solve electrostatic problems relevant to the final application, through the comparison of an analytically computed potential with its simulated counterpart using YALES2BIO.

124

CHAPTER 5. ELECTRICAL RESPONSE

The setup uses a unique and simple geometry with appropriate boundaries, in order to avoid numerical difficulties and to be able to extend this validation to higher dimensions. This geometry is a domain enclosed by concentric spheres, with uniform Dirichlet conditions on both the inner and outer boundaries and angle-independent effective conductivity (see Fig. 5.3).

Figure 5.3: Schematic representation of the test case (slice of the 3D domain) Analytically, the problem is one-dimensional, and the solution reads: V (r) = K + C

Z

r(1−d) dr σ(r)

(5.17)

with d = 2 or 3 the dimension, constants K and C depending on the values of the potential on the internal and external boundaries, Vi and Vo , respectively. From the numerical point of view, however, the problem is solved with an unstructured mesh, and is thus multidimensional, as we do not take advantage of the symmetry of the problem to adopt the mesh. For all the following test cases, the geometry, as well as the potential on the internal and external boundaries are fixed as follows: • R = 0.19 m, with ri = 0.01 m and ro = 0.2 m. • Vi = V (ri ) = −0.1 V and Vo = V (ro ) = 1 V. Only the expression for the conductivity σ(r) changes between the different test cases.

5.3.1

Numerical domain

The numerical domain used in all the following test cases is presented on Fig. 5.4

125

5.3. VALIDATION TEST CASES

Figure 5.4: Slice of the 3D numerical domain used in the test cases.

This mesh contains around 1 million tetrahedral elements, with approximately 20 edges around a perimeter of the inner sphere.

5.3.2

Constant conductivity

In this particular test case, the conductivity σ(r) is considered constant: σ(r) = σ = 0.1 S.m−1 .

(5.18)

Using Eq. (5.17) with d = 3, the analytical profile is easily obtained: Vana (r) = Kc − 







Cc , σr

(5.19)

−Vo o and Cc = σ Vrii −r ri ro . The comparison between Vana and with, Kc = Vi + ro Vrii −V −ro o the simulated potential is presented in Fig. 5.5.

126

Potential (V)

CHAPTER 5. ELECTRICAL RESPONSE

Radial position (m) Figure 5.5: Comparison between the analytical profile and the simulated profile of the potential, in the case of a constant conductivity σ.

A perfect agreement between the analytical and simulated potential profiles is obtained.

5.3.3

Variable conductivty

In this test case, the conductivity σ(r) is considered variable, following: σ(r) =

σ0 . r

(5.20)

Again, using Eq. (5.17) with d = 3, the analytical profile is easily obtained: Vana (r) = Kv + Cv ln( 



−1

σ0 ), r

(5.21) 



−1

and Cv = (Vi − Vo ) ln rroi with, Kv = Vi − (Vi − Vo ) ln( σri0 ) ln rroi comparison between Vana and the simulated potential is presented in Fig. 5.6.

. The

127

Potential (V)

5.3. VALIDATION TEST CASES

Radial position (m) Figure 5.6: Comparison between the analytical profile and the simulated profile of the potential, in the case of a variable conductivity σ(r) = σr0 .

It is seen that there is a good agreement between the analytical and simulated potential profiles.

128

5.4

CHAPTER 5. ELECTRICAL RESPONSE

Computation of the electrical potential in the absence of cells

In this section, the numerical configuration and setup are presented in detail. The computation of the electrical potential inside the full and reduced computational domains are computed and described in the absence of cells.

5.4.1

Full configuration

The electrical potential inside the full configuration is computed using the electrostatic solver presented in the previous section. The conductivity field inside the geometry is constant, it corresponds to the conductivity of the external fluid when a cell is considered. σext = 2.27 S.m−1 . (5.22) Mesh and boundary conditions The full configuration is reminded in Fig. 5.7. The electrodes (anode and cathode) are highlighted in the figure. A Dirichlet boundary condition is imposed on the electrodes: Vcathode = 8 V,

Vanode = 0 V.

(5.23)

Homogeneous Neumann boundary conditions are imposed on all other walls, so that ∇V = 0. These Neumann boundary conditions correspond to the fact that the walls are modeled as insulators. Anode (0 V)

Cathode (+8 V)

1 mm

Figure 5.7: Slice of the full geometry. The electrodes (anode and cathode) are used to impose an electrical field in the micro-orifice. The full configuration is meshed using an unstructured tetrahedral mesh of approximately 18 million elements, with a refined region in the orifice (R/dx ≃ 5 in the orifice). It is presented in Fig. 5.8.

5.4. COMPUTATION OF THE ELECTRICAL POTENTIAL IN THE ABSENCE OF CELLS 129

Figure 5.8: Mesh of the full configuration.

Results Using the electrostatic solver, the electrical potential and electrical field are found in the full configuration. As seen in Fig. 5.9 and 5.10, most of the variation of these fields happen in the micro-orifice.

Figure 5.9: Electrical potential in the full configuration. Most of the variation of the field happens in the micro-orifice.

Figure 5.10: Electrical field in the full configuration. Most of the variation of the field happens in the micro-orifice. A zoom on the electrical potential and electrical field is presented in Fig. 5.11 to

130

CHAPTER 5. ELECTRICAL RESPONSE

better represent their evolution in the orifice. The variation is seen to be large in the micro-orifice and the electrical field is seen to display regions of higher intensity in the corners of the orifice. It is a consequence of the hydrofocalization. The rather homogeneous character of the electrical potential and electrical field suggest that, analogously to the reasoning of Chapter 4 about the velocity field, a reduced configuration can be used to reduce the computation time.

Figure 5.11: Electrical potential and electrical field in the micro-orifice of the full configuration. The higher values of the electrical field on the corners is a consequence of the geometry of the blood analyzer. The electrical current is computed in the full configuration by integrating the electrical field on a section of the micro-orifice and multiplying it by the conductivity of the fluid, which results in, If ull = 0.34 mA. (5.24) This value is in agreement with the experimental measurements provided by Horiba Medical and attests the quality of the computation.

5.4.2

Reduced configuration

In order to reduce computational costs due to the large variety of scales, a reduced configuration is used in all simulations. It is presented in this section. Mesh and boundary conditions The reduced domain was extracted from the full configuration along iso-surfaces of electrical potential in order to be able to simply apply uniform Dirichlet conditions on the boundaries of the reduced domain (left and right boundaries of Fig. 5.12). Homogeneous Neumann conditions are imposed on the other insulating walls. The uniform Dirichlet conditions are set as: Vlef t = 0.3 V Vright = 7.742 V. (5.25) It contains roughly 44 millions tetrahedral elements, which is reasonable since the electrostatic simulations are much faster than the dynamical simulations, because of

5.4. COMPUTATION OF THE ELECTRICAL POTENTIAL IN THE ABSENCE OF CELLS 131

Figure 5.12: Reduced configuration (light grey) compared to the full configuration (dark grey). White lines are iso-surfaces of the electrical potential intersections with the ynormal plane.

the static nature of these simulations. The micro-orifice is more refined than in the full configuration as it will be used to define the conductivity around the membrane and thus needs a finer mesh. The mesh in the configuration and in the micro-orifice, as well as the electrical potential is presented in Fig. 5.13

a)

b)

Figure 5.13: a) Reduced configuration for the electrostatic simulation, with mesh. b) Zoom on the orifice, with mesh.

Results In the case of the reduced configuration, the electrical potential and electrical field have the same variations as in the full configuration. Most of the variation happens in the micro-orifice, as seen in Fig. 5.14.

132

CHAPTER 5. ELECTRICAL RESPONSE

a)

b)

Figure 5.14: a) Electrical potential in the reduced configuration b) Electrical field in the reduced configuration.

As seen in Figure 5.12, the reduced configuration does not match perfectly the iso-surfaces of electrical potential, which may be the source of error in calculations performed with the reduced configuration. In order to quantify this error, Fig. 5.15 presents a comparison between the potential profiles on the centerline in the full and reduced configurations.

Potential (V)

Orifice

Longitudinal position (R) Figure 5.15: Comparison potential profiles along the centerline of the main axis of the full and reduced configurations.

133

5.5. ELECTRICAL PULSE GENERATED BY THE PASSAGE OF A CELL

A good agreement between the potential profile along the centerline of the main axis of the reduced and full configuration is found with a relative error of less than 0.1%, as presented on Fig. 5.15. Thus, all the electrostatic simulations performed in the following sections use this reduced configuration. The electrical current in the reduced configuration is found to be, Ireduced = 0.34 mA.

(5.26)

It is again in agreement with the value of 0.32 mA provided by Horiba Medical. These preliminary results show that the numerical framework and reduced configuration are suitable for the reproduction of the volume measurements in the blood analyzer.

5.5

Electrical pulse generated by the passage of a cell

In this section, the whole procedure followed in order to simulate the electrical pulse created by the passage of a red blood cell in the blood analyzer is detailed. An overview of the phenomenon and its characteristics are laid down in order to preview the approach used in the next section.

5.5.1

Effect of the passage of a cell

As Fig. 5.11 shows, an almost homogeneous electrical field is imposed in the orifice. The principle of the Coulter counter is that the passage of a cell disturbs the electrical field. More precisely, when a red blood cell travels trough the domain, its insulating nature locally disturbs the conductivity of the medium. This perturbation results in a variation of the overall resistance of the domain, through a resistive pulse whose amplitude is directly related to the cell volume. In the following sections, the amplitude of the resistive pulse created by the passage of the cell is denoted by ∆R, with: ∆R ∼ Vcell . (5.27)

5.5.2

Procedure of generation of the electrical pulse

The electrical pulse generated by the passage of a red blood cell is reconstructed from the simulation of the dynamics of this red blood cell in the blood analyzer. For each successive position of the red blood cell in the blood analyzer, the geometry of the red blood cell is used to impose the conductivity field in the analyzer. This results, through the electrostatic solver, in a value for the resistance of the system, which is then compared to the resistance of the system in the absence of red blood cells. The successive positions and shapes of the red blood cell leads to the generation of a value for R along the trajectory of the red blood cell. This results in the reconstruction of the electrical pulse as sketched in Fig. 5.16.

134

CHAPTER 5. ELECTRICAL RESPONSE

∆R (Ω)

Vcell

Time Figure 5.16: Resistive pulse schematic representation (right) and successive positions of a red blood cells in flow. (left)

5.5.3

Range of the study

In reality, the relationship of Eq. (5.27) that states that the electrical pulse amplitude is directly related to the cell volume is true only for a population of cells that all travel through the orifice with the same position, orientation and shape. To clarify, the perturbation of the electrical field depends on the projected surface that is obstructed by the red blood cell in the orifice. This projected surface depends on the orientation and shape of the red blood cell (see Fig. 5.17). The position of the cell can also create undesired effects if it happens to flow in the viscinity of the insulating walls of the orifice.

a)

b)

Figure 5.17: Schematic representation of the influence of the orientation of the cell on the electrical field streamlines. From Is`ebe and N´erin [96].

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

135

Thus, this results in the possibility for a single red blood cell with a fixed volume Vcell to generate more than one value of ∆R when measured in the blood analyzer, depending on its position, orientation and shape. This leads to the necessity of performing numerous simulations with different initial positions and orientations to understand the importance of this effect on the volume measurement in the blood analyzer. This will be discussed in more details in section 5.6.

5.5.4

Test case with a rigid particle

This test case is performed by Horiba Medical in the calibration procedure of their devices. It consists in putting a rigid sphere made of latex in the blood analyzer. Since its volume is known a priori, the device can be calibrated by using this measurement. It is a reliable procedure because of the spherical symmetry of the particle, which ensures the independence of the measurement on the orientation and shape of the particle. This calibration procedure is reproduced using YALES2BIO. The latex sphere is of diameter Dlatex = 5 µm and is found to generate a electrical pulse of amplitude of ∆R = 10.22 Ω traveling through the blood analyzer, which is in agreement with internal reports from Horiba Medical.

5.6

The probability density function of volume measurements

The post-processing step aims at constructing the probability density function of the resistive pulse amplitude generated by a single cell and the volume probability density function for a population of red blood cells. In order to obtain this statistical results, an very large number of results would be required (at least 2 × 104 ). Obviously, full numerical simulations using the YALES2BIO solver for this amount of simulations are not possible because of the potential computational cost. Thus, a fitting procedure is needed to perform the Monte-Carlo simulation in the input parameters (r0 , θ0 , ϕ0 ) domains and increase the number of resistive pulse amplitude results available for the statistical treatment of the results.

5.6.1

Why a statistical approach

A statistical approach is performed to generate the results of the volume measurement of a red blood cell of fixed volume in the blood analyzer, depending on its orientation and position. As stated previously, the amplitude ∆R of the electrical pulse generated by the passage of a red blood cell depends on its orientation and position in the orifice. The orientation and position of the red blood cell in the orifice depends on the initial orientation (θ0 , ϕ0 ) and position r0 of the red blood cell (at z0 = −R). However, the initial orientation and position of the cell at z0 are not known. Thus, the repartition of r0 ,

136

CHAPTER 5. ELECTRICAL RESPONSE

θ0 and ϕ0 must be supposed a priori to compute the probability density function of ∆R. The following sections detail the complete procedure that was used to generate the probability density function of ∆R for a red blood cell of fixed volume. It goes as follows: • The probability density functions of r0 , θ0 and ϕ0 are defined. (Section 5.6.2) • Numerous numerical simulations are performed with (r0 , θ0 , ϕ0 ) as input parameters for a red blood cell of fixed volume. These simulations are used to generate the value of ∆R corresponding to each triplet (r0 , θ0 , ϕ0 ) in YALES2BIO. (Section 5.6.3) • The results of ∆R = f (r0 , θ0 , ϕ0 ) are fitted to allow the drawing of any triplet of values of (r0 , θ0 , ϕ0 ) and the computing of the corresponding value of ∆R. (Section 5.6.4) • Using a Monte-Carlo procedure, numerous triplets (r0 , θ0 , ϕ0 ) are drawn following their respective probability density function. The resulting ∆R are computed using the fitted function, allowing the generation of the probability density function of ∆R. (Section 5.6.5)

5.6.2

Random variables: initial position and orientation

In order to perform the Monte-Carlo simulations, it is necessary to select values for (r, θ, ϕ). An a priori probability density has to be imposed on r, θ and ϕ. Obviously the final results depend on this a priori distribution. In the absence of knowledge of how cells are positioned, we make the assumption of uniform distribution of position and orientation for cells located at z0 = −R and r ≤ 15 µm. We will see that this assumption results in non-trivial distributions or r and θ. Definition of the initial position and orientation A single red blood cell is used in all the red blood cells dynamics simulations, as presented in Section 4.1.4. The dynamics simulations are not presented here as they have been extensively described in Chapter 4. However, the definition of the initial location and orientation is reminded. The initial position of the red blood cell in the numerical domain is still defined by its position z0 along the blood analyzer axis (z-axis), which is fixed: z0 = −R, (5.28) with R = 25 µm, the radius of the aperture (Fig. 5.18). This position was chosen so that most of the acceleration phase is computed at a relatively reduced cost. From this initial coordinate, three varying parameters arise in order to describe the red blood cell position and orientation:

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

137

• r, the initial position along the x-axis. • θ and ϕ, two angles used to define the red blood cell orientation. As it has been presented in Section 4.1.4, the acceptable range for r is defined using the streamlines computed from the exit of the sample injection tube and the symmetries of the problem. Red blood cells are considered to more or less follow the streamlines. Shear-induced diffusion is neglected here, especially because the hematocrit in the sample is very low (0.01% − 0.05%) due to preliminary dilution of the blood sample. Thus, due to the hydrofocalization, the red blood cells will travel in the region bounded by those streamlines, which results in the following range for r, at z = −R (see Fig. 5.18): • r ∈ [0; 15] µm.

Figure 5.18: Reduced configuration around the aperture with streamlines. The local frame introduced in Section 4.1.4 is recalled on Fig. 5.19; it is used to define angles θ and ϕ. The local frame associated to the red blood cell is also briefly recalled: Figure 5.20 shows the ~xRBC and ~zRBC axes which are centered on the center of mass of the red blood cell and aligned with the long and small radii of the cell, respectively. These two axes are chosen to be in the (~xstream , ~zstream ) plane, so that ~yRBC = ~ystream completes the direct (~xRBC , ~yRBC , ~zRBC ) orthonormal basis.

138

CHAPTER 5. ELECTRICAL RESPONSE

~ystream ~zstream ~ystream ~zstream

~xstream

~xstream

Figure 5.19: Definition of axes ~xstream , ~ystream (normal to the plane), ~zstream in the reduced configuration.

Figure 5.20: Representation of the red blood cell local frame and definition of θ and ϕ.

As seen on Fig. 5.20, the angle θ, which is the amount of rotation with respect to the direction of ~zstream , is defined as, cos θ = ~zRBC . ~zstream .

(5.29)

The second angle, ϕ, is the amount of rotation around ~zstream . cos ϕ = ~xp . ~xstream ,

(5.30)

with ~xp , the projection of ~xRBC on the plane formed by ~xstream and ~ystream = ~y . Considering the geometrical symmetries of the blood analyzer itself, as well as the ones of the red blood cell. The resulting ranges for θ and ϕ are: • θ ∈ [−90◦ ; 90◦ ].

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

139

• ϕ ∈ [0◦ ; 90◦ ]. For more details about the initial state of the red blood cells, the reader is referred to Section 4.1.4. Probability density function of the initial position Previously in this thesis, it has been stated that the red blood cells can arrive to the entrance of the aperture at various positions along the x-axis (perpendicular to the blood analyzer axis). Since this problem is 3-dimensional, the actual domain where a red blood cell can pass through at z0 is a disc of center (0, 0, z0 ) and radius Rr = 15 µm (see Fig. 5.21). ψ

Figure 5.21: Schematic representation of the possible values for r. If r is randomly selected between 0 and Rr and ψ randomly selected between 0◦ and the sample points density will not be uniform, as it is seen on Fig. 5.22a. In order to randomly choose a value for r on this disc with a uniform density of selected points on the disc, the probability for a specific value of r to be drawn has to be computed. A uniform distribution means that the surface density should be uniform. It is the probability to choose r inside Ω, a full crown of surface dSΩ on the disc of surface Stot , which can be written as, dSΩ PΩ = , (5.31) Stot with, Stot = πRr2 . (5.32) 360◦

Considering a surface element on the disc and using polar coordinates, dSΩ = rdrdψ,

(5.33)

rdrdψ . πRr2

(5.34)

the probability reads PΩ =

140

CHAPTER 5. ELECTRICAL RESPONSE

Since ψ is uniformly distributed over the considered range, the cumulative distribution function Fr can be found by integrating this probability, Fr =

Z r 2tdt 2 t=0 Rr

=

r2 , Rr2

(5.35)

This allows to perform a uniform random selection on Fr , computing the value for r using: p (5.36) r = R r Fr .

Randomly and uniformly selecting the cumulative distribution function Fr in [0, 1] ensures a uniform density of selected initial positions for the cell, as shown on Fig. 5.22b.

b)

y (µm)

a)

x (µm)

x (µm)

Figure 5.22: Sample points selected on the disc of radius R = 15 µm for a) A random selection of r ∈ [0; 15] µm and ψ in the disc, b) A random selection following Eq. (5.36).

Probability density function of the initial orientation The orientation of the red blood cell can be defined by imposing the direction of the unit vector zRBC . It is defined by a triplet of coordinates lying on a sphere centered on the value of r for the considered case and of radius 1. Due to symmetries, only half a sphere is needed. The classical spherical coordinates system is considered, with θ′ the co-latitude and ϕ′ the azimuth angle. In this coordinates system, the bounds for θ′ and ϕ′ are not the same as the one described earlier, notably, θ′ does not take negative values in this case, and the range of possible values for ϕ′ is greater that the one considered in section 4.1.4. If θ′ and ϕ′ are uniformly sampled in their own ranges, the resulting selected vectors tips density on the surface of the sphere is non-uniform, as seen in Fig. 5.24a. Using the same method as for the selection of r, θ′ and ϕ′ are

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

141

chosen so that the vectors tips on the surface of the half-sphere of center r and radius 1 in the frame (xstream , ystream ) display a uniform distribution.

θ′

ϕ′

Figure 5.23: Definition of angles θ′ and ϕ′ . The probability to choose θ′ and ϕ′ in a surface element of the sphere is found as, PΩ =

sin θ′ dθ′ dϕ′ . 4π

(5.37)

The cumulative distribution function Fθ′ is then computed as: Fθ ′ =

Z t=θ′ sin tdt t=0

2

=

1 − cos θ′ , 2

(5.38)

since ϕ′ is again uniformly distributed over the considered range. Thus, the value of θ′ can be computed from Fθ′ : (5.39) θ′ = arccos(1 − 2Fθ′ ). Uniformly sampling Fθ′ in [0, 1] ensures a uniform density of selected orientations on the sphere describing the possible initial orientations for the cell. After this random selection is performed, a transformation is applied on θ′ and ϕ′ in order to retrieve the more convenient definitions and bounds for θ and ϕ introduced in Section 4.1.4. It is done using:  ′  ϕ = ϕ′ if ϕ′ < 90◦  θ=θ, . (5.40)   ′ ◦ ′ ′ ◦ θ = −θ , ϕ = ϕ − 90 if ϕ > 90 These definitions allow a better visualization of the results and better understanding of the observed symmetries.

142

CHAPTER 5. ELECTRICAL RESPONSE

b)

ystream

a)

xstream

xstream

Figure 5.24: zRBC vectors tips on the surface of the half-sphere, top-view for a) A random selection of θ′ and ϕ′ , b) A random selection following Eq. (5.39) for θ′ and uniformly random for ϕ′ .

The probability density functions presented in this section are chosen so that the initial orientation and position of the cell is uniformly random at z0 = −R. This hypothesis can be discussed, as chapter 4 shows that the red blood cells reorient themselves along their travel inside the blood analysis. Thus, the probability of initial orientations near θ0 = 90◦ (red blood cells’ long axes aligned with the direction of the flow) should be higher than otherwise. However, this probability density cannot be derived from chapter 4 easily and would be arbitrary, thus the choice has been made to consider a uniformly random initial orientation and position.

5.6.3

Raw results

The presented numerical setup allows to perform and post-process numerous simulations in order to find the dependence of ∆R on (r0 , θ0 , ϕ0 ). The performed cases have been determined using a Latin Hypercube Sampling method [95]: each variable from the (r, θ, ϕ) triplet is discretized in a specified number of sub-ranges. In each of these sub-ranges, a value for each parameter is chosen randomly. This allows to probe the resistance pulse amplitude for any possible case in the defined ranges for (r0 , θ0 , ϕ0 ). In this work, r0 , θ0 and ϕ0 have been discretized as follows: • Variable r0 , nr = 3 : r0 ∈ [5n; 5(n + 1)] µm with n = 0, 1 or 2. • Variable θ0 , nθ = 10 : θ0 ∈ [18n◦ ; 18(n + 1)◦ ] with n = −5, −4, ...4. • Variable ϕ0 , nϕ = 5 : ϕ0 ∈ [18n◦ ; 18(n + 1)◦ ] with n = 0, 1, ..., 4.

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

143

90 72 54

ϕ0 ( ) ◦

36 18 0 90 72 54

15

36 18

θ0 (◦ )

0 −18 −36 −54 −72 −90

10 5 0

r0 (µm)

Figure 5.25: Representation of all the cases performed in the (r0 , θ0 , ϕ0 ) space.

Figure 5.25 illustrates the different cases performed. The grid lines indicates the different ranges considered for each parameter for the Latin Hypercube Sampling. More than 500 simulations have been performed over the full range of all parameters. Note that the centered cases (r = 0 µm) have been omitted in Fig. 5.25 for readability, since these cases can be expanded for all values of θ0 and ϕ0 thanks to the geometry and red blood cell symmetries. The obtained results are the resistance pulses amplitude for different values of the triplet (r0 , θ0 , ϕ0 ): ∆R = f (r0 , θ0 , ϕ0 ).

(5.41)

The results are first presented as a function of one parameter, for all values for the other two parameters for the sake of readability. Figures (5.26), (5.27) and (5.28) show the variation of the resistance pulse amplitude ∆R as a function of θ0 , r0 , and ϕ0 , respectively. A strong dependence of ∆R on θ0 is found (Fig. 5.26). On the other hand, no clear dependence of ∆R on both ϕ0 and r0 can be seen from Fig. 5.27 or 5.28.

144

∆Rmax (Ω)

CHAPTER 5. ELECTRICAL RESPONSE

θ0 (◦ )

∆Rmax (Ω)

Figure 5.26: Resistance pulse amplitude versus θ0 for all values of ϕ0 and r0 .

r0 (µm) Figure 5.27: Resistance pulse amplitude versus r0 for all values of θ0 and ϕ0 .

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

∆Rmax (Ω)

145

ϕ0 (◦ ) Figure 5.28: Resistance pulse amplitude versus ϕ0 for all values of θ0 and r0 .

5.6.4

Fitting

Polynomial fitting As seen in the previous section, the value of ∆R mostly depends on the value of θ0 . Only the dependence on θ0 is considered since the other two parameters, ϕ0 and r0 , are found not to significantly influence the results on the value of ∆R. If a polynomial of degree Np is considered, ∆Rpoly = f (θ) reads, ∆Rpoly =

Np X

ai θ0i .

(5.42)

i=0

The convenient degree for the polynomial fit is found using the root mean-square error and the correlation coefficient. Polynomial with degrees Np ∈ [20; 50] produce a satisfying fitting of the results. Even if any value for Np between Np = 16 and Np = 50 correctly predicts the behavior of ∆R, the optimal and chosen value for Np is found to be Np = 49 according to Fig. 5.29, with a correlation coefficient of 0.97 and a root mean-square error of 0.45. In addition, a polynomial fitting of degree Np = 49 allows the capture of higher values of ∆R.

146

RMSE / Correlation coefficient

CHAPTER 5. ELECTRICAL RESPONSE

Np Figure 5.29: Root mean-square error (+) and correlation coefficient (◦) versus polynomial degree.

The scatter plot for the Polynomial fitting is shown on Fig. 5.30. The fitting is seen to correctly predict the behavior of the value of ∆R without the dependency on r0 and ϕ0 . Lower values of the electrical pulse amplitude ∆R are more correctly captured by the polynomial fitting.

5.6.5

Probability density function

From the probability density functions of (r0 , θ0 , ϕ0 ) and the polynomial fitting, N (N ≫ 1) results for ∆R are generated following the supposed distribution of (r0 , θ0 , ϕ0 ) at z0 = −R. The results are then sorted in different classes depending on the value of ∆R N 1 X P (∆Rinf ≤ ∆R ≤ ∆Rsup ) = F (∆Ri ), (5.43) N i=1 where, F (∆Ri ) =

   1 if ∆Rinf ≤ ∆Ri ≤ ∆Rsup  

0 if not

,

(5.44)

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

∆Rpoly (Ω)

147

∆RCF D Figure 5.30: Scatter plot of ∆Rpoly (with Np = 49) versus ∆RCF D .

Polynomial

0,8

Probability YY

0,6

0,4

0,2

0 10

11

12

13

14 XX

15

16

17

18

∆R (Ω) Figure 5.31: Probability density function for the resistance pulse amplitude.

From Table 5.2, the 95-th percentile allows the characterization of the accuracy provided by a blood analyzer with hydrofocalization: 95% of the measurement will show a resulting electrical pulse amplitude of around 13.27 Ω or less. Thus if the actual volume of the cell is considered to be the value of ∆R for a centered cell with θ0 = 90◦

148

CHAPTER 5. ELECTRICAL RESPONSE Statistics on the polynomial approximation Statistical quantity Mean value Median value Standard deviation Skewness Mean error 95-th percentile

Polynomial 11.04 10.71 1.26 4.13 4.22% 13.27 (25%)

Table 5.2: Main statistics of ∆R in the case of a polynomial fitting of the results

(∆R = 10.60 Ω), the error on the volume measurement will be less than 25% in 95% of the cases. It is also found that the mean error on the measurement is around 4%. Application to the generation of a complete volume probability density function In this section, results from the computed simulations are used to mimic the measurement of the volume probability density function in a Coulter counter. The volume probability density function is the distribution of measured volumes of a red blood cell population that has traveled through the blood analyzer. The previous analysis was performed for a unique value of volume. Here, we first assume that in a fictive blood sample, there is a certain volume distribution. The volume distribution is supposed to be Gaussian, with a mean value of 93.54 fL. and a standard deviation of 0.13 [? ]. This probability density function thus reads: −(V −µ)2 1 p.d.f (V ) = √ e 2σ2 , σ 2π

(5.45)

σ = 0.13 and µ = 93.54 fL.

(5.46)

with,

This is the assumed ’real’ distribution in a sample. The aim is now to apply the procedure described earlier for one specific volume to this distribution of volume and generate the probability density function of ∆R for the whole sample: f∆R is now function of (r, θ, ϕ, V ). Thanks to the Monte-Carlo simulations, 20 millions values for (r, θ, φ) are generated and used in this section. The 20 millions values for (r, θ, φ) provide a collection of ∆R0 values generated from the polynomial fit of section 5.6.4. For each value of ∆R0 , a volume is selected following the physiological distribution of red blood cells volumes of Eq. (5.45). A linearity hypothesis is then used to generate the corresponding value of ∆R, it writes: ∆R =

V ∆R0 , V0

(5.47)

5.6. THE PROBABILITY DENSITY FUNCTION OF VOLUME MEASUREMENTS

149

where V0 is the volume of the red blood cell used in all simulations. When this operation is repeated over all of the values of ∆R0 , a probability density function for ∆R of the whole sample is found. Finally, a calibration procedure is introduced to retrieve the measured volume distribution: Vm = α∆R,

with

α=

Vc . ∆Rc

(5.48)

In this equation, the chosen values for ∆Rc and Vc define the calibration. Two options were considered: • Exact calibration: ∆Rc is chosen as the resistive pulse amplitude for a centered cell with θ = 90◦ and Vc = V0 .

Probability

• Latex sphere calibration: it follows the actual calibration process made by Horiba medical experts. ∆Rc is chosen as the resistive pulse amplitude of a latex sphere of diameter Dlatex = 5 µm, ∆Rc = 10.22 Ω, with Vc = Vlatex × f = 97.5 fL, where f is a shape factor to account for the shape difference between the sphere and a cell. Here, f = 1.5 as provided by Horiba Medical internal reports.

Volume (fL) Figure 5.32: Volume probability density function for different calibration values. Input VD is the input volume distribution defined in Eq. (5.45), Observed VD are the observed volume distributions with the two calibration methods described above. From Fig. 5.32, the observed volume distribution depends on the calibration. All the statistical quantities reported in Table. 5.3 show an increase compared to the input volume distribution, which creates a systematic overestimation of the Mean Corpuscular Volume (MCV) and Red Blood Cell Distribution Width (RDW). Obviously, the exact

150

CHAPTER 5. ELECTRICAL RESPONSE

calibration does not create more error than the one measured for a single cell probability density function found in Fig. 5.31, around 4% on the MCV. The latex sphere calibration, where a latex sphere is used as a reference particle as described previously, introduced an error of around 12% on the MCV. Both calibration methods show an increase of around 15% on the observed RDW, compared to the actual input RDW. Note also that the right skewness in the observed volume distributions is in agreement with the literature [28]. Statistical quantity Mean value Median value Standard deviation Skewness Red Blood Cell Distribution Width Mean error

Input VD

OVD (Exact)

OVD (Latex)

93.546 93.538 0.130 0.0 13% 0%

97.106 96.064 0.161 0.816 15% 3.80%

104.95 103.820 0.167 0.817 15% 12.19%

Table 5.3: Quality of the measurement of the volume distribution in the blood analyzer.

5.7

Discussion

These results show that the numerical tool developed in this thesis is able to reproduce the measured volume distribution from an input healthy red blood cells volume distribution and the bank of results used to generate the fitting function of ∆R. Thus, this numerical tool can be used more widely for different configurations or calibration procedures. It is robust enough to handle complicated geometries and flows. It is important to note that the previous results depend on all of the assumptions made for the construction of the probability density functions. Namely, the distribution of θ0 , found by considering a uniform probability for the orientation of the cell is strong, as it neglects the reorientation phenomena that may happen before the deposition at z0 = −R. This kind of assumption can introduce errors that are most likely found in non-hydrofocalized systems.

Chapter

6

Conclusion Chapter contents 6.1 6.2

6.3

Main results . . . . . . . . . . . . Perspectives . . . . . . . . . . . . 6.2.1 Alternative geometries . . 6.2.2 Red blood cell mechanics 6.2.3 Cell-cell interaction . . . . 6.2.4 Electrical modeling . . . . Concluding remark . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

151 153 153 153 154 154 154

This section aims at summarizing the work done in this thesis. Then, it discusses some perspectives of this work and the possibilities of further research offered by the results. Finally, a concluding remark is presented.

6.1

Main results

This thesis presents a numerical and analytical study of the dynamics of red blood cells under flow. This thematic has known a tremendous development in the last years. However, this work is original compared to the literature on the topic, because it focuses on the behavior of red blood cells in flows with a strong elongational character, characteristic of what is found in hydrofocalized blood analyzers based on the Coulter principle. This thesis establishes results on the reorientation of particles in elongational flows through the development of analytical models. The analytical developments of Chapter 3 allow a better understanding of the physical effects taking part in the reorientation process of particles in elongational flows: the Keller & Skalak and Abkarian, Faivre & Viallat models, originally developed for pure shear flows, are re-derived in this context. The former considers a rigid ellipsoidal particle with an imposed membrane velocity. The latter adds the contribution of the 151

152

CHAPTER 6. CONCLUSION

elasticity of the membrane to the Keller & Skalak model. Both are used to predict the orientation of the cell in an elongational flow. It is found that in the limit of infinite capillary number values, the Abkarian, Faivre & Viallat (AFV) model is found to converge towards the Keller & Skalak model (K&S). In the limit of vanishing capillary number, the AFV model converges towards Jeffery’s theory for rigid ellipsoidal particles. Finally, both the AFV and K&S models are found to converge towards Jeffery’s theory for higher values of the viscosity ratio (µext /µint ). This latter analytical model shows a clear reorientation of the cell with the axis of elongation in all cases. The viscosity ratio has an important influence on the reorientation phenomenon, slowing it down for higher values. The contribution of the elasticity introduced through the AFV model is seen to slow down the reorientation as it accounts for the storing of the elastic energy in the membrane and induces changes in the membrane motion. Numerical simulations of the dynamics of red blood cells in a configuration characteristic of hydrofocalized Coulter counters are presented. They show the reorientation phenomenon of centered red blood cells in the blood analyzer. The results clearly show this trend even in the case of a ”late” upstream deposition in the blood analyzer. The first contribution to this reorientation process is the velocity field, the elongational character of the flow forcing the red blood cells to align with the flow. Then, the deformation is seen to contribute in lower angles cases, as the initial deformation has a greater impact on the reorientation. Another important effect is the viscosity ratio: it is seen to play an important role in the reorientation process. Cells with lower viscosity ratios are seen to align themselves with the direction of the flow more rapidly than cells with higher viscosity ratios, in agreement with analytical models. The viscosity ratio is also observed to have more effect on lower angles case, as it modifies the cell ability to resist deformation in the initial moments of the simulation. Finally, a prediction of the electrical measurement of the volume of red blood cells in an industrial blood analyzer is performed. The order of magnitude of the electrical pulse found in experimental conditions is correctly retrieved. The dependence of the electrical pulse on the initial position and orientation of the flowing cell is investigated: the whole range of the input parameters are probed in order to accurately fit the results from an extensive numerical database. This fitting is later used in a statistical approach to produce a probability density function of the electrical pulse amplitude, which is the quantity from which volume is measured. It is found that the fact that red blood cells align themselves with the axis of the blood analyzer results in a homogenization of the electrical response, where a large majority of cells are seen to produce an electrical pulse of similar amplitude. From this probability density function, the measure of a complete healthy red blood cell population is simulated. It is found that discrepancies in the cell alignment would add a right skewness and a tail to the distribution of volumes, thus over-estimating the number of higher volume red blood cells. These results stress the advantages of hydrofocalization in the context of red blood

6.2. PERSPECTIVES

153

cell volumetry: the first effect of hydrofocalization is to force the cells to flow through the center of the orifice, thus avoid them to interact with the complex electrical field that is present near the walls of the orifice. A second effect, probably even more important, is the reorientation of the cells. The latter not only flow in the center of the orifice, but the strong elongational flow makes them accelerate and reorient to be aligned with the flow, whatever their orientation before the acceleration. Hydrofocalization is thus extremely beneficial in terms of robustness and accuracy of volume measurements.

6.2 6.2.1

Perspectives Alternative geometries

It would be of great interest to perform the same study on systems that are not hydrofocalized, as it would allow a characterization of the importance of its effect. Also, the optimization of the current geometry, namely in terms of sheathing flows and their flow rates, would be interesting and preview simple improvements of an actual device. In order to understand the importance of the velocity field and possibly to increase the speed of the reorientation phenomenon, alternative geometries can be used. Notably, the hydrofocalization can be tuned to impose a stiffer velocity field that will flip the red blood cells in alignment with the flow in a shorter time. In addition, the velocity field could be manipulated in order to avoid undesirable effects such as close cells doublets which contribute to the overestimation of the number of higher volume red blood cells in experimental measurements. If the red blood cells membrane was found to break during the measurement, the modification of the velocity field could also help maintain the cells integrity. This kind of studies could help the design of more reliable devices and ensure non-destructive processes.

6.2.2

Red blood cell mechanics

The red blood cell mechanics described and used for the simulations in this work are of satisfying accuracy. As the validation test cases demonstrated it, all of the main features and behaviors of the red blood cell are retrieved. Still, some of the extreme cases are not well reproduced, and some effects could thus be added to the cell: the visco-elasticity of the membrane that can contribute to the dynamics of red blood cells, or the thickness of the membrane that should offer a more natural way of modeling the curvature resistance of the cell membrane are examples. A completely different modeling framework could also be considered, either a less refined one where it would be interesting to see if the features of the reorientation phenomenon are retrieved or a more refined one that might let more complicated effects play a role in the measurement process. For example, the modeling of the bilayer structure of the red blood cells or the viscosity of the membrane viscosity could slow down the deformation and reorientation.

154

6.2.3

CHAPTER 6. CONCLUSION

Cell-cell interaction

The interaction between cells might have an important effect on the orientation phenomenon and in the measurement of the cell volume in the blood analyzer. The interaction of red blood cells happening in the injection tube notably, could have a significant effect on the resulting orientation in the micro-orifice. Aside from this indirect effect, some of the recent experimental measurements show that two red blood cells can stay very close during their travel in the blood analyzer. The numerical simulation of this kind of event is of interest, as it could help its detection from the experimental electrical pulse. Moreover, it could provide insights on how to design the blood analyzer in order to avoid this situation.

6.2.4

Electrical modeling

The current framework for the electrical modeling of red blood cells is simplified, in order to allow fast computations and easier understanding of the results. Still, as it as been discussed in the first part of Chapter 5, many effects can arise from the influence of an electrical field on the cell membrane. It would be of interest to model more complicated electrical effects. Notably, the dynamical modification of the electrical field could be considered and investigated further if the static approximation is no longer made. Also, the interdependence of the generated electrical field and the resulting field in the membrane could have effects on the cell shape, orientation and trajectory.

6.3

Concluding remark

An important result of this thesis is the actual numerical tool and framework for the simulation of the dynamics of red blood cells in blood analyzers and the associated volume distribution measurement. It can be used with different configurations and will certainly allow a faster study of the flow and dynamics of red blood cells in such devices. The reproduction of the electrical measurement also allows the characterization of the principal components of this measurement, and how the future devices should be designed in order to ensure satisfying accuracy. This suggests that YALES2BIO could be an important numerical tool for designing medical devices but also for the understanding of the dynamics of cells in ex-vivo contexts.

Bibliography [1]

P. Aalto, K. H¨ ameri, P. Paatero, M. Kulmala, T. Bellander, N. Berglind, L. Bouso, G. Casta˜ no Vinyals, J. Sunyer, G. Cattani, A. Marconi, J. Cyrys, S. Von Klot, A. Peters, K. Zetzsche, T. Lanki, J. Pekkanen, F. Nyberg, B. Sj¨ovall, and F. Forastiere. Aerosol particle number concentration measurements in five european cities using TSI-3022 condensation particle counter over a three-year period during health effects of air pollution on susceptible subpopulations. J. Air Waste Manag. Assoc., 55:1064–1076, 2005.

[2]

M. Abkarian, M. Faivre, R. Horton, K. Smistrup, C. A. Best-Popescu, and H. A. Stone. Cellular-scale hydrodynamics. Biomed. Mater., 3(034011), 2008.

[3]

M. Abkarian, M. Faivre, and A. Viallat. Swinging of red blood cells under shear flow. Phys. Rev. Lett., 98(188302), 2007.

[4]

M. Abkarian, L. Lanotte, J.-M. Fromental, S. Mendez, D. A. Fedosov, G. Gompper, J. Mauer, and V. Claveria. A new look on blood shear thinning. In 68th Annual Meeting of the APS Division of Fluid Dynamics, Boston, MA., 22-24 November 2015, 2015.

[5]

M. Abkarian and A. Viallat. Vesicles and red blood cells in shear flow. Soft Mat., 4:653–657, 2008.

[6]

C. K. Aidun and J. R. Clausen. Lattice-Boltzmann method for complex flows. Ann. Rev. Fluid Mech., 42:439–472, 2010.

[7]

C. K. Aidun and J. R. Clausen. Supplementary material: Lattice-Boltzmann method for complex flows. Supplemental Material to Annu. Rev. Fluid. Mech. 2010. 42:439–472., 2010.

[8]

K. Asami. Characterization of biological cells by dielectric spectroscopy. Non-Cryst. Sol., 305:268–277, 2002.

[9]

K. Asami. Characterization of heterogeneous systems by dielectric spectroscopy. Prog. Polym. Sci., 27:1617–1659, 2002.

[10]

K. Asami. Effects of membrane disruption on dielectric properties of biological cells. J. Phys. D: Appl. Phys., 39, 2006.

[11]

K. Asami. Simulation of dielectric spectra of erythrocytes with various shapes. J. Phys. D: Appl. Phys., 42:135503, 2009. 155

J.

156

BIBLIOGRAPHY

[12]

P. Bagchi. Mesoscale simulation of blood flow in small vessels. 92:1858—1877, 2007.

Biophys. J.,

[13]

P. Bagchi, P. C. Johnson, and A. S. Popel. Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow. J. Biomech. Eng., 127(4):1070–1080, 2005.

[14]

P. Bagchi and R. M. Kalluri. Dynamics of nonspherical capsules in shear flow. Phys. Rev. E, 80(016307), 2009.

[15]

P. Bagchi and R. M. Kalluri. Rheology of a dilute suspension of liquid-filled elastic capsules. Phys. Rev. E, 81(056320), 2010.

[16]

P. Bagchi and R. M. Kalluri. Dynamic rheology of a dilute suspension of elastic capsules: effect of capsule tank-treading, swinging and tumbling. J. Fluid Mech., 669(498–526), 2011.

[17]

D. Barth`es-Biesel. Motion of a spherical microcapsule freely suspended in a linear shear flow. J. Fluid Mech., 100(4):831–853, 1980.

[18]

D. Barth`es-Biesel. Capsule motion in flow: Deformation and membrane buckling. Comp. Rend. Phys., 10(8):764–774, 2009.

[19]

D. Barth`es-Biesel, A. Diaz, and E. Dhenin. Effect of constitutive laws for twodimensional membranes on flow-induced capsule deformation. J. Fluid Mech., 460:211–222, 2002.

[20]

D. Barth`es-Biesel and J. M. Rallison. The time-dependent deformation of a capsule freely suspended in a linear shear flow. J. Fluid Mech., 113:251–267, 1981.

[21]

D. Barth`es-Biesel and H. Sgaier. Role of membrane viscosity in the orientation and deformation of a spherical capsule suspended in shear flow. J. Fluid Mech., 160:119–135, 1985.

[22]

J. C. Baygents, N. J. Rivette, and Stone H. A. Electrohydrodynamic deformation and interaction of drop pairs. J. Fluid Mech., 368:359–375, 1998.

[23]

J. Beaucourt, F. Rioual, T. S´eon, T. Biben, and C. Misbah. Steady to unsteady dynamics of a vesicle in a flow. Phys. Rev. E, 69(011906), 2004.

[24]

L. I. Berge, T. J¨ossang, and J. Feder. Off-axis response for particles passing through long apertures in Coulter-type counters. Meas. Sci. Technol., 1:471–474, 1990.

[25]

T. Biben and C. Misbah. Tumbling of vesicles under shear flow within an advectedfield approach. Phys. Rev. E, 67(031908), 2003.

BIBLIOGRAPHY

157

[26]

G. Boedec, M. Jaeger, and M. Leonetti. Settling of a vesicle in the limit of quasispherical shapes. J. Fluid Mech., 690:227–261, 2012.

[27]

G. Boedec, M. Leonetti, and M. Jaeger. 3D vesicle dynamics simulations with a linearly triangulated surface. J. Comput. Phys., 230:1020–1034, 2011.

[28]

M. O. Breitmeyer, E. N. Lightfoot, and W. H. Dennis. Model of red blood cell rotation in the flow toward a cell sizing orifice. Biophys. J., 11:146–157, 1971.

[29]

J.-M. Charrier, S. Shrivastava, and R. Wu. Free and constrained inflation of elastic membranes in relation to thermoforming non-axisymmetric problems. The Journal of Strain Analysis for Engineering Design, 24(2):55–74, 1989.

[30]

S. Chen and G. D. Doolen. Lattice-Boltzmann method for fluid flows. Ann. Rev. Fluid Mech., 30:329–364, 1998.

[31]

S. Chien. Shear dependence of effective cell volume as a determinant of blood viscosity. Science, 168:977–979, 1970.

[32]

C. Chnafa, S. Mendez, and F. Nicoud. Image-based large-eddy simulation in a realistic left heart. Comput. Fluids, 94:173–187, 2014.

[33]

C. Chnafa, S. Mendez, F. Nicoud, R. Moreno, S. Nottin, and I. Schuster. Imagebased patient-specific simulation: a computational modelling of the human left heart haemodynamics. Comput. Meth. Biomech. Biomed. Eng., 15(supp1):74–75, 2012.

[34]

A. Chorin. Numerical solution of the Navier–Stokes equations. 22:745–762, 1968.

[35]

J. R. Clausen, D. A. Reasor, and C. K. Aidun. The rheology and microstructure of concentrated non-colloidal suspensions of deformable capsules. J. Fluid Mech., 685:202–234, 2011.

[36]

D. Cordasco and P. Bagchi. Orbital drift of capsules and red blood cells in shear flow. Phys. Fluids, 25(091902), 2013.

[37]

D. Cordasco and P. Bagchi. Intermittency and synchronized motion of red blood cell dynamics in shear flow. J. Fluid Mech., 759:472–488, 2014.

[38]

D. Cordasco, Yazdani, and P. Bagchi. Comparison of erythrocyte dynamics in shear flow under different stress-free configurations. Phys. Fluids, 26(041902), 2014.

[39]

G.-H. Cottet and E. Maitre. A level set method for fluid-structure interactions with immersed surfaces. Math. Mod. Meth. App. Sc., 16(3):415–438, 2006.

Math. Comp.,

158

BIBLIOGRAPHY

[40]

G.-H. Cottet, E. Maitre, and T. Milcent. Eulerian formulation and level set models for incompressible fluid-structure interaction. ESAIM: Math. Mod. Num. Anal., 42:471–492, 2008.

[41]

W. H. Coulter. Means for counting particles suspended in a fluid, Oct. 1953.

[42]

W. H. Coulter. High speed blood cell counter and cell size analyzer. In Proceedings of the National Electronics Conference, volume 12, pages 1034–1042. National Electronics Conference, Inc., 1957. Hotel Sherman, Chicago (IL), 1956, Oct. 3.

[43]

M. Dao, C. T. Lim, and S. Suresh. Mechanics of the human red blood cell deformed by optical tweezers. J. Mech. Phys. Sol., 51:2259–2280, 2003.

[44]

A. Di Biasio and C. Cametti. Effect of the shape of human erythrocytes on the evaluation of the passive electrical properties of the cell membrane. Bioelectrochem., 65:163–169, 2005.

[45]

A. Di Biasio and C. Cametti. Effect of shape on the dielectric properties of biological cell suspensions. Bioelectrochem., 71:149–156, 2007.

[46]

P. Dimitrakopoulos. Effects of membrane hardness and scaling analysis for capsules in planar extensional flows. J. Fluid Mech., 745:487–508, 2014.

[47]

R. Dimova, N. Bezlyepkina, M. D. Jord¨o, R. L. Knorr, K. A. Riske, M. Staykova, P. M. Vlahovska, T. Yamamoto, P. Yang, and R. Lipowsky. Vesicles in electric fields: Some novel aspects of membrane behavior. Soft Mat., 5:3201–3212, 2009.

[48]

S. K. Doddi and P. Bagchi. Effect of inertia on the hydrodynamic interaction between two liquid capsules in simple shear flow. Int. J. Multiph. Flow, 34:375– 392, 2008.

[49]

S. K. Doddi and P. Bagchi. Lateral migration of a capsule in a plane Poiseuille flow in a channel. Int. J. Multiph. Flow, 34:966–986, 2008.

[50]

S. K. Doddi and P. Bagchi. Three-dimensional computational modeling of multiple deformable cells flowing in microvessels. Phys. Rev. E, 79(046318), 2009.

[51]

W. R. Dodson III and P. Dimitrakopoulos. Dynamics of strain-hardening and strain-softening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes. J. Fluid Mech., 641:263– 296, 2009.

[52]

W. R. Dodson III and P. Dimitrakopoulos. Tank-treading of erythrocytes in strong shear flows via a nonstiff cytoskeleton-based continuum computational modeling. Biophys. J., 99:2906–2916, 2010.

BIBLIOGRAPHY

159

[53]

P. N. Dooley and N. J. Quinlan. Effect of eddy length scale on mechanical loading of blood cells in turbulent flow. Ann. of Biomed. Eng., 37(12):2449–2458, 2009.

[54]

V. Doyeux, Y. Guyot, V. Chabannes, C. Prud’homme, and M. Ismail. Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics. J. Comput. Appl. Math., 246:251–259, 2013.

[55]

J. Dupire, M. Abkarian, and A. Viallat. Chaotic dynamics of red blood cells in a sinusoidal flow. Phys. Rev. Lett., 104(168101), 2010.

[56]

J. Dupire, M. Abkarian, and A. Viallat. A simple model to understand the effect of membrane shear elasticity and stress-free shape on the motion of red blood cells in shear flow. Soft Mat., in press, 2015.

[57]

J. Dupire, M/. Socol, and A. Viallat. Full dynamics of e red blood cell in shear flow. Proc. Natl Acad. Sc. USA, 109(51):20808–20813, 2012.

[58]

C. Dupont, A.-V. Salsac, B. Barth`es-Biesel, M. Vidrascu, and P. Le Tallec. Influence of bending resistance on the dynamics of a spherical capsule in shear flow. Phys. Fluids, 27(051902), 2015.

[59]

D. Eberly. Least squares fitting of data. Technical report, Geometric Tools, LLC, 2008.

[60]

C. D. Eggleton and A. S. Popel. Large deformation of red blood cell ghosts in a simple shear flow. Phys. Fluids, 10(8):1834–1845, 1998.

[61]

E. A. Evans and Y. C. Fung. Improved measurements of the erythrocyte geometry. Microv. Res., 4:335–347, 1972.

[62]

A. Farutin, T. Biben, and C. Misbah. 3D numerical simulations of vesicle and inextensible capsule dynamics. J. Comput. Phys., 275:539–568, 2014.

[63]

S. Faure, S. Martin, B. Maury, and T. Takahashi. Towards the simulation of dense suspensions: a numerical tool. ESAIM: Proc., 28:55–79, 2009.

[64]

D. A. Fedosov. Multiscale Modeling of Blood Flow and Soft Matter. PhD thesis, Brown University, 2010.

[65]

D. A. Fedosov, B. Caswell, and G. E. Karniadakis. A multiscale red blood cell model with accurate mechanics, rheology, and dynamics. Biophys. J., 98:2215– 2225, 2010.

[66]

D. A. Fedosov, B. Caswell, A. S. Popel, and G. E. Karniadakis. Blood flow and cell-free layer in microvessels. Microcirc., 17:615–628, 2010.

160

BIBLIOGRAPHY

[67]

D. A. Fedosov, M. Dao, G. E. Karniadakis, and S. Suresh. Computational biorheology of human blood flow in health and disease. Ann. of Biomed. Eng., 42(2):368–387, 2013.

[68]

D. A. Fedosov, W. Pan, B. Caswell, G. Gompper, and G. E. Karniadakis. Predicting human blood viscosity in silico. Proc. Natl Acad. Sc. USA, 108(29):11772– 11777, 2011.

[69]

D. A. Fedosov, M. Peltom¨ aki, and G. Gompper. Deformation and dynamics of red blood cells in flow through cylindrical microchannels. Soft Mat., 10:4258–4267, 2014.

[70]

T. M. Fischer. On the energy dissipation in a tank-treading human red blood cell. Biophys. J., 32:863–868, 1980.

[71]

T. M. Fischer. Shape memory of human red blood cells. Biophys. J., 86:3304– 3313, 2004.

[72]

T. M. Fischer. Tank-tread frequency of the red cell membrane: Dependence on the viscosity of the suspending medium. Biophys. J., 93:2553–2561, 2007.

[73]

T. M. Fischer and R. Korzeniewski. Threshold shear stress for the transition between tumbling and tank-treading of red blood cells in shear flow: dependence on the viscosity of the suspending medium. J. Fluid Mech., 736:351–365, 2013.

[74]

T. M. Fischer, M. St¨ohr-Liesen, and H. Schmid-Sch¨onbein. The red cell as a fluid droplet: Tank tread-like motion of the human erythrocyte membrane in shear flow. Science, 202:894–896, 1978.

´ Foessel, J. Walter, A.-V. Salsac, and D. Barth`es-Biesel. Influence of internal [75] E. viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow. J. Fluid Mech., 672:477–486, 2011. [76]

J. B. Freund. The flow of red blood cells through a narrow spleen-like slit. Phys. Fluids, 25(110807), 2013.

[77]

J. B. Freund. Numerical simulation of flowing blood cells. Ann. Rev. Fluid Mech., 46:67–95, 2014.

[78]

J. B. Freund and M. M. Orescanin. Cellular flow in a small blood vessel. J. Fluid Mech., 671:466–490, 2011.

[79]

G. Fuhr and P. I. Kuzmin. Behavior of cells in rotating electric fields with account to surface charges and cell structures. Biophys. J., 50:789–795, 1986.

[80]

Y. C. Fung. Biomechanics – Mechanical properties of living tissues. SpringerVerlag, 2nd edition, 1993.

BIBLIOGRAPHY

161

[81]

Y. C. Fung. Biomechanics – Circulation. Springer-Verlag, 2nd edition, 1997.

[82]

R. V. Garimella and B. K. Swartz. Curvature estimation for unstructured triangulations of surfaces. LANL Tech. Rep., 2013.

[83]

G. Ghigliotti, T. Biben, and C. Misbah. Rheology of a dilute two-dimensional suspension of vesicles. J. Fluid Mech., 653:489–518, 2010.

[84]

H. L. Goldsmith and J. Marlow. Flow behaviour of erythrocytes. I. rotation and deformation in dilute suspensions. Proc. R. Soc. Lond. B, 182:351–384, 1972.

[85]

M. Gross, T. Kr¨ uger, and F. Varnik. Rheology of dense suspensions of elastic capsules: normal stresses, yield stress, jamming and confinement effects. Soft Mat., 10:4360–4372, 2014.

[86]

J. Guo, T. S. Pui, A. R. A. Rahman, and Y. Kang. 3D numerical simulation of a Coulter counter array with analysis of electrokinetic forces. Electrophoresis, 34:417–424, 2013.

[87]

M. R. Hardeman, P. T. Goedhart, J. G. G. Dobbe, and K. P. Lettinga. Laserassisted optical rotational cell analyser (l.o.r.c.a.); i. a new instrument for measurement of various structural hemorheological parameters. Clin. Hemorheo., 14(4):605–618, 1994.

[88]

Y. Hayashi, I. Oshige, Y. Katsumoto, S. Omori, A. Yasuda, and K. Asami. Dielectric inspection of erythrocyte morphology. Phys. Med. Biol., 53:2553–2564, 2008.

[89]

W. Helfrich. Elastic properties of lipid bilayers: Theory and possible experiments. Z. Naturforsch, 28 c:693–703, 1973.

[90]

S. H´enon. A new determination of the shear modulus of the human erythrocyte membrane using optical tweezers. Biophys. J., 76:1145–1151, 1999.

[91]

J. M. Higgins. Red blood cell population dynamics. Clin. Lab. Med., 35:43–57, 2015.

[92]

B. D. Horne, J. L. Anderson, J. M. John, A. Weaver, T. L. Bair, K. Jensen, D. G. Renlund, and J. B. Muhlestein. Which white blood cell subtypes predict increased cardiovascular risk? J. Am. Coll. Cardio., 45:1638–1643, 2005.

[93]

X.-Q. Hu, A.-V. Salsac, and D. Barth`es-Biesel. Flow of a spherical capsule in a pore with circular or square cross-section. J. Fluid Mech., 705:176–194, 2012.

[94]

X.-Q. Hu, B. S´ev´eni´e, A.-V. Salsac, E. Leclerc, and B. Barth`es-Biesel. Characterizing the membrane properties of capsules flowing in a square-section microfluidic channel: Effects of the membrane constitutive law. Phys. Rev. E, 87:063008, 2013.

162

BIBLIOGRAPHY

[95]

R. L. Iman. Latin Hypercube Sampling. Encyclopedia of Quantitative Risk Analysis and Assessment. Wiley, 2008.

[96]

D. Is`ebe and P. N´erin. Numerical simulation of particle dynamics in an orificeelectrode system. Application to counting and sizing by impedance measurement. Int. J. Numer. Meth. Biomed. Eng., 29(4):462–475, 2013.

[97]

M. Ismail and A. Lefebvre-Lepot. A‘necklace’ model for vesicles simulations in 2D. Int. J. Numer. Meth. Fluids, 76:835–854, 2014.

[98]

A. V. Jagtiani, J. Zhe, J. Hu, and J. Carletta. Detection and counting of microscale particles and pollen using a multi-aperture Coulter counter. Meas. Sci. Technol., 17:1706–1714, 2006.

[99]

G. B. Jeffery. The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. A, 102(715):161–179, 1922.

[100] J Jeˇzek, S. Saic, and K. Segeth. Numerical modeling of the movement of a rigid particle in viscous fluid. App. Math., 44:469–479, 1999. [101] T. B. Jones. Basic theory of dielectrophoresis and electrorotation. Med. Biol. Mag., 6:33–42, 2003. [102] M. Junk and R. Illner. A new derivation of Jeffery’s equation. Mech., 9:455–488, 2007.

IEEE Eng.

J. Math. Fluid

[103] H. Kamada, Y. Imai, M. Nakamura, T. Ishikawa, and T. Yamaguchi. Computational analysis on the mechanical interaction between a thrombus and red blood cells: Possible causes of membrane damage of red blood cells at microvessels. Med. Eng. Phys., 34(10):1411–1420, 2012. [104] S. R. Keller and R. Skalak. Motion of a tank-treading ellipsoidal particle in a shear flow. J. Fluid Mech., 120:27–47, 1982. [105] S. Kessler, R. Finker, and U. Seifert. Swinging and tumbling of elastic capsules in shear flow. J. Fluid Mech., 605:207–226, 2008. [106] K. Khairy, J. Foo, and J. Howard. Shapes of red blood cells: Comparison of 3D confocal images with the bilayer-couple model. Cell. Mol. Bioeng., 1(2-3):173– 181, 2008. [107] T. Kl¨ oppel and W. A. Wall. A novel two-layer, coupled finite element approach for modeling the nonlinear elastic and viscoelastic behavior of human erythrocytes. Biomech. Model. Mechanobiol., 10:445–459, 2011. [108] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky. Fluid vesicles in shear flow. Phys. Rev. Lett., 77(17):3685–3688, 1996.

BIBLIOGRAPHY

163

[109] T. Kr¨ uger, B. Kaoui, and J. Harting. Interplay of inertia and deformability on rheological properties of a suspension of capsules. J. Fluid Mech., 751:725–745, 2014. [110] E. Lac and D. Barth`es-Biesel. Deformation of a capsule in simple shear flow: Effect of membrane prestress. Phys. Fluids, 17(072105), 2005. [111] E. Lac, D. Barth`es-Biesel, N. A. Pelekasis, and J. Tsamopoulos. Spherical capsules in three-dimensional unbounded Stokes flows: effect of the membrane constitutive law and onset of buckling. J. Fluid Mech., 516:303–334, 2004. [112] E. Lac, A. Morel, and D. Barth`es-Biesel. Hydrodynamic interaction between two identical capsules in simple shear flow. J. Fluid Mech., 573:149–169, 2007. [113] G.-B. Lee, C.-C. Chang, S-B. Huang, and R.-J. Yang. The hydrodynamic focusing effect inside rectangular microchannels. J. Micromech. Microeng., 16:1024–1032, 2006. [114] G. Lenormand, S. H´enon, A. Richet, J. Sim´eon, and F. Gallet. Direct measurement of the area expansion and shear moduli of the human red blood cell membrane skeleton. Biophys. J., 81:43–56, 2001. [115] X. Li and K. Sarkar. Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane. J. Comput. Phys., 227:4998–5018, 2008. [116] X. Z. Li, D. Barth`es-Biesel, and A. Helmy. Large deformations and burst of a capsule freely suspended in an elongational flow. J. Fluid Mech., 187:179–196, 1988. [117] Y. Li, E. Jung, W. Lee, H. G. Lee, and J. Kim. Volume preserving immersed boundary methods for two-phase fluid flows. Int. J. Numer. Meth. Fluids, 69:842– 858, 2012. [118] G. H. W. Lim, M. Wortiz, and R. Mukhopadhyay. Stomatocyte-discocyteechinocyte sequence of the human red blood cell: Evidence for the bilayer-couple hypothesis from membrane mechanics. Proc. Natl Acad. Sc. USA, 99(26):16766– 16769, 2002. [119] G. H. W. Lim, M. Wortiz, and R. Mukhopadhyay. Red Blood Cell Shapes and Shape Transformations: Newtonian Mechanics of a Composite Membrane, volume Lipid Bilayers and Red Blood Cells of Soft Matter, chapter 2, pages 94–269. WILEY-VCH Verlag GmbH & Co. KGaA, 2008. [120] W. K. Liu, S. Jun, and Y. F. Zhang. Reproducing kernel particle methods. Int. J. Numer. Meth. Fluids, 20:1081–1106, 1995.

164

BIBLIOGRAPHY

[121] V. Lleras. Mod´elisation, analyse et simulation de probl`emes de contact en m´ecanique des solides et des fluides. PhD thesis, Universit´e de Franche Comt´e – Besan¸con, 2009. [122] J. M. L´opez-Herrera, A. Popinet, and M.A. Herrada. A charge-conservative approach for simulating electrohydrodynamic two-phase flows using volume-of-fluid. J. Comput. Phys., 230:1939—1955, 2011. [123] Z. Y. Luo, L. He, and B. F. Bai. Deformation of spherical compound capsules in simple shear flow. J. Fluid Mech., 775:77–104, 2015. [124] Z. Y. Luo, S. Q. Wang, L. He, F. Xu, and B. F. Bai. Inertia-dependent dynamics of three-dimensional vesicles and red blood cells in shear flow. Soft Mat., 9:9651– 9660, 2013. [125] M. MacMeccan, J. R. Clausen, G. P. Neitzel, and C. K. Aidun. Simulating deformable particle suspensions using a coupled lattice-Boltzmann and finiteelement method. J. Fluid Mech., 618:13–39, 2009. [126] E. Maitre. Equations de transport, level set et m´ecanique eul´erienne. application au couplage fluide-structure. Technical report, HDR report, Universit´e de Grenoble, 2008. [127] E. Maitre, T. Milcent, G.-H. Cottet, A. Raoult, and Y. Usson. Applications of level set methods in computational biophysics. Math. Comput. Mod., 49(11– 12):2161–2169, 2009. [128] M. Malandin, N. Maheu, and V. Moureau. Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines. J. Comput. Phys., 238:32–47, 2013. [129] J. L. McWhirter, H. Noguchi, and G. Gompper. Flow-induced clustering and alignment of vesicles and red blood cells in microcapillaries. Proc. Natl Acad. Sc. USA, 106(15):6039–6043, 2009. [130] J. L. McWhirter, H. Noguchi, and G. Gompper. Deformation and clustering of red blood cells in microcapillary flows. Soft Mat., 7:10967–10977, 2011. [131] S. Mendez, C. Chnafa, E. Gibaud, J. Sig¨ uenza, V. Moureau, and F. Nicoud. YALES2BIO: a computational fluid dynamics software dedicated to the prediction of blood flows in biomedical devices. In 5th International Conference on Biomedical Engineering in Vietnam, IFMBE Proceedings Series, 2014. [132] S. Mendez, E. Gibaud, and F. Nicoud. An unstructured solver for simulations of deformable particles in flows at arbitrary Reynolds numbers. J. Comput. Phys., 256(1):465–483, 2014.

BIBLIOGRAPHY

165

[133] T. Milcent. Une approche eul´erienne du couplage fluide-structure, analyse math´ematique et applications en biom´ecanique. PhD thesis, Universit´e de Grenoble. Joseph Fourier, 2009. [134] J. P. Mills, L. Qie, M. Dao, C. T. Lim, and S. Suresh. Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers. Mech. Chem. Biosys., 1(3):169–180, 2004. [135] C. Misbah. Vacillating breathing and tumbling of vesicles under shear flow. Phys. Rev. Lett., 96(028104), 2006. [136] N. Mohandas and E. A. Evans. Mechanical properties of the red cell membrane in relation to molecular structure and genetic defects. Ann. Rev. Biophys. Biomol. Struct., 23:787–818, 1994. [137] V. Moureau, P. Domingo, and L. Vervisch. Design of a massively parallel CFD code for complex geometries. Comp. Rend. M´ec., 339(2-3):141–148, 2011. [138] V. Narsimhan, A. Spann, and E. S. G. Shaqfeh. The mechanism of shape instability for a vesicle in extensional flow. J. Fluid Mech., 750:144–190, 2014. [139] V. Narsimhan, H. Zhao, and E. S. G. Shaqfeh. Coarse-grained theory to predict the concentration distribution of red blood cells in wall-bounded Couette flow at zero Reynolds number. Phys. Fluids, 25(061901), 2013. [140] H. Noguchi and G. Gompper. Fluid vesicles with viscous membranes in shear flow. Phys. Rev. Lett., 93(258102), 2004. [141] H. Noguchi and G. Gompper. Dynamics of fluid vesicles in shear flow: Effect of membrane viscosity and thermal fluctuations. Phys. Rev. E, 72(011901), 2005. [142] H. Noguchi and G. Gompper. Shape transitions of fluid vesicles and red blood cells in capillary flows. Proc. Natl Acad. Sc. USA, 102(40):14159–14164, 2005. [143] H. Noguchi and G. Gompper. Swinging and tumbling of fluid vesicles in shear flow. Phys. Rev. Lett., 98:128103, 2007. [144] W. Pan, B. Caswell, and G. E. Karniadakis. A low-dimensional model for the red blood cell. Soft Mat., 6:4366–4376, 2010. [145] Z. Peng, R. J. Asaro, and Q. Zhu. Multiscale simulation of erythrocyte membranes. Phys. Rev. E, 81(031904), 2010. [146] Z. Peng, R. J. Asaro, and Q. Zhu. Multiscale modelling of erythrocytes in Stokes flow. J. Fluid Mech., 686:299–337, 2011.

166

BIBLIOGRAPHY

[147] Z. Peng, A. Mashayekh, and Q. Zhu. Erythrocyte responses in low-shear-rate flows: effects of non-biconcave stress-free state in the cytoskeleton. J. Fluid Mech., 742:96–118, 2014. [148] C. S. Peskin. The immersed boundary method. Acta Num., 11:479–517, 2002. ˇ s. The prolate-to-oblate shape transition of [149] P. Peterlin, S. Svetina, and B. Zekˇ phospholipid vesicles in response to frequency variation of an AC electric field can be explained by the dielectric anisotropy of a phospholipid bilayer. J. Phys.: Condens. Matter, 19:136220, 2007. [150] A. Pinelli, I. Z. Naqavi, U. Piomelli, and J. Favier. Immersed-boundary methods for general finite-difference and finite-volume Navier–Stokes solvers. J. Comput. Phys., 229:9073–9091, 2010. [151] I. V. Pivkin and G. E. Karniadakis. Accurate coarse-grained modeling of red blood cells. Phys. Rev. Lett., 101(118105), 2008. [152] C. Pozrikidis. The axisymmetric deformation of a red blood cell in uniaxial straining Stokes flow. J. Fluid Mech., 216:231–254, 1990. [153] C. Pozrikidis. Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press., 1992. [154] C. Pozrikidis. Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow. J. Fluid Mech., 297:123–152, 1995. [155] C. Pozrikidis. Effect of membrane bending stiffness on the deformation of capsules in simple shear flow. J. Fluid Mech., 440:269–291, 2001. [156] C. Pozrikidis. Modeling and Simulation of Capsules and Biological Cells. Boca Raton: Chapman & Hall/CRC, 2003. [157] C. Pozrikidis. Numerical simulation of the flow-induced deformation of red blood cells. Ann. of Biomed. Eng., 31:1194–1205, 2003. [158] C. Pozrikidis. Axisymmetric motion of a file of red blood cells through capillaries. Phys. Fluids, 17(031503), 2005. [159] A. R. Pries, T. W. Secomb, and P. Gaehtgens. Biophysical aspects of blood flow in the microvasculature. Cardiov. Res., 32:654–667, 1996. [160] Z. Qin, J. Zhe, and G.-X. Wang. Effects of particle’s off-axis position, shape, orientation and entry position on resistance changes of micro Coulter counting devices. Meas. Sci. Technol., 22:1–10, 2011. [161] C. Qu´eguiner and D. Barth`es-Biesel. Axisymmetric motion of capsules through cylindrical channels. J. Fluid Mech., 348:349–376, 1997.

BIBLIOGRAPHY

167

[162] S. Ramanujan and C. Pozrikidis. Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities. J. Fluid Mech., 361:117–143, 1998. [163] D. A. Reasor, J. R. Clausen, and C. K. Aidun. Coupling the lattice-Boltzmann and spectrin-link methods for the direct numerical simulation of cellular blood flow. Int. J. Numer. Meth. Fluids, 68:767–781, 2012. [164] A. M. Roma, C. S. Peskin, and M. J. Berger. An adaptive version of the immersed boundary method. J. Comput. Phys., 153:509–534, 1999. [165] R. Roscoe. On the rheology of a suspension of viscoelastic spheres in a viscous liquid. J. Fluid Mech., 28:273–293, 1967. [166] U. Seifert. Morphology and dynamics of vesicles. Current Opinion in Colloid & Interface Science, 1(3):350–357, 1996. [167] U. Seifert. Configurations of fluid membranes and vesicles. Adv. Phys., 46(1):13– 137, 1997. [168] U. Seifert. Fluid membranes in hydrodynamic flow fields: Formalism and an application to fluctuating quasispherical vesicles in shear flow. Eur. Phys. J. B, 8:405–415, 1999. [169] U. Seifert, K. Berndl, and R. Lipowsky. Shape transformations of vesicles: Phase diagram for spontaneous-curvature and bilayer-coupling models. Phys. Rev. A, 44(2):1182–1202, 1991. [170] U. Seifert and R. Lipowsky. Morphology of Vesicles: ch.8 of Handbook of Biological Physics, volume 1. Elsevier Science B. V., 1995. [171] L. Shi, T.-W. Pan, and R. Glowinski. Lateral migration and equilibrium shape and position of a single red blood cell in bounded Poiseuille flows. Phys. Rev. E, 86(056308), 2012. [172] J. Sig¨ uenza. Numerical simulation of the red blood cells. Master’s thesis, Polytech’Montpellier, 2013. [173] J. Sig¨ uenza, S. Mendez, and F. Nicoud. Characterisation of a dedicated mechanical model for red blood cells: numerical simulations of optical tweezers experiment. Comput. Meth. Biomech. Biomed. Eng., 17(supp. 1):28–29, 2014. [174] R. Skalak, A. Tozeren, R. P. Zarda, and S. Chien. Strain energy function of red blood cell membranes. Biophys. J., 13:245–264, 1973. [175] J. M. Skotheim and T. W. Secomb. Red blood cells and other nonspherical capsules in shear flow: Oscillatory dynamics and the tank-treading-to-tumbling transition. Phys. Rev. Lett., 98(078301), 2007.

168

BIBLIOGRAPHY

[176] A. J. Smolen, L. L. Wright, and T. J. Cunningham. Neuron numbers in the superior cervical sympathetic ganglion of the rat: a critical comparison of methods for cell counting. J. Neurocyto., 12:739–750, 1983. [177] R. Specogna and F. Trevisan. A discrete geometric approach to solving time independent schr¨odinger equation. J. Comput. Phys., 230:1370–1381, 2011. [178] D. Spencer and H. Morgan. Positional dependence of particles in microfludic impedance cytometry. Lab. Chip, 11:1234–1239, 2011. [179] Y. Sui, X. B. Chen, Y. T. Chew, P. Roy, and H. T. Low. Numerical simulation of capsule deformation in simple shear flow. Comput. Fluids, 39:242–250, 2010. [180] Y. Sui, Y. T. Chew, P. Roy, Y. P. Cheng, and H. T. Low. Dynamic motion of red blood cells in simple shear flow. Phys. Fluids, 20(112106), 2008. [181] Y. Sui, H. T. Low, Y. T. Chew, and P. Roy. Tank-treading, swinging, and tumbling of liquid-filled elastic capsules in shear flow. Phys. Rev. E, 77(016310), 2008. [182] T. Sun and H. Morgan. Single-cell microfluidic impedance cytometry: a review. Microfluid. Nanofluid., 8:423–443, 2010. [183] S. Suresh, J. Spatz, J. P. Mills, A. Micoulet, M. Dao, C. T. Lim, M. Beil, and T. Seufferlein. Connections between single-cell biomechanics and human disease states: gastrointestinal cancer and malaria. Acta Biomaterialia, 1:15–30, 2005. [184] G. Tomaiuolo. Biomechanical properties of red blood cells in health and disease towards microfluidics. Biomicrofluid., 8(5):1–20, 2014. [185] G. Tomar, D. Gerlach, G. Biswas, N. Alleborn, A. Sharma, F. Durst, S. W. J. Welch, and A. Delgado. Two-phase electrohydrodynamic simulations using a volume-of-fluid approach. J. Comput. Phys., 227:1267—1285, 2007. [186] L. A. Torre, F. Bray, R. L. Siegel, J. Ferlay, J. Lortet-Tieulent, and A. Jemal. Global cancer statistics, 2012. CA Cancer J. Clin., 108:65–87, 2015. [187] R. Tran-Son-Tay, S. P. Sutera, and P. R. Rao. Determination of red blood cell membrane viscosity from rheoscopic observations of tank-treading motion. Biophys. J., 46:65–72, 1984. [188] R. Trozzo, G. Boedec, M. Leonetti, and M. Jaeger. Axisymmetric boundary element method for vesicles inacapillary. J. Comput. Phys., 289:62–82, 2015. [189] K.-I. Tsubota, S. Wada, and H. Liu. Elastic behavior of a red blood cell with the membrane’s nonuniform natural state: equilibrium shape, motion transition under shear flow, and elongation during tank-treading motion. Biomech. Model. Mechanobiol., 13(4):735–746, 2014.

BIBLIOGRAPHY

169

[190] S. O. Unverdi and G. Tryggvason. A front-tracking method for viscous, incompressible, multi-fluid flows. J. Comput. Phys., 100:25–37, 1992. [191] A. Valero, T. Braschler, and P Renaud. A unified approach to dielectric single cell analysis: Impedance and dielectrophoretic force spectroscopy. Lab. Chip, 10:2216–2225, 2010. [192] B. P. Van Poppel, O. Desjardins, and J. W. Daily. A ghost fluid, level set methodology for simulating multiphase electrohydrodynamic flows with application to liquid fuel injection. J. Comput. Phys., 229:7977–7996, 2010. [193] S. K. Veerapaneni, D. Gueyffier, G. Biros, and D. Zorin. A numerical method for simulating the dynamics of 3d axisymmetric vesicles suspended in viscous flows. J. Comput. Phys., 228:7233–7249, 2009. [194] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros. A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D. J. Comput. Phys., 228:2334–2353, 2009. [195] S. K. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin. A fast algorithm for simulating vesicle flows in three dimensions. J. Comput. Phys., 230:5610–5634, 2011. [196] O. Vizika and D. A. Saville. The electrohydrodynamic deformation of drops suspended in liquids in steady and oscillatory electric field. J. Fluid Mech., 239:1–21, 1992. [197] P. M. Vlahovska, R. Serral Graci` a, S. Aranda-Espinoza, and R. Dimova. Electrohydrodynamic model of vesicle deformation in alternating electric fields. Biophys. J., 96:4789–4803, 2009. [198] P. M. Vlahovska, Y.-N. Young, G. Danker, and C. Misbah. Dynamics of a nonspherical microcapsule with incompressible interface in shear flow. J. Fluid Mech., 678:221–247, 2011. [199] B. Vozarova, N. Stefan, R. S. Lindsay, A. Saremi, R. E. Pratley, C. Bogardus, and P. A. Tataranni. High alanine aminotransferase is associated with decreased hepatic insulin sensitivity and predicts the development of type 2 diabetes. Diabetes, 51:1889–1895, 2002. [200] J. Walter, A.-V. Salsac, and D. Barth`es-Biesel. Ellipsoidal capsules in simple shear flow: prolate versus oblate initial shapes. J. Fluid Mech., 676:318–347, 2011. [201] J. Walter, A.-V. Salsac, D. Barth`es-Biesel, and P. Le Tallec. Coupling of finite element and boundary integral methods for a capsule in a Stokes flow. Int. J. Numer. Meth. Eng., 83:829–850, 2010.

170

BIBLIOGRAPHY

[202] X. Wang and W. K. Liu. Extended immersed boundary method using FEM and RKPM. Comput. Meth. Appl. Mech. Eng., 193:1305–1321, 2004. [203] J. C. Weaver and Y. A. Chizmadzhev. Theory of electroporation: A review. Bioelectro. Bioenerg., 41:135–160, 1996. [204] J. Wu and C. K. Aidun. A numerical study of the effect of fibre stiffness on the rheology of sheared flexible fibre suspensions. J. Fluid Mech., 662:123–133, 2010. [205] A. Z. K. Yazdani and P. Bagchi. Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method. Phys. Rev. E, 85(056308), 2012. [206] A. Z. K. Yazdani and P. Bagchi. Influence of membrane viscosity on capsule dynamics in shear flow. J. Fluid Mech., 718:569–595, 2013. [207] A. Z. K. Yazdani, R. M. Kalluri, and P. Bagchi. Tank-treading and tumbling frequencies of capsules and red blood cells. Phys. Rev. E, 83(046305), 2011. [208] J. Zhang, P. C. Johnson, and A. S. Popel. Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method. J. Biomech., 41:47–55, 2008. [209] J. Zhang, P. C. Johnson, and A. S. Popel. Effects of erythrocyte deformability and aggregation on the cell free layer and apparent viscosity of microscopic blood flows. Microv. Res., 77:265–272, 2009. [210] H. Zhao, A. H. G. Isfahani, L. N. Olson, and J. B. Freund. A spectral boundary integral method for flowing blood cells. J. Comput. Phys., 229:3726–3744, 2010. [211] H. Zhao and E. S. G. Shaqfeh. The dynamics of a vesicle in simple shear flow. J. Fluid Mech., 674:578–604, 2011. [212] H. Zhao and E. S. G. Shaqfeh. Shear-induced platelet margination in a microchannel. Phys. Rev. E, 83(061924), 2011. [213] H. Zhao and E. S. G. Shaqfeh. The dynamics of a non-dilute vesicle suspension in a simple shear flow. J. Fluid Mech., 725:709–731, 2013. [214] H. Zhao and E. S. G. Shaqfeh. The shape stability of a lipid vesicle in a uniaxial extensional flow. J. Fluid Mech., 719:345–361, 2013. [215] H. Zhao, E. S. G. Shaqfeh, and V. Narsimhan. Shear-induced particle migration and margination in a cellular suspension. Phys. Fluids, 24(011902):1–21, 2012. [216] H. Zhao, A. Spann, and E. S. G. Shaqfeh. The dynamics of a vesicle in a wallbound shear flow. Phys. Fluids, 23(121901), 2011.

BIBLIOGRAPHY

171

[217] M. Zhao and P. Bagchi. Dynamics of microcapsules in oscillating shear flow. Phys. Fluids, 23(111901), 2011. [218] O. Y. Zhong-can and W. Helfrich. Bending energy of vesicle membranes: General expressions for the first, second, and third variation of the shape energy and applications to spheres and cylinders. Phys. Rev. A, 39(10):5280–5288, 1989. [219] L. Zhu, C. Rorai, D. Mitra, and L. Brandt. A microfluidic device to sort capsules by deformability: a numerical study. Soft Mat., 10:7705–7713, 2014.