УНИВЕРЗИТЕТ У КРАГУЈЕВЦУ ФАКУЛТЕТ ИНЖЕЊЕРСКИХ НАУКА МОДЕЛИРАЊЕ КРЕТАЊА ДЕФОРМАБИЛНОГ ТЕЛА У ФЛУИДУ И ПРИМЕНА У БИОМЕДИЦИНСКОМ ИНЖЕЊЕРИНГУ Докторска дисертација Ментор Кандидат др Ненад Филиповић, Тијана Ђукић, редовни професор маст. инж. маш. Крагујевац, 2014. I. Аутор Име и презиме: Тијана Ђукић Датум и место рођења: 01.04.1988. Крагујевац Садашње запослење: Истраживачко-развојни центар за биоинжењеринг, Крагујевац II. Докторска дисертација Наслов: Моделирање кретања деформабилног тела у флуиду и примена у биомедицинском инжењерингу Број страница: 168 Број слика: 100 Број библиографских података: 114 Установа и место где је рад израђен: Факултет инжењерских наука, Крагујевац Научна област (УДК): Примењена информатика и рачунарско инжењерство (519.87) Ментор: др Ненад Филиповић, редовни професор III. Оцена и одбрана Датум пријаве теме: 04.02.2014. Број одлуке и датум прихватања докторске дисертације: Комисија за оцену подобности теме и кандидата: 1. др Ненад Филиповић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 2. др Вељко Милутиновић, ред. проф., Електротехнички факултет у Београду, Научне области: Рачунарска техника и информатика 3. др Радован Славковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 4. др Мирослав Живковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 5. др Гордана Јовичић, ванр. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 6. др Милован Матовић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Нуклеарна медицина, Медицинска информатика 7. др Мирко Росић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Физиологија Комисија за оцену докторске дисертације: 1. др Ненад Филиповић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 2. др Вељко Милутиновић, ред. проф., Електротехнички факултет у Београду, Научне области: Рачунарска техника и информатика 3. др Радован Славковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 4. др Мирослав Живковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 5. др Гордана Јовичић, ванр. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 6. др Милован Матовић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Нуклеарна медицина, Медицинска информатика 7. др Мирко Росић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Физиологија Комисија за одбрану докторске дисертације: 1. др Ненад Филиповић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 2. др Вељко Милутиновић, ред. проф., Електротехнички факултет у Београду, Научне области: Рачунарска техника и информатика 3. др Радован Славковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 4. др Мирослав Живковић, ред. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 5. др Гордана Јовичић, ванр. проф., Факултет инжењерских наука у Крагујевцу, Научне области: Примењена механика, Примењена информатика и рачунарско инжењерство 6. др Милован Матовић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Нуклеарна медицина, Медицинска информатика 7. др Мирко Росић, ред. проф., Факултет медицинских наука у Крагујевцу, Научне области: Физиологија Датум одбране докторске дисертације: Резиме рада Кардиоваскуларни систем је један од најважнијих система у људском организму. Истраживања струјања крви су клинички веома значајна због могућности побољшане дијагностике болести, планирања одговарајућих превентивних мера, анализе транспорта лекова итд. Ове појаве се тешко експериментално испитују и због тога у овој области нумеричке симулације могу да допринесу добијању много нових и корисних информација. У оквиру овог рада представљен је комплетан нумерички модел који успешно симулира тродимензионално струјање крви кроз крвне судове у људском организму. Ово је показано на неколико примера, где је моделирано струјање крви кроз људску аорту, каротидну и коронарну артерију. Такође, представљен је нумерички модел који симулира кретање црвених крвних зрнаца и сферних честица кроз комплексне геометријске домене, као што су мали крвни судови са стенозом, бифуркацијом, као и кроз микрофлуидне чипове за сепарацију канцер ћелија, али и кроз реалне домене добијене експерименталним снимањем протока крви кроз посебну врсту рибе, тзв. зебра рибу. С обзиром да један овакав софтвер за научне симулације захтева много рачунарских ресурса и да извршавање једне симулације траје доста дуго, приликом имплементације нумеричке методе посебна пажња посвећена је техникама паралелизације. Компјутерски програм је развијан тако да поред рада на стандардним PC рачунарима, може да се извршава и на компјутерским кластерима са великим бројем процесора, као што су Tesla суперкомпјутери, које производи светска компанија NVIDIA. Постигнута су значајна убрзања развијеног софтвера, тако да се нумеричке симулације могу извршавати за само неколико минута, наспрам пар сати колико је потребно за извршавање исте симулације на стандардном PC рачунару, што омогућава практично интерактивно праћење кретања деформабилних тела у реалном времену. Због нумеричке методе која доста реално осликава услове у људском организму, велике тачности која је показана кроз поређење са експерименталним резултатима и бројним другим решењима из литературе, као и због брзог извршавања на компјутерским кластерима, које је омогућено паралелизацијом, нумеричка метода представљена у овом раду и развијени софтвер имају велики потенцијал када је реч о практичној примени у симулацијама реалних проблема у биомедицинском инжењерингу. Кључне речи: нумерички модел, компјутерска симулација, паралелизација, капиларни крвни судови, комплексно струјање крви Abstract Cardiovascular system is one of the most important systems in human organism. Investigation of blood flow is clinically very relevant due to the possibility of improved diagnostic of diseases, planning of appropriate preventive measures, analisys of drug transport etc. Experimental investigation of these phenomena is difficult and therefore numerical simulation can contribute to the acquisition of manu new and useful information. Within this dissertation a complete numerical model is presented, that is capable of simulating three-dimensional blood flow through the blood vessels in human organism. This is demonstrated on several examples, where blood flow through human aorta, carotid and coronary arteries is modelled. Also, numerical model that simulates solid-fluid interaction is presented and this model is used to simulate motion of red blood cells and spherical particles through complex geometric domains, such as small blood vessels with stenosis and bifurcation, as well as microfluidic chips for cancer cell separation and real domains obtained by experimental recording of blood flow through a specific type of fish, so called zebrafish. This type of software for numerical simulation requires a lot of computer resources and the execution of one simulation lasts quite long. Therefore, during the implementation of the numerical method, special attention was dedicated to parallelisation techniques. The software was developed such that beside execution on standarad PC computer, it is possible to execute the program on computer clusters with large number of processors, such as Tesla supercomputers, that are manufactured by the world known company NVIDIA. Significant speed-ups were obtained and numerical simulation can now be executed in several minutes, instead of several hours, which was the case when the same simulation was executed on standard PC computer. This enables interactive tracking of motion of deformable bodies in real time. Due to the numerical method that describes the conditions in human organism quite realistically, big accuracy of the solutions that was demonstrated through the comparison with experimental results and other results presented in literature and fast program execution on computer clusters that is enabled with parallelisation, the numerical model that is presented in this dissertation and the developed software have a great potential when it comes to practical application in modeling real phenomena in biomedical engineering. Key words: numerical model, computer simulation, parallelisation, capillary blood vessels, complex blood flow Тијана Ђукић I Садржај Садржај ................................................................................................................................ I Списак слика .................................................................................................................... IV Списак табела .................................................................................................................... X Списак коришћених ознака .......................................................................................... XI Коришћене скраћенице ............................................................................................. XVI 1 Увод .................................................................................................................................. 1 2 Моделирање струјања флуида .................................................................................... 5 2.1 Boltzmann-ова једначина ............................................................................................ 6 2.2 Извођење Навије-Стоксове једначине ...................................................................... 7 2.3 Дискретизација Boltzmann-ове једначине ............................................................... 18 2.3.1 Дискретизација поља брзина ............................................................................ 18 2.3.2 Временска и просторна дискретизација .......................................................... 22 2.4 Израчунавање макроскопских величина коришћењем функције расподеле ...... 25 3 Моделирање понашања деформабилног тела ....................................................... 28 3.1 Моделирање реакције на површинско напрезање мембране ............................... 28 3.2 Моделирање реакције на промену запремине ....................................................... 42 3.3 Моделирање реакције на промену површине спољашње мембране ................... 43 3.4 Моделирање реакције на савијање мембране ........................................................ 44 4 Моделирање солид-флуид интеракције .................................................................. 46 4.1 Immersed boundary метода ....................................................................................... 47 5 Моделирање струјања вишекомпонентног флуида ............................................. 52 5.1 Моделирање струјања вишекомпонентног флуида применом Shan-Chen-овог модела ................................................................................................................................ 52 Докторска дисертација Садржај Тијана Ђукић II 5.2 Моделирање струјања више флуида различите вискозности раздвојених мембраном солида ............................................................................................................ 54 6 Паралелизација развијеног софтвера ..................................................................... 57 6.1 CUDA архитектура за програмирање GPU уређаја ............................................... 57 6.2 Паралелизација софтвера на бази lattice Boltzmann методе.................................. 59 6.3 Паралелизација на кластерима са више GPU уређаја ........................................... 63 6.4 Поређење брзине рада секвенцијалне и паралелизоване верзије LB солвера .... 67 6.4.1 Покретање паралелизоване верзије на једном GPU уређају ......................... 67 6.4.2 Покретање паралелизоване верзије на више GPU уређаја, уз примену OpenMP приступа ......................................................................................................... 71 6.4.3 Покретање паралелизоване верзије на кластеру са великим бројем GPU уређаја, уз примену MPI приступа .............................................................................. 72 7 Резултати симулација ................................................................................................. 75 7.1 Конверзија бездимензионих величина у физичке величине ................................ 75 7.2 Симулације струјања крви кроз кардиоваскуларни систем ................................. 78 7.2.1 Моделирање струјања крви кроз људску аорту .............................................. 79 7.2.2 Моделирање струјања крви кроз каротидну бифуркацију ............................ 82 7.2.3 Моделирање струјања крви кроз коронарну артерију ................................... 85 7.3 Симулација мешања две течности .......................................................................... 88 7.4 Симулације струјања крви кроз капиларне крвне судове ..................................... 90 7.4.1 Моделирање понашања сферне честице под дејством смичућег профила брзине ............................................................................................................................ 91 7.4.2 Моделирање понашања црвеног крвног зрнца под дејством смичућег профила брзине ............................................................................................................. 99 7.4.3 Моделирање кретања сферне честице у правом каналу .............................. 104 7.4.4 Моделирање кретања црвеног крвног зрнца у правом каналу .................... 106 7.4.5 Моделирање кретања сферне честице кроз артерију са сужењем .............. 108 7.4.6 Моделирање кретања црвеног крвног зрнца и канцер ћелије кроз микрофлуидни чип за сепарацију ............................................................................. 114 Докторска дисертација Садржај Тијана Ђукић III 8 Поређење решења нумеричких симулација са експерименталним резултатима ..................................................................................................................... 122 8.1 Кретање црвених крвних зрнаца кроз стаклену цев in-vitro ............................... 122 8.2 Кретање црвених крвних зрнаца кроз каротидну бифуркацију пацова ex-vivo 125 8.3 Моделирање протока крви кроз кардиоваскуларни систем зебра рибе ............ 129 9 Закључак ..................................................................................................................... 136 Додатак А ......................................................................................................................... 138 Додатак Б ......................................................................................................................... 140 Додатак Ц ......................................................................................................................... 142 Литература ...................................................................................................................... 144 Тијана Ђукић IV Списак слика Слика 3.1 - Скициране величине које су потребне за дефинисање сила реакције на површинско напрезање .......................................................................................................... 29 Слика 3.2 - Троугаони коначни елемент у глобалном Декартовом и природном координатном систему ........................................................................................................... 35 Слика 3.3 - Троугаони коначни елемент у глобалном и локалном Декартовом координатном систему ........................................................................................................... 38 Слика 3.4 - Скициране величине које су потребне за израчунавање запремине унутар мембране .................................................................................................................................. 42 Слика 3.5 - Скициране величине које су потребне за дефинисање сила реакције на промену површине ................................................................................................................. 44 Слика 3.6 - Скициране величине које су потребне за дефинисање сила реакције на савијање мембране ................................................................................................................. 45 Слика 4.1 – Размена информација приликом моделирања солид-флуид интеракције .... 47 Слика 4.2 - Међусобни утицај тачака домена флуида и солида ........................................ 48 Слика 4.3 – Дијаграм промене вредности Dirac-ове функције у зависности од растојања између чворова ........................................................................................................................ 49 Слика 4.4 - Алгоритам моделирања солид-флуид интеракције применом Immersed boundary методе ...................................................................................................................... 51 Слика 5.1 – Шематски приказ два флуида различите вискозности раздвојених мембраном солида .................................................................................................................. 55 Слика 6.1 - Илустрација мреже блокова и нити [43] .......................................................... 58 Слика 6.2 - Алгоритам рада LB солвера са означеним деловима које треба паралелизовати ....................................................................................................................... 60 Слика 6.3 – Алгоритам извршавања паралелизованог LB солвера – моделирање кретања деформабилног тела у флуиду .............................................................................................. 62 Слика 6.4 – Декомпозиција домена симулације дуж x осе ................................................. 64 Слика 6.5 – Компоненте функције расподеле које се размењују између два суседна GPU уређаја – случај симулације дводимензионалног струјања флуида .................................. 65 Докторска дисертација Списак слика Тијана Ђукић V Слика 6.6 – Компоненте функције расподеле које се размењују између два суседна GPU уређаја – случај симулације тродимензионалног струјања флуида .................................. 65 Слика 6.7 – Размена података између GPU уређаја када се примењује OpenMP приступ .................................................................................................................................................. 66 Слика 6.8 – Размена података између GPU уређаја када се примењује MPI приступ ..... 66 Слика 6.9 – Упоредни приказ времена извршавања за различит број чворова LB мреже; случај дводимензионалног струјања .................................................................................... 68 Слика 6.10 - Упоредни приказ времена извршавања за различит број чворова LB мреже; случај тродимензионалног струјања ..................................................................................... 68 Слика 6.11 - Дијаграм промене убрзања паралелизованог LB солвера у зависности од броја чворова LB мреже ......................................................................................................... 69 Слика 6.12 - Упоредни приказ убрзања програма за симулирање кретања деформабилног тела у флуиду за различит број чворова мреже солида .......................... 70 Слика 6.13 - Дијаграм промене убрзања паралелизованог LB солвера за симулирање вишекомпонентног струјања у зависности од броја чворова LB мреже ........................... 70 Слика 6.14 - Упоредни приказ времена извршавања за различит број чворова LB мреже и различит број GPU уређаја на којима се апликација извршава ...................................... 71 Слика 6.15 - Дијаграм промене убрзања паралелизованог LB солвера на више GPU уређаја у зависности од броја чворова LB мреже ............................................................... 72 Слика 6.16 - Strong scaling убрзање LB солвера извршаваног на више GPU уређаја, у поређењу са извршавањем на једном GPU уређају ............................................................. 73 Слика 6.17 - Weak scaling убрзање LB солвера извршаваног на више GPU уређаја, у поређењу са извршавањем на једном GPU уређају ............................................................. 74 Слика 7.1 – Положај аорте и срчаних артерија у људском организму ............................. 80 Слика 7.2 – Софтвер за генерисање геометрије на основу DICOM снимака ................... 80 Слика 7.3 – Расподела притиска у људској аорти ............................................................... 81 Слика 7.4 – Струјне линије и расподела брзина у људској аорти ..................................... 81 Слика 7.5 – Положај каротидне артерије у људском организму ....................................... 82 Слика 7.6 – Геометријски подаци за модел каротидне артерије (равански пресек) ........ 83 Слика 7.7 – Расподела притиска у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) .............................................................................. 84 Слика 7.8 – Струјне линије и расподела брзина у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) ...................................................... 84 Докторска дисертација Списак слика Тијана Ђукић VI Слика 7.9 – Расподела смичућих напона у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) ........................................................................ 85 Слика 7.10 – Положај коронарних артерија у људском организму ................................... 85 Слика 7.11 - Расподела притиска у коронарној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) .............................................................................. 86 Слика 7.12 - Струјне линије и расподела брзина у коронарној артерији – симулација применом LB методе .............................................................................................................. 87 Слика 7.13 - Расподела смичућих напона у коронарној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) ........................................................................ 87 Слика 7.14 – Моделирање процеса мешања две течности, у више временских тренутака .................................................................................................................................................. 89 Слика 7.15 – Пример 1 – Шема геометрије домена флуида ............................................... 92 Слика 7.16 – Пресек сферне честице у току деформисања и полуосе добијене елипсе .. 93 Слика 7.17 – Пример 1 – Поље брзине флуида и тренутни положај сферне честице ..... 94 Слика 7.18 – Пример 1 – Струјне линије.............................................................................. 94 Слика 7.19 – Пример 1 – Поређење коначног облика сферне честице у LB симулацији (црвена линија) са резултатима из литературе [42] (плава линија) ................................... 95 Слика 7.20 – Пример 1 – Промена Taylor-овог параметра облика током времена, за различите вредности параметра G (пуне линије представљају решења из литературе [42], а тачке решења LB симулације) ................................................................................... 95 Слика 7.21 - Пример 1 – Промена угла закошења честице током времена, за различите вредности параметра G (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) ......................................................................................................... 96 Слика 7.22 – Пример 1 - Промена Taylor-овог параметра облика током времена, за различите вредности параметра G, 2.0 (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) ............................................................... 97 Слика 7.23 - Пример 1 - Промена Taylor-овог параметра облика током времена, за различите вредности параметра G, 5 (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) ............................................................... 97 Слика 7.24 – Пример 1 – Деформација честице за G=0.1, за различите вредности коефицијента  : црвена линија - 5 ; зелена линија - 1 ; плава линија - 2.0 98 Слика 7.25 - Пример 1 - Промена Taylor-овог параметра облика током времена, за случај растерећења честице ................................................................................................... 98 Докторска дисертација Списак слика Тијана Ђукић VII Слика 7.26 - Пример 1 – Деформација честице током растерећења за G=0.075; плава линија - 0dt ; зелена линија - 1dt ; црвена линија - 5dt ........................................... 99 Слика 7.27 – Дискретизован модел и попречни пресек црвеног крвног зрнца .............. 100 Слика 7.28 - Пример 2 - Промена Taylor-овог параметра облика током времена, 1 ................................................................................................................................................ 100 Слика 7.29 - Пример 2 - Промена Taylor-овог параметра облика током времена, 5 ................................................................................................................................................ 101 Слика 7.30 - Пример 2 – Промена угла закошења црвеног крвног зрнца током времена, за различите вредности почетног угла (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) ................................................................................. 101 Слика 7.31 - Пример 2 – Промена облика црвеног крвног зрнца током времена, када је 1.0Ca ; плава линија – резултати из литературе [90]; црвена линија – LB симулација ................................................................................................................................................ 102 Слика 7.32 - Пример 2 – Промена облика црвеног крвног зрнца током времена, када је 5.0Ca ; плава линија – резултати из литературе [90]; црвена линија – LB симулација ................................................................................................................................................ 102 Слика 7.33 - Пример 2 – Струјне линије за случај када је 1.0Ca ................................. 103 Слика 7.34 – Пример 2 – Струјне линије за случај када је 5.0Ca ................................ 103 Слика 7.35 - Пример 3 – Шема геометрије домена флуида.............................................. 104 Слика 7.36 – Пример 3 - Поље брзине флуида и тренутни положај сферне честице .... 105 Слика 7.37 – Пример 3 – Промена облика сферне честице током времена, која је у центру правог канала у почетном тренутку; плава линија – резултати из литературе [42]; црвена линија – LB симулација .................................................................................. 105 Слика 7.38 - Пример 3 – Промена облика сферне честице током времена, која је у почетном тренутку померена за 063.0ry ; плава линија – резултати из литературе [42]; црвена линија – LB симулација .................................................................................. 106 Слика 7.39 - Пример 3 – Промена облика сферне честице током времена, која је у почетном тренутку померена за 13.0ry ; плава линија – резултати из литературе [42]; црвена линија – LB симулација ........................................................................................... 106 Слика 7.40 - Пример 4 - Поље брзине флуида и тренутни положај крвног зрнца – Први случај ..................................................................................................................................... 107 Слика 7.41 - Пример 4 - Поље брзине флуида и тренутни положај крвног зрнца – Други случај ..................................................................................................................................... 107 Докторска дисертација Списак слика Тијана Ђукић VIII Слика 7.42 - Пример 4 – Промена облика крвног зрнца током времена – Први случај; плава линија – резултати из литературе [91]; црвена линија – LB симулација.............. 108 Слика 7.43 - Пример 4 – Промена облика крвног зрнца током времена – Други случај; плава линија – резултати из литературе [91]; црвена линија – LB симулација.............. 108 Слика 7.44 – Шематски приказ нормалне и артерије са стенозом .................................. 109 Слика 7.45 – Пример 5 – Шема геометрије домена флуида ............................................. 110 Слика 7.46 – Пример 5 - Поље брзине флуида и почетни положај сферне честице ...... 110 Слика 7.47 - Пример 5 - Поље притиска флуида и почетни положај сферне честице ... 110 Слика 7.48 – Пример 5 – Поређење промене пада притиска током времена добијене применом LB солвера са резултатима из литературе [26] ................................................ 111 Слика 7.49 – Пример 5 - Симулација проласка сферне честице кроз артерију са стенозом ................................................................................................................................. 113 Слика 7.50 – Микроскопски снимак микрофлуидног чипа [94] ...................................... 115 Слика 7.51 - Пример 6 – Шема геометрије домена флуида.............................................. 117 Слика 7.52 - Пример 6 - Поље брзине флуида и почетни положај канцер ћелије ......... 117 Слика 7.53 - Пример 6 - Поље притиска флуида и почетни положај црвеног крвног зрнца....................................................................................................................................... 118 Слика 7.54 - Пример 6 - Симулација проласка канцер ћелије кроз микрофлуидни чип ................................................................................................................................................ 119 Слика 7.55 - Пример 6 - Симулација проласка црвеног крвног зрнца кроз микрофлуидни чип ............................................................................................................... 121 Слика 8.1 – Микроскопски снимак проласка црвених крвних зрнаца кроз стаклену цев [9] ........................................................................................................................................... 123 Слика 8.2 – Микроскопски снимак струјања крви кроз крвни суд пацова [9] ............... 123 Слика 8.3 – Упоредни приказ облика црвеног крвног зрнца добијених нумеричком симулацијом применом LB солвера (слике лево) и експериментом (слике десно); релативна почетна позиција крвног зрнца у односу на центар стаклене цеви одозго на доле - 00 y ; my 5.00  ; my 10  ................................................................................. 124 Слика 8.4 – Поље притиска флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку постављено на средину стаклене цеви .............................................. 124 Слика 8.5 - Поље брзине флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку померено за m5.0 у односу на средину стаклене цеви ................. 125 Слика 8.6 - Поље брзине флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку померено за m1 у односу на средину стаклене цеви .................... 125 Докторска дисертација Списак слика Тијана Ђукић IX Слика 8.7 – Кретање црвених крвних зрнаца кроз артерију са бифуркацијом .............. 126 Слика 8.8 - Шема геометрије артерије са бифуркацијом ................................................. 126 Слика 8.9 – Поље притиска флуида и почетни положај црвеног крвног зрнца ............. 127 Слика 8.10 - Симулација проласка црвеног крвног зрнца кроз артерију са бифуркацијом ................................................................................................................................................ 128 Слика 8.11 – Промена облика црвеног крвног зрнца током времена .............................. 129 Слика 8.12 – Фазе ембрионалног развоја зебра рибе [110] .............................................. 130 Слика 8.13 – Одрасла зебра риба и њен кардиоваскуларни систем ................................ 131 Слика 8.14 – Кардиоваскуларни систем у 3 дана старој зебра риби [111] ...................... 131 Слика 8.15 – Сплет репних вена зебра рибе [114] ............................................................. 132 Слика 8.16 – Поређење резултата добијених из експеримента и LB симулације за први почетни положај крвног зрнца; лево – микроскопски снимак зебра рибе, са означеним конкретним крвним зрнцем; у средини – издвојени облици конкретног крвног зрнца (црвено – експеримент, плаво – симулација); десно резултати добијени применом LB симулације ............................................................................................................................. 134 Слика 8.17 - Поређење резултата добијених из експеримента и LB симулације за други почетни положај крвног зрнца; лево – микроскопски снимак зебра рибе, са означеним конкретним крвним зрнцем; у средини – издвојени облици конкретног крвног зрнца (црвено – експеримент, плаво – симулација); десно резултати добијени применом LB симулације ............................................................................................................................. 135 Тијана Ђукић X Списак табела Табела 3.1 - Координате чворних тачака у природном координатном систему .............. 35 Табела 3.2 – Gauss-ова нумеричка интеграција за троугаоне коначне елементе ............. 41 Табела 7.1 – Подаци за LB солвер за пример 1 .................................................................... 93 Табела 7.2 – Подаци за LB солвер за пример 3 .................................................................. 105 Табела 7.3 – Подаци за LB солвер за пример 5 .................................................................. 109 Табела 7.4 - Вредности параметара за канцер ћелију дојке ............................................. 116 Табела 7.5 – Подаци за LB солвер за пример 6 .................................................................. 116 Табела Ц.1 – Карактеристике NVIDIA GeForce GT 630 .................................................... 142 Табела Ц.2 – Карактеристике NVIDIA Tesla C1060 ........................................................... 143 Табела Ц.3 – Карактеристике NVIDIA Tesla M2090 .......................................................... 143 Тијана Ђукић XI Списак коришћених ознака iA - тренутна површина посматраног елемента солида iA 0 - почетна површина посматраног елемента солида  naα - коефицијент експанзије n -тог степена α - векторски индекс B - леви Cauchy-Green-ов тензор деформације B0 - почетна конфигурација солида Bt - тренутна конфигурација солида hB - матрица извода интерполационих функција Ca - капиларни број c - микроскопска брзина sc - брзина звука у јединицама мреже xyD - Taylor-ов параметар облика солида  - фактор скалирања за густину t - фактор скалирања за време x - фактор скалирања за дужину uA - виртуални рад унутрашњих сила e - виртуална деформација  r - Dirac-ова делта функција u - виртуално померање ij - Kronecker-ов делта симбол e - вектор деформације ξ - вектор брзине iξ - вектор апсцисе Gauss-Hermite-ове квадратуре iξ - јединични вектор апсцисе Gauss-Hermite-ове квадратуре Докторска дисертација Списак коришћених ознака Тијана Ђукић XII  - Knudsen-ов број e - унутрашња енергија F t 0 - градијент деформације iF - дискретизован израз за спољашњу силу A F - сила реакције на промену површине солида B F - сила реакције на савијање мембране солида S F - сила реакције на површинско напрезање мембране солида V F - сила реакције на промену запремине солида S uF - унутрашње силе реакције на деформацију у равни мембране солида pF - фактор скалирања за конверзију бездимензионих величина у физичке величине f - функција вероватноће расподеле маса честица  0f - равнотежна функција расподеле  1f - ванравнотежна функција расподеле Nf - скраћена функција расподеле if - дискретизована функција расподеле  0 if - дискретизована равнотежна функција расподеле if - дискретизована функција расподеле у нумеричкој имплементацији  1 if - дискретизована ванравнотежна функција расподеле 21ii  - тренутни угао између два елемената солида 21 0 ii - почетни угао између два елемената солида G - бездимензиони коефицијент смицања G - Green-ова функција G - параметар јачине интеракције између два флуида  - коефицијент смицања g - вектор спољашњих сила  rH - Heaviside-ова функција ih - вектор који дефинише удаљеност тежишта елемента од координатне равни h - растојање између 2 тачке фиксне (Euler-ове) мреже Докторска дисертација Списак коришћених ознака Тијана Ђукић XIII kh - интерполационa функцијa  n αH - Hermite-ов полином n -тог степена  n iH - дискретизовани Hermite-ов полином n -тог степена  kji ,, - јединични вектори глобалног Декартовог координатног система  111 ,, kji - јединични вектори локалног Декартовог координатног система 1I - прва инваријанта 2I - друга инваријанта '1I - измењена прва инваријанта '2I - измењена друга инваријанта J - Јакобијан трансформације између глобалног и локалног система J - детерминанта градијента деформације j - количина кретања Bk - Boltzmann-ова константа AlK - параметар отпора локалној промени површине солида AtK - параметар отпора промени укупне површине солида BK - параметар отпора савијању мембране солида SK - параметар отпора површинском напрезању мембране солида VK - параметар отпора промени запремине солида sk - параметар отпора на смичућу деформацију у равни мембране солида k - параметар отпора на површинску деформацију у равни мембране солида iL - карактеристична дужина посматраног домена флуида у правцу i -те осе 12L - тренутно растојање између два чвора солида 12 0L - почетно растојање између два чвора солида l - средњи слободни пут 12l - јединични вектор који спаја два чвора солида  - однос вискозности унутрашњег и спољашњег флуида 1 - главно издужење у правцу осе x 2 - главно издужење у правцу осе y 3 - главно издужење у правцу осе z Докторска дисертација Списак коришћених ознака Тијана Ђукић XIV M - број елемената солида m - маса честице  - динамичка вискозност in - динамичка вискозност унутрашњег флуида out - динамичка вискозност спољашњег флуида N - број чворова коначног елемента lbN - резолуција, односно број чворова мреже флуида у правцу једне осе n - вектор нормале на зид in - јединични вектор нормале посматраног елемента солида  - кинематска вискозност  - оператор судара  - тежинска функција i - дискретни тежински коефицијент Gauss-Hermite-ове квадратуре P - тензор напона p - притисак флуида  tsr ,, - осе природног координатног система  - густина флуида S - тензор деформације σ - Cauchy-ев тензор напона T - апсолутна температура у Келвинима T - матрица трансформације  n T - вектор напона у тачки imT - компонента Piola-Kirchhoff-овог тензора напона прве врсте t - време τ - смичући напон  - време релаксације  - модификовано време релаксације k U - померање k -тог чвора коначног елемента у глобалном координатном систему k U - померање k -тог чвора коначног елемента у локалном координатном систему u - средња брзина флуида Bu - брзина граничне (Lagrange-ове) тачке Докторска дисертација Списак коришћених ознака Тијана Ђукић XV lbu - карактеристична брзина (у бездимензионим јединицама) V - тренутна запремина солида V0 - почетна запремина солида SW - густина деформационе енергије jw - јединични вектор који спаја тежиште са чвором посматраног елемента солида  - величина која је једнака mTkB k X - вектор положаја k -тог чвора коначног елемента  tBX - координате граничних (Lagrange-ових) тачака x - вектор положаја  CCC zyx ,, - координате тежишта солида dx - физичка величина у бездимензионом облику px - физичка величина у стандардном систему јединица Докторска дисертација Списак коришћених ознака Тијана Ђукић XVI Коришћене скраћенице BGK - Bhatnagar-Gross-Krook модел CT - Computed Tomography CUDA - Compute Unified Device Architecture GPU - Graphics Processing Unit GPGPU - general purpose GPU Kn - Knudsen-ов број LB - lattice Boltzmann Ma - Mach-ов број МКЕ - метод коначних елемената PRACE - Partnership for Advanced Computing in Europe RBC - red blood cell Re - Reynolds-ов број SRT - апроксимација једноструког времена релаксације (енг. single relaxation time approximation) Тијана Ђукић 1 1 Увод Кардиоваскуларни систем је један од најважнијих система у људском организму. Разумевање бројних појава које се дешавају у оквиру овог система је значајно да би се лакше анализирали процеси настанка и развоја многих болести, као што су настанак плака, формирање тромбозе, развој атеросклерозе итд. Поред наведених примена у анализи болести које су примарно везане за кардиоваскуларни систем, моделирање струјања крви је битно и због разумевања многих других процеса као што су на пример снабдевање организма кисеоником, циркулисање малигних ћелија кроз организам и спровођење интравенозних лекова, матичних ћелија и генетског материјала до циљаних органа. Због недовољног разумевања ових појава, ограничене су могућности дијагностике болести, планирања и предузимања одговарајућих превентивних мера, анализе транспорта лекова и дизајнирање ефикаснијих терапијских процедура. Због тога су истраживања струјања крви клинички веома значајна, а компјутерске методе се веома успешно могу применити управо у овој области. Математички и нумерички модели могу да допринесу добијању много нових и корисних информација, а при томе је драстично јефтиније и једноставније спровести нумеричку симулацију, него спроводити комплетан експеримент. Нумеричке симулације су до сада успешно примењене у многим областима, укључујући и моделирање и анализу процеса као што су испорука нанолекова или стварање атеросклерозе [1,2]. Када се разматра струјање крви кроз кардиоваскуларни систем, потребно је раздвојити математичке моделе који разматрају струјање флуида на макроскопском и микроскопском нивоу. Наиме, природа струјања се значајно разликује у зависности од пречника крвног суда. Ако је крвни суд пречника већег од 200 микрометара, струјање се може довољно тачно моделирати ако се претпостави да је крв хомогени флуид. Међутим, када се разматра струјање крви кроз капиларне крвне судове, чији је пречник од 4 до 100 микрометара, потребно је узети у обзир и структуру крви на микроскопском нивоу. Људска крв је сложен флуид који се састоји од ћелија (формираних ћелијских елемената) и крвне плазме. Постоји неколико типова крвних ћелија, а то су еритроцити, леукоцити и тромбоцити. Око 40-45% запремине крви чине еритроцити (црвена крвна зрнца), која имају велики утицај на карактеристике струјања крви на микроскопском нивоу. Механички гледано, црвено крвно зрнце се може посматрати као капсула, чија је унутрашњост испуњена нестишљивим вискозним флуидом (хемоглобином), а ограничена је јако флексибилном мембраном, коју чине два слоја липида и танак цитоскелет сачињен од међусобно повезаних протеина [3]. Густина унутрашњег флуида (хемоглобина) и околне крвне плазме је иста, али је зато вискозност унутрашњег флуида 5 до 6 пута већа од вискозности крвне плазме [4,5]. За разлику од већине биолошких ћелија, еритроцити имају веома редуковану унутрашњу структуру и немају Докторска дисертација Поглавље 1 Тијана Ђукић 2 ћелијско једро. Њихова улога у људском организму је да узимају кисеоник током проласка крви кроз плућа, да транспортују кисеоник до свих осталих делова тела и да омогућавају размену кисеоника и угљен диоксида у околном ткиву преко капиларних крвних судова. Због своје структуре, улоге и различитих услова струјања којима су изложена црвена крвна зрнца током кретања кроз кардиоваскуларни систем, она морају да буду веома деформабилна. Та способност деформације највише долази до изражаја приликом проласка кроз мале крвне судове и сужења која постоје у капиларном циркулаторном систему. Пречник ових ћелија је у микроциркулацији истог реда величине као и пречник крвног суда и због тога у нумеричким симулацијама оваквих појава није довољно третирати крв као хомогени флуид, већ се црвена крвна зрнца морају третирати као независне ћелије које су уроњене у крвну плазму. Анализа динамике кретања ових ћелија појединачно је управо због њихове велике улоге један од најважнијих проблема у физиологији и биомеханици. Многи аутори су анализирали понашање црвених крвних зрнаца, како експериментално, тако и теоријски. Експериментално је посматрано понашање вештачких капсула [6,7], али и крвних зрнаца [8,9]. Механичке особине еритроцита, са посебним фокусом на карактеристике ћелијске мембране, проучавали су многи аутори, како у прошлом веку [10,11], тако и у претходној деценији [12,13,14]. Када је у питању нумеричко моделирање понашања црвених крвних зрнаца у флуиду, у литератури је примењено неколико различитих метода. Pozrikidis [15] је применио методу граничних елемената (енг. Boundary Element Method – BEM) да моделира кретање црвених крвних зрнаца. Ову методу су користили и други аутори [16,17]. Ramanujan и Pozrikidis [18] су проширили овај приступ и разматрали неколико различитих материјалних модела којима се дефинише однос силе и деформације честица. Поред тога, у овом раду је поред ефекта спољашњег флуида, анализиран и ефекат флуида унутар честице, који има различиту вискозност у односу на спољашњи флуид. Слични утицај су разматрали и Diaz и други [19]. За моделирање утицаја флуида на солид и обрнуто, Eggleton и Popel [20] су применили методу уроњене границе (енг. Immersed boundary method). Sui и други [21] су унапредили овај приступ тако што су додали побољшања густине мреже, да би се прецизније моделирало кретање капсуле са већом густином мреже у делу где је честица. Такође су коришћене и дискретне честичне методе за моделирање ових појава. Boryczko и други [22] су моделирали солид применом класичне механике континуума, док су флуид моделирали као скуп честица. Са друге стране, по угледу на [23], Tsubota и други [24] су и флуид и солид представили као скуп честица и анализирали њихове међусобне интеракције и на тај начин моделирали кретање крвних зрнаца кроз крвну плазму. Симулације већег броја црвених крвних зрнаца и њихова међусобна интеракција су биле тема одређеног броја радова [25,26,27]. Међутим, у већини ових радова је разматрано кретање честица и крвних зрнаца у једноставнијим геометријским условима, где не долазе до изражаја особине крвних зрнаца да драстично мењају свој облик. Циљ овог рада је да моделира кретање крвних зрнаца и сферних честица кроз комплексне геометријске домене, као што су артерије са стенозом, бифуркацијом, као и кроз микрофлуидне чипове за сепарацију канцер ћелија, али и кроз реалне домене добијене експерименталним снимањем протока крви кроз посебну врсту рибе, тзв. зебра рибу. У овом раду ће бити примењене три нумеричке методе, које су међусобно повезане у једну целину. За моделирање струјања флуида коришћена је lattice Boltzmann (скраћено LB) метода [28,29]. Ово је релативно нова метода која се користи за решавање Навије- Докторска дисертација Поглавље 1 Тијана Ђукић 3 Стоксових једначина. Она је одабрана због ефикасности у симулацијама струјања вишекомпонентних флуида, као и струјања у микро и нано каналима [30,31], али и због релативно једноставне имплементације и паралелизације. Такође, ова метода је већ успешно примењена за моделирање кретања честица LDL-а (енг. low-density lipoprotein) кроз крвоток [32], а примењена је и за анализу и предвиђање кретања нанолекова кроз крвоток [33]. За моделирање понашања деформабилних тела узети су у обзир утицаји различитог типа деформација. За реакције честице на деформацију која изазива промену запремине, површине и савијање мембране коришћене су дискретне формуле за одређивање силе реакције, по угледу на [25]. За моделирање површинског напрезања мембране коришћена је метода коначних елемената [2,34,35], и посебни материјални модел мембране црвеног крвног зрнца, који је коришћен у литератури [36] и који је развио Skalak [37]. За моделирање интеракције између флуида и солида коришћена је Immersed Boundary метода, коју је први представио Peskin [38]. Она је одабрана зато што се теоријске основе ове методе добро уклапају са начином на који је математички развијена LB метода и заједно чине један складан спрегнут систем који се успешно примењује за моделирање како чврстих [39,40,41], тако и деформабилних честица које су уроњене у флуид [36,42]. Главна предност ове методе је што приликом њеног коришћења није потребно вршити додатне модификације Навије-Стоксових једначина, већ се ефекат уроњене честице узима у обзир преко дејства спољашње силе. Поред тога, ова метода се успешно примењује када треба моделирати танке и деформабилне мембране које раздвајају два флуида, као што је то случај који се разматра у овом раду. Још једна предност је што су код ове методе домени флуида и ћелијске мембране раздвојени што се тиче нумеричког решавања, па се може користити различита дискретизација за ова два домена. Нумеричке симулације као што је моделирање кретања деформабилних честица у флуиду, јако су захтевне када су у питању компјутерски ресурси. Због тога је приликом имплементације нумеричког модела у овом раду посебна пажња посвећена техникама паралелизације. Компјутерски програм је развијан тако да поред рада на стандардним PC рачунарима, може да се извршава и на компјутерским кластерима са великим бројем процесора, као што су Tesla суперкомпјутери, које производи светска компанија NVIDIA. Овим суперкомпјутерима је приступано у оквиру PRACE (енг. Partnership for Advanced Computing in Europe) међународних пројеката, на којима је био или је тренутно ангажован аутор овог рада. Коришћени су принципи GPU програмирања и CUDA архитектура коју је развила NVIDIA [43]. Постигнута су значајна убрзања развијеног софтвера, тако да се нумеричке симулације на тај начин могу извршавати за само неколико минута, наспрам пар сати колико је потребно за извршавање исте симулације на стандардном PC рачунару, што омогућава практично интерактивно праћење резултата у реалном времену. Поређењем резултата спроведених нумеричких симулација са резултатима који су објављени у литератури, као и са експерименталним подацима, биће показано да се развијена нумеричка метода и софтвер могу успешно примењивати у симулацијама реалних проблема у биомедицинском инжењерингу. Овај рад је организован на следећи начин. У поглављу 2 описана је lattice Boltzmann нумеричка метода која је коришћена за моделирање струјања флуида. У поглављу 3 описан је поступак моделирања реакције солида на деформације. У овом поглављу посебна пажња посвећена је методи коначних елемената, која је коришћена за Докторска дисертација Поглавље 1 Тијана Ђукић 4 моделирање реакције деформабилног тела на површинско напрезање, при чему је примењен посебан материјални модел који се користи за хипереластичне материјале, какав је и мембрана црвеног крвног зрнца. Спрезање домена солида и флуида, односно моделирање солид-флуид интеракције, остварено је помоћу Immersed boundary методе, која је описана у поглављу 4. У поглављу 5 описани су приступи који су коришћени за моделирање струјања вишекомпонентних флуида. Паралелизација развијеног софтвера описана је у поглављу 6. У поглављу 7 приказани су примери у којима су поређени резултати добијени применом развијеног софтвера са резултатима из литературе. Поређење резултата добијених нумеричким симулацијама са експерименталним резултатима приказано је у поглављу 8. У поглављу 9 дат је кратак закључак целог рада. Тијана Ђукић 5 2 Моделирање струјања флуида Постоји велики број нумеричких метода за моделирање струјања флуида. Неке од стандардних метода које се користе у рачунској динамици флуида (енг. Computational Fluid Dynamics – CFD) су метода коначних елемената, коначних разлика, коначних запремина. Суштина свих метода је решавање система једначина (Навије-Стоксове једначине и једначине континуитета), да би се као резултат добило поље макроскопских физичких величина – брзине, притиска, смичућег напона. Поред стандардних метода које се базирају на једначинама које су изведене за континуум, постоје и дискретне методе, као што су lattice Boltzmann (скраћено LB), DPD (енг. Dissipative Particle Dynamics), SPH (енг. Smoothed-Particle Hydrodynamics). У овом раду је коришћена LB метода. Она је изабрана због ефикасности у симулацијама струјања у микро и нано каналима [30,31] и у симулацијама кретања солида у домену флуида [25,40,44]. Поред тога, велика предност ове методе је што нема потребе за решавањем диференцијалних једначина или система једначина, па је имплементација релативно једноставна и могуће је извршити паралелизацију. Lattice Boltzmann метода припада класи проблема такозваних Cellular Automata (CA). Код свих оваквих проблема физички систем се посматра на идеализован начин, тако да су време и простор дискретизовани, односно цео простор је сачињен од великог броја идентичних ћелија [45]. Флуид се посматра као скуп фиктивних честица и проучава се динамика кретања ових честица, кроз њихове међусобне сударе и даљу пропагацију у посматраном домену. Дефинише се посебна функција пропагације, која зависи од стања суседних честица и идентичног је облика за све честице, односно све чворове дискретизоване мреже (енг. lattice mesh). Стање у свим чворовима се ажурира синхронизовано кроз низ итерација, у дискретним временским корацима. Једначине LB методе се изводе из континуалне Boltzmann-ове једначине. Поступак извођења и трансформације које је потребно спровести да би се добио облик погодан за нумеричко решавање објашњене су у наставку овог поглавља. Такође, може се доказати да се моделирањем флуида на микроскопском нивоу и праћењем динамике дискретних честица може моделирати и струјање флуида на макроскопском нивоу. Поступак којим се показује да се из једначина LB методе могу извести Навије-Стоксове једначине такође је описан у овом поглављу. Докторска дисертација Поглавље 2 Тијана Ђукић 6 2.1 Boltzmann-ова једначина У оквиру LB методе флуид се посматра као скуп честица и очигледно је да је потребно посматрати велики број честица, које међусобно интерагују и сударају се у посматраном простору. Уколико би се посматрала свака честица појединачно, то би значило да се флуид у потпуности анализира на микроскопском нивоу. Овакав приступ није погодан јер је овакво праћење тешко изводљиво, а са друге стране, за проблеме динамике флуида, који су у овом случају предмет интересовања, најбитније је пратити промене у систему на макроскопском нивоу. Зато је у LB методи основна величина функција расподеле f , која описује понашање и кретање честица у простору. Функција расподеле је дефинисана тако да  tf ,,ξx представља вероватноћу да се честица налази унутар елемента dxd око тачке  ξx, , у временском тренутку t , при чему x и ξ означавају вектор положаја и вектор брзине честице, респективно. Основна једначина од које почињу сва даља извођења је Boltzmann-ова једначина, парцијална диференцијална једначина која важи за континуум и која заправо представља једначину баланса функције расподеле. Ако су узете у обзир спољашње силе, Boltzmann-ова једначина гласи:          ξ g x ξ f m f t f (2.1) где је са g означен вектор спољашњих сила, а са  је означен оператор судара (енг. collision operator), који представља промену функције расподеле услед судара честица. Оператор судара  се представља функцијом јако комплексног облика, па је уведен поједностављени модел, који су предложили Bhatnagar, Gross и Krook [46]. Уведена је претпоставка да је ефекат судара честица такав да доводи флуид „ближе“ равнотежном стању. Овај модел је познат као апроксимација једноструког времена релаксације (енг. single relaxation time approximation) или Bhatnagar-Gross-Krook (BGK) модел. Оператор  је дефинисан на следећи начин:   01 ff   (2.2) где је  време релаксације, односно просечно време које протекне између 2 судара, а  0f равнотежна функција расподеле, тзв. Maxwell-Boltzmann-ова функција расподеле. Као што је већ поменуто, LB метода се заснива на принципима који важе за гасове, али су уведена нека ограничења да би се ови принципи применили у динамици флуида. Maxwell-Boltzmann-ова функција расподеле вероватноће је преузета из статистичке Докторска дисертација Поглавље 2 Тијана Ђукић 7 механике, где се користи да би се одредиле брзине молекула. Детаљи везани за дефинисање и извођење Maxwell-Boltzmann-ове функције расподеле су овде изостављени, а могу се наћи у литератури [47,48]. Коначан израз равнотежне функције расподеле гласи:                        t t t t tf D ,2 , exp ,2 , ,, 2 2/ 0 x ξxu x x ξx   (2.3) У овом изразу: D је број физичких димензија ( 3D за тродимензионални простор) mTkB K J kB 231038,1  је Boltzmann-ова константа T је апсолутна температура у Келвинима (K) m је маса честице (у даљем тексту, сматра се да је 1m ) Коначно, BGK модел Boltzmann-ове једначине, који важи за континуум, дефинисан је следећом једначином:   01 ffff t f          ξ g x ξ (2.4) 2.2 Извођење Навије-Стоксове једначине Потребно је применити одређени поступак да би се Boltzmann-ова једначина која је дефинисана у претходном поглављу трансформисала у систем једначина који се користи у класичној динамици флуида – једначину континуитета и Навије-Стоксову једначину. Постоје два начина да се добију тражене једначине. Један је да се једначина (2.4) пројектује на Hermite-ову базу (енг. Hermite polynomials basis projection) и да се онда даље трансформишу добијени изрази. Ову процедуру је први предложио и спровео Grad [49,50], па се и назива Grad-ових 13 моментних једначина (енг. Grad’s 13- moments equations). Други приступ су предложили Chapman и Cowling [51] и назива се Chapman-Enskog експанзија на више скала (енг. Chapman-Enskog multiscale expansion). У оквиру тог поступка се једначине раздвајају на више скала тако што се сви изводи, функција расподеле и све величине развијају у функцији различитих степена Кнудсеновог броја. Тако добијени изрази се даље анализирају за различите степене Кнудсеновог броја и добијају се тражене једначине. У овом раду је примењен Докторска дисертација Поглавље 2 Тијана Ђукић 8 једноставнији поступак, који при томе задржава основну идеју. Овај поступак је применио Huang [52] и објашњен је у литератури [41,53,54]. Најпре је потребно пројектовати цео систем на Hermite-ову полиномску базу. Пројекција функције расподеле дата је следећом једначином:          ta n tf n n n , ! 1 ,, 0 xξξξx αα    H (2.5) У изразу (2.5):  ξ је тежинска функција  n αH и  naα представљају Hermite-ов полином n -тог степена и коефицијент експанзије n -тог степена, респективно α је векторски индекс;  n ,, 21α , при томе важи Einstein-ова конвенција о сабирању по поновљеним индексима. Rodrigues-ова формула дефинише Hermite-ов полином n -тог степена, у D- димензионом простору:        ξ ξ      nn n 1H (2.6) Тежинска функција  ξ се дефинише:             2 exp 2 1 2 2/    D ξ (2.7) Коефицијент експанзије се дефинише:        ξξxξxα nn tfdta H ,,, (2.8) Као што се може видети из једначине (2.5), пројекцијом на Hermite-ову базу функција расподеле се развија у бесконачни ред. Међутим, у макроскопским једначинама (законима одржања масе и количине кретања) фигурише само 3 величине - густина флуида  , количина кретања u и енергија e . Стога се може извести закључак да се могу занемарити чланови после одређеног степена. Тиме се функција расподеле „скраћује“, односно губе се одређене информације које је оригинална функција Докторска дисертација Поглавље 2 Тијана Ђукић 9 расподеле садржала. У литератури је показано да те информације заправо и нису неопходне за симулације струјања флуида [53]. Ако се занемаре сви чланови изнад четвртог реда (ако је 4N ), из BGK једначине се могу извести Навије-Стоксове једначине за стишљив термални флуид. Слично, ако се изврши скраћивање до трећег реда, добијају се идентичне једначине одржања масе и количине кретања, као за 4N . Али, губитак информација из функције расподеле доводи до тога да једначина одржања енергије није више идентична и више не може довољно тачно да симулира топлотна струјања. Оваква формулација ипак сасвим реално може да репрезентује изотермална струјања стишљивог флуида. Скраћивањем функције расподеле и занемаривањем чланова изнад другог степена (за 2N ), добија се најчешће коришћен lattice Boltzmann модел. Овај модел важи само за изотермалне и слабо стишљиве флуиде. Најчешће се користи због тога што репрезентује и најчешће коришћен облик Навије-Стоксових једначина у рачунској динамици флуида. У наставку овог рада извођење ће бити спроведено управо за овај модел – за изотермално струјање нестишљивог флуида. Због свега наведеног, у поступку који је предложио Grad, довољно је израчунати дефинисане коефицијенте до другог степена. Али најпре је потребно дефинисати макроскопске величине преко функције расподеле. Пошто LB метода описује флуид на молекуларном нивоу, карактеристике флуида са аспекта континуума су имплицитно садржане у моделу. Функција расподеле се због тога користи за дефинисање макроскопских величина. Густина се дефинише интегралом:     tfdt ,,, ξxξx (2.9) Количина кретања се дефинише интегралом:       tfdtt ,,,, ξxξξxuxj  (2.10) Интеграција се врши на целом домену брзине. У изразу (2.10) u представља поље средње брзине. Остале макроскопске величине се рачунају као функције вишег степена, па се ради прегледности не пишу у функцији средње брзине, већ преко одступања од средње брзине. Због тога се уводи појам микроскопске брзине, као поље одступања од средње брзине    tt ,,, xuξξxc  (2.11) Докторска дисертација Поглавље 2 Тијана Ђукић 10 Густина унутрашње енергије (енг. internal energy density) се сада може изразити:   tfde ,,2 1 2 ξxcc (2.12) где је 22 cc  . У складу са једначином која важи за једноатомски гас, и у случају флуида који се разматра на овај начин, унутрашња енергија се следећом једначином може повезати са температуром:  2 D e  (2.13) Односно:  2 D e  (2.14) Поред ових величина, потребно је дефинисати још једну величину која фигурише у законима одржања. Ова величина је тензор напона, који се дефинише на следећи начин:   tfd ,,ξxcccP (2.15) где је cc тензорски производ c и c . Хидростатички притисак се може дефинисати преко везе са унутрашњом енергијом: D e p 2  (2.16) Коришћењем израза (2.14) долази се до једначине стања идеалног гаса и то је уједно и начин на који се у LB методи уводи притисак у систем једначина: m Tk p B  (2.17) Докторска дисертација Поглавље 2 Тијана Ђукић 11 Пошто су дефинисане макроскопске величине преко функције расподеле, сада се на основу формуле (2.6) могу израчунати Hermite-ови полиноми до другог степена:   10 H (2.18)     1H (2.19)      2H (2.20) Ако се у општој једначини за коефицијент експанзије (2.8) замене макроскопске величине (дефинисане једначинама (2.9), (2.10) и (2.15)) добијају се следеће једначине за коефицијенте експанзије до другог степена:   0a (2.21)    ua  1 (2.22)       uuPa 2 (2.23) У складу са овим закључцима, функција расподеле се може записати на следећи начин:            ta n tftf n N n nN , ! 1 ,,,, 0 xξξξxξx αα   H (2.24) при чему се сви чланови реда вишег од N занемарују. Сада се може дефинисати и равнотежна функција расподеле, развијена у ред:      n N n n a n f αα 0 0 0 ! 1    H (2.25) где су коефицијенти  na α0 дефинисани једначином:      nn fda αα H 0 0  (2.26) Докторска дисертација Поглавље 2 Тијана Ђукић 12 У овом изразу (и у свим наредним изразима), сматра се да је 2N , као што је већ наглашено. Коефицијенти  na0 се могу лако одредити применом Гаусове интеграције и заменом израза за макроскопске величине:   00a (2.27)    ua  1 0 (2.28)       1 2 0  uua (2.29) Како се у овом раду разматра струјање нестишљивог, изотермалног флуида, температура се не узима у обзир, па закон одржања енергије нема никакав утицај и сматра се да је 1 . Због тога коефицијент  20a постаје:     uua  2 0 (2.30) Заменом вредности коефицијената  na0 , као и Hermite-ових полинома дефинисаних једначинама (2.18)-(2.20), добија се коначни израз за равнотежну функцију расподеле:             uuf i 20 2 1 1 Huξ (2.31) Сада се Boltzmann-ова једначина (2.4) може пројектовати на Hermite-ову базу. Најпре ће бити пројектована само лева страна једначине, а за десну ће потом бити примењено Chapman-Enskog растављање да би се добио жељени систем једначина.                                                 ξξgξx ξ ξξξxξ x ξξξx ξξξx ξ g x ξ dtfdtfdtf t dtf t nnn n HHH H ,,,,,, ,, (2.32) Трећи члан израза (2.32) може се трансформисати применом правила о изводу производа:                          ξξ ξ gξxξξgξx ξ ξξgξx ξ dtfdtfdtf nnn HHH ,,,,,, (2.33) Докторска дисертација Поглавље 2 Тијана Ђукић 13 Први члан израза (2.33) се даље трансформише применом Green-ове формуле и постаје једнак нули:               0,,,, dtfdtf nn nξgξxξξgξx ξ HH (2.34) На основу дефиниције Hermite-ових полинома (једначина (2.6) и (2.7)) могу се извести следеће релације [53,55]:      nnn HHH ξ ξ    1 (2.35)    1-nHH    n ξ (2.36) где се подразумева сумирање по индексу  . Ако се изрази (2.33), (2.35) и (2.36) примене на израз (2.32), добија се:                           ξξgξxξξxξ x ξξξx dtfdtfdtf t nn 1-n1-n HHHH  ,,,,,, 1 (2.37) Ако се сада примени и израз за коефицијент експанзије (2.8), добија се коначни израз за леву страну Boltzmann-ове једначине пројектован на Hermite-ову базу:        111           nnnn aaaa t  g xx (2.38) Сада се примењује Chapman-Enskog растављање на величине различитог реда да би се добиле познате једначине динамике флуида. Ако се најпре претпостави да је функција расподеле једнака равнотежној вредности функције расподеле   0ff  , десна страна једначине (2.4) једнака је нули, па BGK једначина пројектована на Hermite-ову базу гласи:         010 1 0 1 00            n nnn a aa t a   g xx (2.39) Докторска дисертација Поглавље 2 Тијана Ђукић 14 Заменом вредности коефицијената експанзије за 0n и 1n добијају се две основне једначине рачунске динамике флуида – закони одржања масе и количине кретања:     0      u x  t (2.40)     guuI x u        t (2.41) Добијене једначине (2.40) и (2.41) су Ојлерове једначине за флуид, а циљ извођења је да се добију Навије-Стоксове једначине. Зато се дефиниција функције расподеле f мења и сматра се да није једнака  0f , већ се може разложити на равнотежни (енг. “equilibrium” part) и ванравнотежни део (енг. “off-equilibrium” part), при чему је ванравнотежни део много мањи и сматра се малим одступањем од равнотежног стања. Ова дефиниција се може представити једначином:    10 fff  (2.42) при чему је са  означен Кнудсенов број ( Kn ) који се дефинише као однос средњег слободног пута честица и карактеристичне дужине. L l  (2.43) У изразу (2.43) са L је означена карактеристична дужина посматраног домена, а scl  је средњи слободни пут честица (просечни пут који честица пређе између 2 судара). Вредност sc се дефинише као карактеристична брзина честица, која се још назива и брзина звука у јединицама мреже (енг. speed of sound in lattice units). Вредност овог параметра се дефинише у оквиру поступка дискретизације, што ће бити тема наредног поглавља. Јасно је из дефиниције ванравнотежног дела функције расподеле да мора да важи    01 ff  . Ако би Кнудсенов број био близу 1 или већи од 1, слободан пут честица би био истог реда величине као карактеристична дужина посматраног домена, а самим тим претпоставка да се симулира флуид као континуум не би била валидна. Пошто је циљ да се симулира флуид као континуум, изводи се закључак да целокупни приступ LB методе и BGK модел заправо важе само у области малог Knudsen-овог броја, па је неопходно да важи 1 . Докторска дисертација Поглавље 2 Тијана Ђукић 15 Ако се узме у обзир израз за функцију расподеле (2.42), десна страна једначине (2.4) се такође може пројектовати на Hermite-ову базу:        nn adtf 1 1 1,, 1   ξξξx H (2.44) Целокупна BGK једначина пројектована на Hermite-ову базу сада гласи:          nn nnn aa aa t a 1 1 0 1 0 1 00 1               g xx (2.45) при чему је коефицијент експанзије за неравнотежну функцију расподеле дефинисан:      nn fda αα H 1 1  (2.46) а ванравнотежна функција расподеле развијена у ред се дефинише:      n N n n a n f αα 1 0 1 ! 1    H (2.47) Посредним диференцирањем парцијални извод коефицијента  na0 по времену може се трансформисати (при томе се узимају у обзир претходно добијене Ојлерове једначине (2.40) и (2.41)):                                               uuI x g jx u u j      nn nnn aa t a t a t a 00 000 (2.48) На основу већ наведених израза за коефицијенте  na0 (једначине (2.27)-(2.30)), лако је одредити парцијалне изводе по  и j :   1 0 0     a ;   0 1 0     a ;          20a (2.49)   0 0 0    j a ;         j a 10 ;       uu j a    20 (2.50) Докторска дисертација Поглавље 2 Тијана Ђукић 16 Сада је једноставно израчунати вредности коефицијената  na1 :   001 a ;   011 a ;   021 a ;    Sa 2 2 1  (2.51) где је са S означен тензор деформације:   ,, 2 1 uuS  (2.52) У овом изразу са  ,u је означен парцијални извод   x u   . Треба напоменути једну важну апроксимацију која је примењена приликом израчунавања коефицијента  2 1a . Уколико би се израчунавао поступно, израз би гласио:      321 Ma22 1 OSuuu x Sa          (2.53) тако да се занемарују чланови вишег реда. У овом изразу, Ma је Mach-ов број, који се дефинише као однос карактеристичне брзине флуида и „брзине звука“ sc . sc u Ma (2.54) Целокупна LB метода и BGK модел важе само у области малог Mach-овог броја. Коришћењем ове апроксимације, чланови вишег реда у једначини (2.53) се могу занемарити и на тај начин се добија већ написани израз (2.51). Коришћењем ове апроксимације и свих претходно одређених величина, може се израчунати коначан израз за ванравнотежну функцију расподеле  1f :      Sf 21 H (2.55) У Додатку А је описана процедура којом се Boltzmann-ова једначина трансформише у једначине одржања масе и количине кретања. Та процедура је врло слична овој која је већ спроведена приликом пројектовања на Hermite-ову базу, па је због тога Докторска дисертација Поглавље 2 Тијана Ђукић 17 изостављена из главног дела овог рада. Као крајњи резултат ове процедуре, добијају се следеће једначине:   0      u x   t (2.56)     guuP x u        t (2.57) Да би се добио потпун систем једначина, потребно је одредити тензор напона P . Применом истог принципа Chapman-Enskog растављања, који је већ коришћен приликом израчунавања функције расподеле, тензор напона се може раставити на 2 компоненте:    10 PPP  (2.58) На основу једначине (2.15), добија се:                ξuuuξξuξξ ξuξuξP dff dff 10 10     (2.59) Даље се могу трансформисати обе компоненте појединачно:                ξdfuuuuP nn   0HHHH   112 (2.60) Заменом 0n и 1n , имајући у виду претходно израчунате вредности коефицијената експанзије (2.27)-(2.30) и (2.51), добија се:     0P (2.61)    SP 2 1  (2.62) Укупна вредност тензора напона износи: SP  2 I (2.63) где је са   означена динамичка вискозност. Докторска дисертација Поглавље 2 Тијана Ђукић 18 Ако се сада добијени израз за тензор напона замени у једначину (2.57), добија се познати систем једначина – једначина континуитета и Навије-Стоксова једначина за нестишљив флуид:   0      u x   t (2.64)     gSIuu x u        2p t (2.65) 2.3 Дискретизација Boltzmann-ове једначине BGK Boltzmann-ова једначина која је дефинисана у поглављу 2.1 се односи на континуум и подразумева да постоји континуално поље брзина. Овакав облик једначине није погодан за нумеричко решавање, па је потребно извршити дискретизацију. При томе се мора водити рачуна да систем једначина остане „нетакнут“ и да и даље буду задовољене Навије-Стоксове једначине, да би цела метода могла да се користи за симулације струјања флуида. Поступак дискретизације се спроводи у два корака. Први корак подразумева дискретизацију поља брзина, а затим се примењује просторна и временска дискретизација. 2.3.1 Дискретизација поља брзина Једначина којом су дефинисани коефицијенти експанзије (једначина (2.8)) може се трансформисати на следећи начин:            trdtfdta nNn ,,,,, ξxξξξξxξxα   H (2.66) где је са  tr ,,ξx означен полином степена не већег од N2 , који зависи од брзине ξ . Интеграл (2.66) се може израчунати тако што се трансформише у тежинску суму, коришћењем поступка Gauss-Hermite квадратуре [56,57]. Докторска дисертација Поглавље 2 Тијана Ђукић 19 У овом конкретном случају, за полином  tr ,,ξx , Gauss-ова квадратура тежи да што приближније одреди вредност интеграла и то тако што се дефинише оптималан скуп апсциса qii ,,2,1, ξ , при чему важи:         q i ii rdr 1 ξξξξ  (2.67) где је qii ,,2,1,  скуп константних тежинских коефицијената. За ово правило се каже да има алгебарски степен тачности m , ако за било који полином r , максимално m -тог степена, важи егзактна једнакост горње једначине. Алгебарски степен тачности m је значајан зато што се помоћу њега дефинише број апсциса које ће се користити за тачно одређивање Hermite-ових коефицијената. У зависности од типа струјања које треба да се симулира, потребан је различит број коефицијената. За топлотно струјање стишљивог флуида узима се 9m , за изотермално струјање стишљивог флуида узима се 7m , а за случај који се разматра у овом раду – изотермално струјање нестишљивог флуида, алгебарски степен тачности је 5m . Наравно, што је веће m , већи су и меморијски захтеви и потребно је више времена за извршавање симулације. У Додатку Б дати су Gauss-Hermite-ови тежински коефицијенти и апсцисе за алгебарски степен тачности 5, који се користи за решавање дводимензионалних и тродимензионалних проблема струјања нестишљивог флуида. Једначина (2.66) се своди на:            i n i N q i i i q i ii n tftrta ξξx ξ ξxx αα H,,,,, 1 0 1 0          (2.68) На овај начин, дискретизована функција расподеле   1,,0,,,  qitf i N ξx потпуно одређује  tf N ,,ξx и све особине величина које се израчунавају на основу Nf остају непромењене. Као што се може видети из Додатка Б, вектори који дефинишу апсцисне осе немају јединичне координате, а пожељно је да то буде случај. Пошто су поједине вредности ирационални бројеви 3ij , неопходно је увођење корекционог фактора да би се добиле тражене вредности. Управо због тога се дефинише константа sc (тзв. брзина звука мреже), која је већ поменута у једначини (2.54). Ова константа зависи од избора дискретних брзина. За моделе који су разматрани у овом раду и који су приказани у Додатку Б, узима се да је: 3 12 sc (2.69) Докторска дисертација Поглавље 2 Тијана Ђукић 20 Да би се добиле јединичне вредности апсцисних координата, уводи се смена: isi c ξξ  (2.70) Већ је поменуто да се sc често у литератури назива брзина звука мреже (енг. lattice speed of sound), али заправо представља само фактор скалирања. Да би се исправно дискретизовала и функција расподеле, потребно је увести дискретне Hermite-ове полиноме:   10 iH (2.71)    ii  1H (2.72)     22 siii cH (2.73) Увођење фактора скалирања доводи и до промена у изразима за коефицијенте експанзије, који сада гласе:   0a (2.74)    ua  1 (2.75)       uucPa s22 (2.76) Коефицијенти експанзије који се односе на равнотежну функцију расподеле после увођења фактора скалирања могу се изразити следећим једначинама:   00a (2.77)    ua  1 0 (2.78)       1 22 0  scuua (2.79) Докторска дисертација Поглавље 2 Тијана Ђукић 21 Пошто се у овом раду занемарује топлотно струјање, у једначини (2.79) треба узети у обзир да је 1 :     uua  2 0 (2.80) Равнотежна функција расподеле после дискретизације поља брзине има облик:                               2 2 2 42 2 420 2 0 2 0 2 1 2 1 1 2 1 1 ! 1 uuξ uξ uξ s i ss i i i ss i i n i n n n s ii ccc uu cc a nc f    HH (2.81) Ванравнотежни део функције расподеле у дискретизованом облику гласи:        S c f i s ii 2 2 1 H (2.82) Аналогно претходном, може се дискретизовати и израз за спољашњу силу:                         242 1 0 2 1 2 , ! 1 ss ii s i i n i n n s ii ccc at cn uguξgξgξ xgF   H (2.83) Када је у питању једначина за тензор напона, и овде је неопходно узети у обзир фактор скалирања: SIP  22 2 ss cc  (2.84) На исти начин је потребно изменити и израз за динамичку вискозност и једначину која дефинише везу између притиска флуида и густине флуида:  2sc (2.85) 2scp  (2.86) Коначно, BGK Boltzmann-ова једначина дискретизована по пољу брзина гласи:    iiiiii ff f t f F x ξ       01  (2.87) Докторска дисертација Поглавље 2 Тијана Ђукић 22 Макроскопске величине се израчунавају преко функције расподеле, с тим што после дискретизације поља брзина ове величине могу да се израчунају као тежинске суме, док је у случају континуалног поља брзина било неопходно да се врши интеграљење.  i if (2.88)  i ii fξu  1 (2.89)  i ii fe 2 2 1 ξ  (2.90)  i iii fξξP (2.91) 2.3.2 Временска и просторна дискретизација У наставку ће бити спроведен поступак просторне и временске дискретизације, као што је то објашњено у литератури [53,58,59], при чему је заправо најважније дискретизовати конвективни члан једначине (2.87). BGK Boltzmann-ова једначина се може трансформисати на следећи начин:  eqiii ff dt df   1 (2.92) где је   ii eq i Fff  0 (2.93) Ако се једначина (2.92) интеграли, добија се:            t i eq iiiiii dsstsfstsftftttf 0 ,, 1 ,, ξxξxxξx  (2.94) где је t временски корак за интеграцију. Докторска дисертација Поглавље 2 Тијана Ђукић 23 Десна страна ове једначине се може апроксимирати трапезним правилом:                2 0 ,,,, 2 ,, 1 tOtftftttftttf t dsstsfstsf eq iii eq iii t i eq iii       xxξxξx ξxξx   (2.95) У даљем тексту се за временски корак узима 1t , што је погодно за нумеричко решавање једначина и имплементацију на рачунару. Једначина (2.95) је у имплицитном облику, а за нумеричко решавање је погодније користити експлицитни облик. Због тога се уводи смена:         tftftftf eqiiii ,, 2 1 ,, xxxx   (2.96) Ако се изрази (2.95) и (2.96) замене у једначини (2.94) и добијена једначина потом среди, добија се:         tftftftf eqiiiii ,, 1 ,1, xxxξx   (2.97) где је  модификовано време релаксације, које се дефинише: 2 1  (2.98) На овај начин је добијен експлицитни облик једначине (2.95), тако да се у нумеричким симулацијама користи управо if и модификовано време релаксације. У једначини (2.95) је са  2tO  назначено да се чланови вишег реда занемарују. Наравно, овим занемаривањем је направљена одређена грешка и управо да би се кориговала та грешка уведен је корекциони фактор 21 у изразу (2.98). На тај начин је постигнуто да LB метода нема нумеричке нестабилности или грешке. Још једна предност увођења модификованог времена релаксације јесте и то што омогућава примену експлицитних временских корака, дакле бољу и прикладнију временску дискретизацију. Треба нагласити и да је величина eqif у једначини (2.93) дефинисана као функција немодификованог времена релаксације  и физичких величина густине и брзине. Докторска дисертација Поглавље 2 Тијана Ђукић 24 Имајући у виду ову чињеницу, десна страна једначине (2.97) се може трансформисати на следећи начин:                uuxuux , 2 1 1,, 1 ,,, 1 00       iiiiii FftfFftf        (2.99) Сада се може написати једначина која репрезентује целу нумеричку шему LB методе:            uuxxξx , 2 1 1,, 1 ,1, 0     iiiiii Fftftftf        (2.100) Ако се направи осврт на почетак овог поглавља, наглашено је да се поступком дискретизације мора очувати својство BGK једначине да се из ње могу извести Навије- Стоксове једначине. Уколико би се овако дискретизована BGK једначина развила у Taylor-ов ред, добио би се исти облик Навије-Стоксових једначина, који је добијен трансформацијом континуалне BGK једначине, што је показано у литератури [60]. Применом овог поступка може се утврдити колика је грешка занемаривања чланова вишег реда у једначини (2.95) [61]. Такође, применом Chapman-Enskog растављања на дискретизовани облик BGK једначине добио би се поменути облик Навије-Стоксових једначина [55], што је такође поменуто на почетку поглавља 2.2. На овај начин се може доказати да је овакав вид дискретизације исправан. Када се имплементира LB метода, најчешће се једначина (2.100) раздваја на два корака – корак судара (енг. collision step) и корак пропагације (енг. propagation step). Могу се дефинисати две вредности функције расподеле - inif и out if , које представљају вредности дискретизоване функције расподеле пре и после судара. Наведени кораци се могу описати једначинама: Корак судара:           iiiniiniouti Fftftftf           2 1 1,, 1 ,, 0 uxxx (2.101) Корак пропагације:    tftf outiiini ,1, xξx  (2.102) Ова два корака се понављају у итерацијама, при чему сваки од ова два корака мора да се примени на све чворове мреже, пре него што почне следећи корак. Алгоритам рада једног програма базираног на LB методи приказан је на слици 2.1. Докторска дисертација Поглавље 2 Тијана Ђукић 25 Слика 2.1 - Алгоритам рада програма базираног на LB методи 2.4 Израчунавање макроскопских величина коришћењем функције расподеле Макроскопске величине, као што су густина, притисак, брзина, смичући напон израчунавају се коришћењем израчунатих вредности компоненти функције расподеле. Изрази помоћу којих се неке од ових величина израчунавају већ су написани у овом поглављу, али ће у наставку бити наведени коначни изрази, који ће бити коришћени у симулацијама. Због увођења модификованог времена релаксације (које је дефинисано једначином (2.98)), потребно је променити и изразе за израчунавање макроскопских величина симулације, да би се обезбедила конзистентност модела. Докторска дисертација Поглавље 2 Тијана Ђукић 26 Уместо  треба написати 2 1  у изразу којим се дефинише вискозност:        2 12  sc (2.103) Такође, услед увођења смене (2.96), треба проверити да ли долази до промена у изразима за макроскопске величине. Када је у питању густина, нема никаквих промена [41,53,54]:  i if (2.104) Међутим, израз за брзину сада гласи: 2 1 g uξu   i ii f  (2.105) Из ове једнакости је јасно да брзина која се у симулацијама израчунава и која је означена са u , не представља „физичку брзину“ струјања. Да би се одредила физичка брзина (која је једна од најважнијих карактеристика струјања на макроскопском нивоу), потребно је извршити следеће прерачунавање: 2 g uu  (2.106) Наравно, у случају да не постоји утицај спољашњих сила (када је 0g ) брзина је иста без обзира на коришћење поменуте смене у симулацији ( uu  ). Притисак се израчунава преко густине флуида у посматраној тачки, применом једнакости: 2scp  (2.107) Потребно је објаснити и како се у оквиру симулације применом LB методе може одредити смичући напон. Наиме, тензор напона је по дефиницији једнак збиру притиска и напону који се јавља као последица вискозних сила: ' ijijij p   (2.108) где је ij Kronecker-ов делта симбол. Докторска дисертација Поглавље 2 Тијана Ђукић 27 Применом Cauchy-еве формуле може се израчунати вектор напона у некој тачки:   jiji nT  n (2.109) где је jn вектор нормале на зид. Смичући напон се може израчунати као разлика између укупног напона и његове пројекције на нормалу на зид:   ikjkjjiji nnnn   (2.110) Заменом израза за тензор напона (2.108) у једначину (2.110), узимајући у обзир да је 1kjjk nn , добија се:   ikjkjjiji nnnn ''   (2.111) Сада је потребно одредити вредност напона који потиче од вискозних сила. За нестишљив, Њутновски флуид, овај тензор напона линеарно зависи од тензора деформације и може се изразити на следећи начин: ijij S 2 '  (2.112) Коначни израз за смичући напон гласи:    ikjijjiji nnnSnS   2 (2.113) Тензор напона се у LB методи може изразити преко функције расподеле и ова веза је дефинисана једначином (2.91). Већ је дефинисана веза између ванравнотежног дела тензора напона и тензора деформације и она се може записати следећом једначином:   SP 21 2 sc (2.114) Заменом једначина (2.91) и (2.114) у једначину (2.113), добија се израз којим се на основу функције расподеле и познатих вектора нормале одређује смичући напон:     lijljikjkjj s i nnnff c      0 2 (2.115) Тијана Ђукић 28 3 Моделирање понашања деформабилног тела Деформабилна тела која се разматрају у овом раду имају структуру као црвена крвна зрнца и веома су подложна великим деформацијама. У унутрашњости црвених крвних зрнаца налази се нестишљив Њутновски флуид (хемоглобин), који је окружен двослојном мембраном липида, која садржи и танак скелет сачињен од протеина [3]. Управо због такве структуре еритроцита, и услед удруженог дејства флуида који се налази у унутрашњости и вискоеластичне мембране, ћелије имају велику способност деформације [62,63]. У овом раду се сматра да је мембрана деформабилног тела занемарљиве дебљине и да је сачињена од одређеног броја тачака које су међусобно повезане. Због тога је извршена таква дискретизација да је цео домен солида подељен на одређени број троуглова и на тај начин је креирана мрежа. Тада је могуће за сваки елемент и за сваки чвор мреже посебно одредити силу отпора којом се тај чвор супротставља дефинисаној спољашњој деформацији. Постоје четири параметра који утичу на понашање деформабилне честице, а то су: запремина унутар мембране, површина спољашње мембране, површинско напрезање мембране и савијање мембране. У наставку поглавља биће описано како се моделира реакција деформабилне честице на промену сваког од четири поменута параметра. 3.1 Моделирање реакције на површинско напрезање мембране Постоје два начина моделирања површинског напрезања мембране. Један је да се мембрана третира као скуп чворова који су међусобно повезани опругама. Овај приступ су користили Dupin и други [25]. Тада свако померање чворова изазива силу у опругама и на тај начин се одређује реакција деформабилне честице на овај вид деформације. Овај приступ је једноставнији и брже се могу нумерички одредити вредности сила у чворовима. Ако се са 12L и 12 0L означе тренутно растојање између чворова 1j и 2j и почетно растојање између посматраних чворова, респективно, као што је то скицирано на слици Докторска дисертација Поглавље 3 Тијана Ђукић 29 3.1, онда се услед померања чворова јавља сила која се може одредити на следећи начин: 12 12 0 12 0 12 21 lFF L LL K SSj S j   (3.1) где је са SK означен параметар отпора, који се дефинише за одговарајући тип деформабилне честице, а са 12l је означен јединични вектор чији се правац поклапа са правом која спаја чворове 1j и 2j и овај вектор је такође скициран на слици 3.1. Слика 3.1 - Скициране величине које су потребне за дефинисање сила реакције на површинско напрезање Други приступ моделирању површинског напрезања мембране је да се мембрана посматра као хипереластичан материјал, при чему се однос напона и деформације може одредити на основу густине деформационе енергије. Развијено је неколико хипереластичних материјалних модела. За моделирање везе између напона и деформације код црвених крвних зрнаца који су главна тема овог рада, користе се Mooney-Rivlin-ов материјални модел, као и neo-Hook-ов материјални модел [18,21,64]. Међутим, Skalak и други [37] су анализирали поменуте моделе и дошли су до закључка да они не могу са довољном тачношћу да симулирају понашање црвених крвних зрнаца, па су предложили свој материјални модел за моделирање понашања мембране црвеног крвног зрнца, који ће бити коришћен у овом раду. Густина деформационе енергије се у Skalak-овом моделу деформабилне мембране дефинише као функција две инваријанте Cauchy-Green-овог тензора деформације:   222121 ' 12 '2'2' 12 I k III k W sS  (3.2) при чему су са '1I и '2I означене измењене инваријанте као што је предложено у [37]. Са sk и k означени су параметри отпора на смичућу и површинску деформацију у Докторска дисертација Поглавље 3 Тијана Ђукић 30 равни мембране, респективно. Ови параметри се дефинишу за одговарајући тип деформабилне честице. Пре него што се дефинишу инваријанте и измењене инваријанте из једначине (3.2), потребно је дефинисати још две величине помоћу којих се инваријанте дефинишу, а то су тензор градијента деформације и леви Cauchy-Green-ов тензор деформације. Када се посматра неко деформабилно тело, у почетном тренутку ово тело се налази у почетној конфигурацији која се означава са B0 и тада одабрана материјална тачка P има координате x0 . У тренутно посматраном временском тренутку t , тело се налази у тренутној конфигурацији Bt и одабрана материјална тачка има координате xt . Веза између координата у тренутно посматраном временском тренутку и почетном временском тренутку јесте градијент деформације, који се може написати као: j k t kj t x x F 00    (3.3) Ако се померање материјалне тачке дефинише као: xxu 0t (3.4) и ако се овај израз примени у једначини (3.3), добија се: x u IF 00   t (3.5) где је I јединични тензор. На сличан начин се може дефинисати градијент деформације, који се посматра у конфигурацији Bt , а и мери се у истој конфигурацији B t (израз је дат у индексној нотацији): j t i kjkj t t x u F    (3.6) где је kj Kronecker-ов делта симбол. Леви Cauchy-Green-ов тензор деформације се сада може дефинисати: T FFB  (3.7) или у индексној нотацији: Докторска дисертација Поглавље 3 Тијана Ђукић 31 jmimij FFB  (3.8) Главне инваријанте тензора су величине које не зависе од координатних трансформација и могу се изразити у функцији Cauchy-Green-овог тензора деформације B , на следећи начин [65]: nnBBBBI  3322111 trB (3.9)   222 trtr 2 1 BB I (3.10) Такође, инваријанте се могу изразити у функцији главних издужења 1 , 2 и 3 : 2 3 2 2 2 11  I (3.11) 2 1 2 3 2 3 2 2 2 2 2 12  I (3.12) Уколико се посматра танка деформабилна мембрана, у локалном координатном систему се може сматрати да је у питању напрезање у једној равни, па се због тога може занемарити издужење у правцу z осе, односно 03  . Измењени облици инваријанти, као што су предложили Skalak и други [37], гласе: 22' 1 2 2 2 11  II  (3.13) 11' 2 2 2 2 12  II  (3.14) На основу дефиниција за измењене инваријанте, могуће је одредити везу између главних инваријанти и густине деформационе енергије код Skalak-овог модела:    12 12 222 12 2 2 221 2 1  II k III k W sS  (3.15) Сада је потребно одредити везу између напона и деформације, користећи једначину густине деформационе енергије (3.15). Cauchy-ев тензор напона σ се може изразити на следећи начин [66]: Докторска дисертација Поглавље 3 Тијана Ђукић 32 jmimij FT J 1  (3.16) где су са imT означене компоненте Piola-Kirchhoff-овог тензора напона прве врсте, а J је детерминанта градијента деформације:  FdetJ (3.17) Piola-Kirchhoff-ов тензор напона прве врсте се може представити као парцијални извод густине деформационе енергије по градијенту деформације: im S im F W T    (3.18) Заменом једначине (3.18) у једначину (3.16) добија се: jm im S ij F F W J    1  (3.19) Парцијални извод из једначине (3.19) се може израчунати посредним диференцирањем: im S im S im S im S F J J W F I I W F I I W F W                  2 2 1 1 (3.20) Парцијални извод детерминанте по градијенту деформације једнак је: 1   mi im JF F J (3.21) Парцијални извод прве инваријанте по градијенту деформације једнак је:   im im np np im npnp im F F F F F FF F I 221          (3.22) Парцијални извод друге инваријанте по градијенту деформације једнак је: Докторска дисертација Поглавље 3 Тијана Ђукић 33        kmikimpn im np im pn np im im pnnp imimim FBFIB F B F B B F I I F BB FFF I 222 2 1 tr tr2 2 1 trtr 2 1 1 1 1 22 2                                     B B BB (3.23) Заменом израза (3.21)-(3.23) у израз (3.20), а затим и у једначину (3.19), добија се:   jmmi S kmikim S im S ij FJF J W FBFI I W F I W J                 11 21 222 1  (3.24) Ако се узме у обзир једначина (3.8) и једнакост 1 jkijik FF (3.25) једначина за Cauchy-ев тензор напона постаје:   kjik S ij S ij SS ij S kjikij S ij S ij BB I W JJ W B I W I I W J J J W BBBI I W B I W J 22 1 1 1 21 22 222 1                                   (3.26) На основу једначине (3.15) могуће је израчунати парцијалне изводе густине деформационе енергије по детерминанти и инваријантама:  22 12 1 1    I k I W s S ;  22 126 2 2    I kk I W s S  ; 0   J W S (3.27) Коначни израз којим се повезују напон и деформација гласи:                        kjik s ij s ij BBI kk BII kk J 1 33 1 33 1 221  (3.28) Докторска дисертација Поглавље 3 Тијана Ђукић 34 Ова изведена релација важи за континуум. Да би се она могла користити у нумеричким симулацијама, потребно је да се примени метода коначних елемената (МКЕ). Наиме, цео домен, односно мембрана деформабилне честице се дискретизује – дели на поддомене, који се називају коначни елементи. Овом дискретизацијом је омогућена нумеричка симулација реакције континуума на спољашње утицаје, у овом случају симулација реакције мембране деформабилне честице на силе којима околни и унутрашњи флуид делују на ту мембрану. Све физичке величине које су обухваћене симулацијом, у овом случају то су унутрашње силе које су изазване деформацијом мембране, методом коначних елемената биће израчунате у дискретном облику, односно у свим тачкама дискретизоване мреже. Ове тачке у којима се рачунају тражене величине се називају чворне тачке мреже коначних елемената или скраћено чворови. Дакле, у оквиру решавања применом МКЕ потребно је на основу задатих померања чворних тачака одредити унутрашње силе које се јављају као последица дефинисаних деформација. На основу померања чворних тачака могуће је одредити градијент деформације, тензор деформације, као и напон, применом претходно написаних једначина. Сада треба дефинисати једначину која повезује напон и унутрашње силе. Ова једначина се може дефинисати применом принципа виртуалног рада [34,67]. Ако се са u означе виртуална померања, а са e означе виртуалне деформације, тада се виртуални рад унутрашњих сила SuF може представити као рад напона σ на виртуалним деформацијама:  V TS u T u dVA σeFu  (3.29) Треба напоменути да принцип виртуалног рада важи и за велика и за мала померања, велике и мале деформације, а применљив је и за било који тип конститутивних релација [34] (за еластично тело, еласто-пластично тело, круто тело, као и за хипереластично тело, које се разматра у овом раду). Већ је напоменуто да се све величине у симулацијама применом МКЕ израчунавају само у чворним тачкама. У свим осталим тачкама дискретизованог домена, величине се могу израчунати интерполацијом величина у чворним тачкама. Ако се координате k - тог чвора означе са k iX , где i означава индекс координатне осе ( 1i одговара x оси, 2i одговара y оси, а 3i одговара z оси), онда се координата произвољне тачке унутар посматраног коначног елемента може одредити интерполацијом:    N k k iki Xhx 1 (3.30) где је N укупан број чворова коначног елемента (на пример 4N за изопараметарски дводимензионални коначни елемент, 8N за изопараметарски тродимензионални коначни елемент итд.), а kh су интерполационе функције које се дефинишу за одговарајући тип коначног елемента. Докторска дисертација Поглавље 3 Тијана Ђукић 35 Слично једначини (3.30), на основу померања чворних тачака могуће је одредити померање произвољне тачке посматраног елемента:    N k k iki Uhu 1 (3.31) У оквиру симулација применом МКЕ најпре се све тражене величине одређују у тзв. природном (локалном) координатном систему, који је везан за сваки елемент посебно, затим се једначине за сваки елемент посебно групишу у једну целину, у систем једначина целог домена, које се потом решавају и на крају се добијене величине трансформишу назад у глобални координатни систем. Природни (локални) координатни систем се дефинише координатним осама r , s и t , које су међусобно ортогоналне, као и осе глобалног координатног система, а средиште локалног координатног система се везује за центар посматраног коначног елемента. Заправо, елемент се из простора дефинисаног x , y и z координатама пресликава на елемент чије координате чворова задовољавају услове: 11  ir , 11  is , 11  it . Слика 3.2 - Троугаони коначни елемент у глобалном Декартовом и природном координатном систему Табела 3.1 - Координате чворних тачака у природном координатном систему Редни број чвора ir is 1 1 0 2 0 1 Докторска дисертација Поглавље 3 Тијана Ђукић 36 3 0 0 У уводном делу поглавља речено је да је поред утицаја површинског напрезања, потребно разматрати и друге утицаје на деформабилну честицу, који ће бити разматрани у наставку поглавља. Али, како је за одређивање реакција на друге поменуте утицаје најпогодније користити дискретизацију на троугаоне елементе, и за симулације применом МКЕ ће бити коришћени троугаони коначни елементи. Ово није један од „класичних“ елемената који се користи у анализама применом МКЕ, али је овде одабран због коришћења јединствене дискретизоване мреже за одређивање свих реакција. На слици 3.2 је приказан троугаони коначни елемент, са означеним осама природног координатног система. Троугаони коначни елемент је дводимензионални елемент, тако да се трећа координата у локалном координатном систему не разматра (сматра се да је 0t ). У табели 3.1 наведене су координате чворова коначног елемента у природном координатном систему (такође је примењена локална нумерација чворова). За овај тип коначног елемента интерполационе функције се дефинишу на следећи начин: rh 1 sh 2 srh 13 (3.32) Вектор деформације e (који има 6 компоненти: 123123332211 ,,,,, eeeeee ) се у одређеној материјалној тачки унутар коначног елемента може одредити помоћу једначине: UBe h (3.33) где је hB матрица која садржи изводе интерполационих функција. Ова матрица, исто као и интерполационе функције, има другачији облик у зависности од типа коначног елемента. За троугаони коначни елемент, за k -ти чвор, матрица извода интерполационих функција има облик:            1, 2, 2, 1, 0 0 k k k k k h h h h h B (3.34) при чему је извод интерполационе функције по координати j k jk x h h   , . Докторска дисертација Поглавље 3 Тијана Ђукић 37 Наравно, треба напоменути како се израчунавају изводи интерполационих функција по глобалним координатама. Пошто су интерполационе функције дефинисане у функцији локалних координата, може се применити посредно диференцирање: j k j k j k x s s h x r r h x h             (3.35) Уводи се Јакобијан трансформације између глобалног и локалног система:                            s y s x r y r x r x J (3.36) а може се дефинисати и инверзни Јакобијан:                              y s y r x s x r x r J 1 (3.37) Сада се изводи из једначине (3.35) могу написати: s h J r h J x h k j k j j k          1 2 1 1 (3.38) Изводи интерполационих функција по природним координатама се лако могу одредити из једначине (3.32). Елементи матрице Јакобијана се могу одредити заменом једначине (3.30) у једначину (3.36): k j N k i k ij X r h J      1 (3.39) Ако се израз за вектор деформације (3.33) замени у једначини виртуелног рада (3.29), добија се: Докторска дисертација Поглавље 3 Тијана Ђукић 38 S u T V T h T V T u dVdVA FUσBUσe    (3.40) Следи да се унутрашње силе могу израчунати на основу напона применом релације:  V T h S u dVσBF (3.41) Да би се на основу напона израчунале унутрашње силе, потребно је одредити напон. У овом раду је коришћен Skalak-ов модел мембране деформабилне честице, тако да се напон израчунава применом већ изведене једначине (3.28). Треба само напоменути да се након спроведене дискретизације на коначне елементе градијент деформације на нивоу елемента може израчунати на следећи начин:      N k j kk iijij x h UF 1  (3.42) Слика 3.3 - Троугаони коначни елемент у глобалном и локалном Декартовом координатном систему Међутим, у Skalak-овом моделу мембране деформабилне честице сматра се да је мембрана занемарљиво мале дебљине и да се може сматрати да се напрезање јавља само у равни мембране. Са друге стране, померања чворова коначних елемената су дефинисана у глобалном Декартовом координатном систему и могу имати компоненте у правцу све 3 осе. Због тога је неопходно пре било каквих израчунавања извршити трансформацију вектора померања из глобалног Декартовог у локални Декартов координатни систем. На слици 3.3 приказан је један троугаони коначни елемент и означене су осе глобалног и локалног Декартовог система. Да би се дефинисале координатне осе локалног система, најпре се одреде два вектора који леже у равни Докторска дисертација Поглавље 3 Тијана Ђукић 39 елемента. Правац вектора 1i поклапа се са једном страницом троугла, овај вектор је јединични вектор и користи се да се дефинише правац координатне осе x локалног система. За правац вектора 2i се може узети да се поклапа са другом страницом троугла, као што је приказано на слици 3.3. Векторски производ ова два вектора дефинише вектор нормале на раван елемента 21 21 1 ii ii k    (3.43) и овај вектор се користи да се дефинише правац координатне осе z локалног система. Ознака p се користи да се означи интензитет вектора p . Сада је лако одредити вектор којим се дефинише правац координатне осе y , поново применом векторског производа: 11 11 1 ik ik j    (3.44) Јединични вектори 1i , 1j и 1k дефинишу осе локалног Декартовог координатног система. Сваки од ових вектора се може представити на следећи начин: kjii 1111 coscoscos   (3.45) при чему су 1 , 1 и 1 углови између јединичних вектора 1i и i , 1i и j , 1i и k , респективно. Могу се увести ознаке: 11 cosl ; 11 cosm ; 11 cosn (3.46) Треба напоменути да важи једнакост: 121 2 1 2 1  nml (3.47) Исте једнакости важе и за друга два јединична вектора. Ако се са kU означе померања k -тог чвора у локалном координатном систему (померање у правцу осе z се може занемарити пошто је у оквиру материјалног модела уведена претпоставка да се ради о напрезању у равни мембране), онда постоји веза између померања у локалном и глобалном систему: Докторска дисертација Поглавље 3 Тијана Ђукић 40                                                      N N N N N U U U U U U nml nml nml nml U U U U 3 2 1 1 3 1 2 1 1 222 111 222 111 2 1 1 2 1 1 000 000 000 000 000 000 000 000  (3.48) или у матричном облику, за један елемент: TUU  (3.49) где је са T означена матрица трансформације. Лако се може утврдити да важи једнакост: ITT T (3.50) при чему је I јединична матрица димензија 6x6. За конкретан случај који се разматра у овом раду, нумерички је најефикасније да се све величине израчунавају у локалном систему, а да се коначни резултат (вредности сила у чворним тачкама) на крају трансформишу у глобални систем. Трансформација (3.49) важи и за вектор сила у чворовима коначног елемента: S u S u TFF  (3.51) или S u TS u FTF  (3.52) Вектор сила у чворовима у локалном систему има по две компоненте за сваки чвор, укупно 6 компоненти за један троугаони елемент, док вектор у глобалном систему има по три компоненте за сваки чвор, укупно 9 компоненти за један троугаони елемент.  NSuNSuSuSuTSu FFFF 211211 F  NSuNSuNSuSuSuSuTSu FFFFFF 321131211 F (3.53) Интеграли на нивоу коначног елемента се нумерички израчунавају применом Gauss-ове интеграције. У овом раду се разматра проблем напрезања у равни мембране, па ће због Докторска дисертација Поглавље 3 Тијана Ђукић 41 тога и интеграли бити рачунати у дводимензионалном домену. За тај случај се применом Gauss-ове интеграције добија:          i jn i n j jiji wwsrfdrdssrf 1 1 1 1 1 1 ,, (3.54) где је in и jn број интеграционих тачака у правцу оса природног координатног система, iw су тежински коефицијенти, а са  ji sr , је дефинисана позиција Gauss-ове тачке. У табели 3.2 су дате координате Gauss-ових тачака и одговарајући тежински коефицијенти, за троугаони коначни елемент. Табела 3.2 – Gauss-ова нумеричка интеграција за троугаоне коначне елементе Ред интеграције Интеграционе тачке r координата s координата Тежински коефицијенти 3 тачке 16666667.01 r 66666667.02 r 13 rr  11 rs  12 rs  23 rs  33333333.01 w 12 ww  13 ww  7 тачака 10128651.01 r 79742698.02 r 13 rr  47014206.04 r 45 rr  05971587.06 r 33333333.07 r 11 rs  12 rs  23 rs  64 rs  45 rs  46 rs  77 rs  12593918.01 w 12 ww  13 ww  13239415.04 w 45 ww  46 ww  225.07 w Докторска дисертација Поглавље 3 Тијана Ђукић 42 3.2 Моделирање реакције на промену запремине Деформабилна честица тежи томе да унутрашња запремина остане константна током времена, без обзира на деформације којима је честица изложена. Унутрашња запремина се може израчунати на следећи начин: i i iiAV hn  (3.55) где је сумирање извршено по свим елементима мреже Mi ,,2,1  , а M је број елемената. У једначини (3.55) са ih је означен вектор који дефинише удаљеност тежишта C посматраног елемента од одабране координатне равни (у овом раду је одабрано да то буде раван x-y, пошто је то раван која се поклапа са дном флуидног домена, па ће сигурно све честице бити увек изнад ове равни), са in означен је јединични вектор нормале посматраног елемента, а iA је површина посматраног елемента. Све ове величине су означене на слици 3.4. Слика 3.4 - Скициране величине које су потребне за израчунавање запремине унутар мембране Уколико се са V и V0 означе тренутна запремина деформабилне честице после одређене деформације и почетна запремина, респективно, онда таква промена Докторска дисертација Поглавље 3 Тијана Ђукић 43 запремине изазива силу у сваком елементу, као што је наведено у [25], и ова сила је једнака: ii VV i A V VV K nF 0 0  (3.56) где је са VK означен параметар отпора, који се дефинише за одговарајући тип деформабилне честице. У овом раду је разматрана дискретизација на троугаоне елементе, па се због тога сила израчуната за елемент равномерно распоређује на сва 3 чвора посматраног елемента: V i V i V i V i FFFF 3 1 321  (3.57) 3.3 Моделирање реакције на промену површине спољашње мембране Уколико се са iA и iA 0 означе тренутна површина i -тог елемента спољашње мембране после одређене деформације и почетна површина i -тог елемента мембране, респективно, онда таква промена површине изазива силу у сваком чвору сваког елемента, као што је наведено у [25] и ова сила је једнака: j i iiAA ij A AA K wF 0  (3.58) где је са i означен индекс елемента, а j је локални индекс чвора у оквиру тог елемента  3,2,1j , док је AK параметар отпора који се дефинише за сваки тип деформабилне честице. Са jw означен је јединични вектор чији се правац поклапа са правом која спаја тежиште посматраног елемента са посматраним чвором елемента. Вектори jw су означени на слици 3.5. Докторска дисертација Поглавље 3 Тијана Ђукић 44 Слика 3.5 - Скициране величине које су потребне за дефинисање сила реакције на промену површине Међутим, потребно је узети у обзир и утицај који промена укупне површине мембране има на сваки елемент посебно. Уколико се разматра случај када је деформабилна честица црвено крвно зрнце, онда укупна промена површине има много већи утицај на силу реакције на деформацију, од локалне промене површине [25]. Овај услов се у нумеричким симулацијама узима у обзир кроз следећу једначину: j At i iiAlA ij A AA K A AA K wF           00 (3.59) где горњи индекс l означава локалну промену површине, а горњи индекс t означава укупну (тоталну) промену површине. Услов који је претходно поменут се узима у обзир тако што се приликом дефинисања параметара AlK и AtK води рачуна о томе да важи: AlAt KK  (3.60) 3.4 Моделирање реакције на савијање мембране Да би се дефинисало савијање мембране, посматра се промена угла између два елемента. Ако се са 21ii  и 21 0 ii означе тренутни и почетни угао између елемената 1i и 2i , респективно, онда се услед промене угла јавља сила у сва три чвора оба елемента Докторска дисертација Поглавље 3 Тијана Ђукић 45 [25], при чему је сила у чвору који није заједнички за ова два елемента (чвор обележен са 31N на слици 3.6) једнака: 1 21 2121 3 0 0 i ii iiiiBB j K nF     (3.61) док је сила у заједничним чворовима једнака: B j B j B j 32112 2 1 FFF  (3.62) Слика 3.6 - Скициране величине које су потребне за дефинисање сила реакције на савијање мембране Укупна сила којом се солид супротставља деформацији се израчунава у сваком чвору, као сума сила израчунатих за сва четири параметра деформације.   Bi A i V i S ii t FFFFF  (3.63) Тијана Ђукић 46 4 Моделирање солид-флуид интеракције Постоје одређени проблеми код којих је неопходно истовремено решавање два или више физичких система који међусобно интерагују. Пример таквог проблема постоји и у области моделирања струјања флуида, када се унутар домена флуида налазе честице, односно солиди. При томе флуид делује одређеним силама на солид и проузрокује његово кретање и деформисање, а важи и обрнуто – солид се супротставља утицају флуида и на тај начин утиче на струјање. То практично значи да солид и флуид формирају један спрегнути механички систем. Да би се моделирао један овакав систем потребно је разматрати оба домена истовремено и моделирати солид-флуид интеракцију. Методе које се користе приликом моделирања понашања солида и флуида су веома различите, тако да је потребно користити неки приступ који ће омогућити њихово спрезање и уклапање у једну целину. За моделирање солид-флуид интеракције користе се два приступа: слабо и јако спрезање. У случају слабог спрезања солид и флуид се практично моделирају појединачно, а сви потребни параметри добијени из солвера за један домен се прослеђују солверу за други домен. Посебан спољашњи програм управља решавањем оба домена и обезбеђује правилну размену информација. У случају јаког спрезања оба домена се моделирају истовремено, као да је у питању један механички систем. Формира се систем једначина који се решава у сваком временском кораку и тако се обезбеђује да се све величине (како за солид, тако и за флуид) мењају истовремено. Слабо спрезање је једноставније за имплементацију, али један од главних недостатака овог приступа је чињеница да је неопходно користити исти временски корак за нумеричко решавање оба домена, што није увек могуће пошто су карактеристике солида и флуида различите. Јако спрезање се примењује када је потребно прецизно моделирање кретања солида у флуиду. Главни недостатак овог другог приступа је што је солвер знатно комплекснији и извршавање оваквог програма је спорије. Јако спрезање је коришћено у овом раду због веће прецизности и могућности примене различитог временског корака и различите дискретизације домена солида и флуида. Суштина принципа јаког спрезања приказана је блок дијаграмом на слици 4.1. Моделирање струјања флуида се спроводи применом CFD солвера (енг. computational fluid dynamics – CFD), а моделирање понашања солида се спроводи применом CSD солвера (енг. computational solid dynamics – CSD). Да би се повезали домени солида и флуида, потребно је да се одређене информације добијене из домена флуида користе за прорачун у домену солида и обрнуто. CFD солверу се прослеђују подаци о тренутној геометрији додирне површине солида и флуида. Са друге стране, CSD солверу се прослеђују оптерећења (у овом случају силе) којима флуид делује на додирне површине. На тај начин се моделира солид-флуид интеракција. Докторска дисертација Поглавље 4 Тијана Ђукић 47 Слика 4.1 – Размена информација приликом моделирања солид-флуид интеракције Постоји неколико приступа којима се нумерички моделира кретање солида у домену флуида. Неки од приступа су коришћење координатне трансформације [68], коришћење покретне мреже коначних елемената [69], примена тзв. ALE формулације (енг. Arbitrary Lagrangian Eulerian) [70], примена Immersed boundary методе (IBM) [39]. Сваки од ових приступа има своје предности и недостатке. У овом раду је коришћена Immersed boundary метода, због тога што се код ове методе систем једначина за флуид не мења, већ се дејство солида на флуид остварује преко спољашње силе. Такође, могуће је користити мреже различите густине за дискретизацију оба домена, пошто се спрезање врши интерполацијом у околини сваке тачке оба домена. У наставку овог поглавља биће описани детаљи ове методе. 4.1 Immersed boundary метода Immersed boundary методу (IBM), прво је осмислио и представио Peskin [38]. Код ове методе се за представљање домена флуида користи фиксна Декартова мрежа, која се састоји од Euler-ових тачака. Овакав опис флуида је заправо еквивалентан репрезентацији домена флуида која се користи у LB методи, па је јасно да се IBM може применити када се струјање флуида моделира LB методом. Са друге стране, IBM посматра солид као тело уроњено у домен флуида, отуда и назив методе (енг. immersed boundary значи уроњена граница). Солид је дакле издвојени део флуида, са границом која се представља скупом Lagrange-ових тачака. Основна идеја је да се физичка граница између два домена третира као лако деформабилна граница, велике крутости [40]. Флуид делује на солид, односно на граничну површину преко силе, која деформише границу између два домена. Ова метода се успешно примењује и у случају Докторска дисертација Поглавље 4 Тијана Ђукић 48 чврстих тела уроњених у флуид [39,40,41], а и у случају деформабилних тела [36,42]. У зависности од тога да ли је солид деформабилан или не, зависи каква ће реакција солида на деформацију бити. Услед овако насталих деформација, у солиду се јавља унутрашња сила отпора, којом ће солид деловати на околни флуид. Цео домен флуида се решава применом Навије-Стоксових једначина, узимајући у обзир дејство уроњеног солида, које се у домен флуида уводи преко спољашње силе. Гранична површина солида уроњеног у флуид се представља скупом Lagrange-ових тачака, при чему се координате l -те тачке означавају са  tlBX , Gl ,,2,1  , а G је број тачака на граничној површини. Све величине које су значајне за симулацију интеракције између солида и флуида се рачунају интерполацијом преко тачака на граничној површини. Пошто постоји могућност различите дискретизације домена флуида и солида, проблем интерполације са различитом дискретизацијом се решава тако што се за једну граничну тачку која припада солиду разматра утицај већег броја тачака из фиксиране мреже домена флуида, као што је приказано на слици 4.2. Слика 4.2 - Међусобни утицај тачака домена флуида и солида Интерполација се спроводи применом Dirac-ове делта функције, која је апроксимирана на следећи начин:         lBijlBijlBijijB YyXxDt   XxXx (4.1) Вредност  r је дефинисана у литератури [38]:                       2,0 2, 2 cos1 4 1 r r r hr   (4.2) где је са h означено растојање између 2 тачке фиксне мреже (у овом раду, пошто је у питању LB метода, 1h ). Докторска дисертација Поглавље 4 Тијана Ђукић 49 Овако дефинисана Dirac-ова делта функција обезбеђује да интерполирана вредност буде максимална у чворовима две мреже који се апсолутно поклапају и да се утицај интерполиране величине смањује са повећањем растојања између посматраних чворова, као што се може видети на слици 4.3. Слика 4.3 – Дијаграм промене вредности Dirac-ове функције у зависности од растојања између чворова У оквиру симулације струјања флуида ефекат солида који је уроњен у флуид се симулира преко граничне површине, односно преко спољашње запреминске силе која делује на флуид у околини граничне површине. Та сила се може изразити на следећи начин:          dsttt BXxFxg , (4.3) где је са  tF сила којом се солид супротставља дејству флуида. Израчунавање ове силе је детаљно објашњено у поглављу 3. Ако се примени интерполација применом Dirac-ове делта функције, постиже се да се утицај силе у једној граничној Lagrange-овој тачки солида преноси на више тачака Euler-ове мреже флуида (што је илустровано на слици 4.2) и добија се:          m l l Bijijlij tDtt 1 , XxFxg (4.4) Докторска дисертација Поглавље 4 Тијана Ђукић 50 Основна предност IB методе је што се ефекат солида уроњеног у флуид симулира само преко дејства спољашње силе, а није потребно мењати једначине које важе за флуид, у овом случају основну једначину LB методе (једначину (2.100)), која је овде поновљена:            uuxxξx , 2 1 1,, 1 ,1, 0     iiiiii Fftftftf        (4.5) У овој једначини треба обратити пажњу на израз за спољашњу силу iF . Да би се исправно имплементирало дејство спољашње силе у LB методи, мора се водити рачуна да се нешто не поремети у систему једначина изведеном у поглављу 2. Детаљну анализу поступка укључивања спољашње силе спровели су Guo и други [71] и на основу њихових закључака је имплементирана сила и у овом раду. Израз за силу има следећи облик (то је заправо једначина (2.83) из поглавља 2, која је овде поновљена ради лакшег праћења извођења):    t cc s ii s i ii ,42 xg ξuξuξ F            (4.6) Такође је важно нагласити да се дејство спољашње силе мора узети у обзир и приликом рачунања физичког поља брзине струјања, тако да се физичка брзина израчунава применом једначине: 2 g ξu  i ii f (4.7) До сада је објашњено како се моделира ефекат солида у флуиду. Флуид такође делује на солид, а утицај флуида на солид се моделира тако што се сматра да нема проклизавања на граници између ова два домена. Да би био задовољен услов да на граници између солида и флуида нема проклизавања, брзина у домену флуида у граничним тачкама мора бити једнака брзини солида у тој истој граничној тачки. Брзина у Lagrange-овим тачкама се интерполира на основу околних Euler-ових тачака које се налазе у флуиду (што је илустровано на слици 4.2), што се може представити следећом једначином:             S BB B dtttt t t xXxxuXu X ,, (4.8) Уколико се примени већ дефинисана Dirac-ова делта функција, брзина у свакој Lagrange-овој тачки се може израчунати:        ji l Bijijij l B l B Dtt , ,, XxxuXu (4.9) Докторска дисертација Поглавље 4 Тијана Ђукић 51 Применом Euler-ове интеграције на једначину (4.9), добија се нова позиција граничних тачака солида: tlB l B tl B tt  uXX (4.10) На основу ових података одређује се деформација солида и тада се примењује поступак описан у поглављу 3, да би се одредиле силе реакције које је нова деформација изазвала у солиду. Те силе се затим преносе на флуид, а флуид потом изазива нове деформације солида и то је поступак који се понавља, док се не достигне жељени услов или задати број итерација. На слици 4.4 је приказан алгоритам поступка спрезања домена солида и флуида применом Immersed boundary методе. Слика 4.4 - Алгоритам моделирања солид-флуид интеракције применом Immersed boundary методе Тијана Ђукић 52 5 Моделирање струјања вишекомпонентног флуида Једна од предности LB методе наспрам класичних метода рачунске динамике флуида је једноставно моделирање струјања вишекомпонентног флуида. Постоје многи проблеми који се моделирају применом нумеричких симулација вишефазних флуида или више различитих флуида. Неки од оваквих примера су симулације струјања воде и водене паре или воде и уља. Оваква струјања су комплексна зато што је потребно моделирати и интеракцију између компоненти на микроскопском нивоу. Ако се разматра на пример струјање два флуида, у класичним методама потребно је написати диференцијалне једначине за обе компоненте, као и једначину која симулира интеракцију између компоненти. Имплементација и нумеричко решавање једног оваквог система једначина је сложено и захтевно са софтверске тачке гледишта. То није случај са LB методом, пошто се код LB методе разматра кретање честица на микроскопском нивоу и на основу тога се одређују тражене макроскопске величине. Због тога се обе компоненте могу симулирати применом истих једначина и једноставно се одређеном скупу честица додељују карактеристике једног, односно другог флуида. Потребно је само да се у систем једначина дода сила интеракције између две компоненте. Сва израчунавања се врше за све честице истовремено, што значи упоредо за обе компоненте. У наставку овог поглавља биће детаљније описани поступци моделирања струјања вишекомпонентног флуида. 5.1 Моделирање струјања вишекомпонентног флуида применом Shan-Chen-овог модела У овом раду је коришћен Shan-Chen-ов модел за моделирање мултикомпонентних и мултифазних струјања [72]. Овај модел је у литератури добро прихваћен и често примењиван, због једноставне имплементације и јако добрих резултата који се постижу његовом применом. У оквиру овог модела укључена је локална спољашња сила у једначинама за обе компоненте флуида, која узима у обзир међумолекуларне интеракције између компоненти. Овај модел има и одређена ограничења, као што је потенцијална нумеричка нестабилност симулација у случају да постоји веома велика разлика у густини између два флуида, али су и ова ограничења превазиђена у неким побољшаним верзијама овог модела [73]. Ипак и поред поменутих ограничења, овај Докторска дисертација Поглавље 5 Тијана Ђукић 53 модел је успешно примењиван у симулацијама двофазног струјања крви [74], термалног струјања [75], струјања у микроканалима [76] и многим другим. Код симулација струјања вишекомпонентног флуида применом Shan-Chen-овог модела користе се исте једначине LB методе за симулацију сваке компоненте посебно.               eqieqiiiii Fftftftf         uuxxξx , 2 1 1,, 1 ,1, 0        (5.1) где је са горњим експонентом  означена компонента флуида. Нумеричким решавањем једначине (5.1) могу се израчунати густина и брзина сваке компоненте флуида:    q i if 1  ;    q i ii f 1 1     ξu (5.2) Да би се симулирала интеракција између две компоненте, најпре се дефинише потенцијал интеракције:        '',', xxxxxx  GV  (5.3) где је са  означена друга компонента флуида, а x' се такође односи на тачку простора у домену у коме се налази друга компонента флуида. У једначини (5.3) са G је означена Green-ова функција, која се може апроксимирати на следећи начин:       h h G ', ',0 xx xx   G (5.4) где је са h означено растојање између 2 тачке фиксне мреже (као што је већ дефинисано, пошто је у питању LB метода, 1h ). Величина означена са  у једначини (5.3) заправо је једнака густини флуида –    . Вредност G дефинише јачину интеракције између компоненти, а знак овог коефицијента указује на то да ли се ради о привлачној или одбојној врсти интеракције. Докторска дисертација Поглавље 5 Тијана Ђукић 54 Пошто је дефинисан потенцијал интеракције, сада се може одредити сила интеракције између два флуида:             q i ii 1 ξξxxF G (5.5) Још једна корекција коју је потребно увести је везана за израчунавање равнотежне функције расподеле, за обе компоненте флуида. Наиме, равнотежна функција расподеле се израчунава применом тзв. композитне брзине:              1 1 ' 1 q i ii fξ u (5.6) Као што је раније већ наглашено (у једначинама (2.106) и (4.7)), потребно је укључити и ефекат спољашње силе приликом израчунавања макроскопске брзине, тако да се равнотежна брзина (која се користи за израчунавање равнотежне функције расподеле) за сваку компоненту израчунава на следећи начин:        F uu  'eq (5.7) Ове измене не утичу на основни део солвера за струјање флуида, а омогућавају јако ефикасну симулацију струјања вишекомпонентног флуида. Такође, све измене је једноставно имплементирати у LB солвер, а притом је и даље могуће паралелизовати цео софтвер, с обзиром да је задржана и та карактеристика LB методе да су сва израчунавања везана за појединачне чворове, односно да су међусобно независна и неспрегнута, па се може применити техника паралелизације која је описана у поглављу 6. 5.2 Моделирање струјања више флуида различите вискозности раздвојених мембраном солида Тема овог рада је моделирање понашања црвених крвних зрнаца у малим крвним судовима. Као што је већ објашњено, структура црвених крвних зрнаца је таква да она имају танку мембрану и унутрашњост је испуњена нестишљивом течношћу, цитоплазмом [3]. Густина цитоплазме унутар ћелија и крвне плазме у крвном суду је слична, тако да се за мале Reynolds-ове бројеве та разлика може занемарити. Али, код Докторска дисертација Поглавље 5 Тијана Ђукић 55 нормалних крвних зрнаца, вискозност цитоплазме је око 5 пута већа од вискозности крвне плазме у коју су ћелије уроњене [4,5]. Због тога би у нумеричким симулацијама требало узети у обзир и овај фактор. Међутим, за овакве симулације није применљив претходно описани Shan-Chen-ов модел, пошто овде не долази до интеракције између два флуида. С обзиром да је једина разлика у вискозности флуида, поново је могуће користити исте једначине, а чак није потребно ни да се решавају два независна флуида, већ је довољно да се у одређеним тачкама домена честицама додели друга вискозност. Наиме, вискозност је у LB методи повезана са временом релаксације  (што је дефинисано релацијом (2.103)), тако да је довољно да одређени чворови мреже имају другу вредност времена релаксације и то аутоматски значи да ће у тим чворовима флуид имати другу вредност вискозности. Ово је шематски приказано на слици 5.1. Током кретања ћелије потребно је кориговати вредности вискозности у чворовима који се налазе унутар и ван ћелијске мембране. Слика 5.1 – Шематски приказ два флуида различите вискозности раздвојених мембраном солида Постоји више начина на који се ово ажурирање вредности вискозности у чворовима може извршити. Може се направити индексно поље, да би се одредила релативна позиција сваког чвора домена флуида у односу на ћелијску мембрану. Tryggvasson и други [77] су у свом раду решавали Poisson-ову једначину у свакој итерацији да би одредили поменуто индексно поље. Francois и други [78] су предложили да се карактеристике флуида ажурирају у складу са нормалним растојањем од мембране. У овом раду је примењен овај други приступ, пошто је компјутерски ефикаснији. Са r се дефинише вектор који повезује чвор из домена флуида (из Euler-ове мреже) са њему најближим чвором мембране солида (из Lagrange-ове мреже). Оно што је карактеристично јесте да се уз интензитет овог вектора додаје и знак, који указује на то да ли је чвор флуида унутар мембране солида (када се узима да је r позитивно) или је чвор флуида ван мембране солида (када се узима да је r негативно). Докторска дисертација Поглавље 5 Тијана Ђукић 56 Дефинише се и Heaviside-ова функција:                     h hh h h r h r H d i ii 2 22 2 ,1 , 2 sin 1 2 1 2 1 ,0 1 r r r r   (5.8) Променљива карактеристика флуида, у овом случају вискозност, се може повезати са растојањем флуидног чвора од мембране солида, применом овако дефинисане Heaviside-ове функције:       xrx Hinoutout   (5.9) где се доњи индекси in и out односе на вредности вискозности унутар и ван мембране солида, респективно. Из једначине (5.8) се јасно види да се током кретања солида могу мењати карактеристике само оних чворова флуида чије је растојање мање од 2, тако да је то додатно олакшање и убрзање при имплементацији, пошто није потребно ажурирати карактеристике свих чворова из домена флуида, већ само оних који се налазе врло близу солида. Тијана Ђукић 57 6 Паралелизација развијеног софтвера У претходним поглављима описане су теоријске основе метода које су коришћене за моделирање кретања деформабилних тела у флуиду. На основу математичког модела развијен је програм који је коришћен за извршавање симулација које су приказане у поглављу 7, а симулације обухватају тродимензионално струјање флуида у комплексном домену, као и кретање деформабилних честица кроз овако комплексан домен. Овај софтвер је веома захтеван са аспекта компјутерских ресурса и једна симулација може да траје неколико сати. Због тога је извршена паралелизација софтвера за рад на GPU уређајима. GPU (енг. Graphics Processing Unit) уређаји су савремена верзија графичких картица, имају велики број процесора и могу паралелно да обрађују податке. То практично значи да се савремене графичке картице такође могу користити за извршавање апликација које су прилагођене за извршавање на GPU уређајима. GPU уређаји су посебно намењени решавању проблема који се другачије могу описати као извршавање великог броја аритметичких операција над неким (такође великим) скупом података. Укратко, када треба исти низ операција извршити на великом броју елемената, паралелно, са малим бројем меморијских операција и великим бројем аритметичких операција – пожељно је користити GPU уређаје. На стандардном рачунару, овакав програм би се извршавао секвенцијално, док се на GPU уређају подаци преносе у меморију уређаја и деле се на хиљаде нити (енг. threads) које се паралелно процесирају. Апликације које користе низове или матрице великих димензија су посебно погодне за паралелизацију применом овог принципа и тада се израчунавања драстично убрзавају. У овом раду су принципи GPU програмирања примењени на симулацију струјања флуида и кретања деформабилних тела. Поступак паралелизације и добијено убрзање ће бити тема овог поглавља. 6.1 CUDA архитектура за програмирање GPU уређаја NVIDIA је један од највећих светских произвођача графичких картица. Током свог развоја, графичке картице су постајале све програмабилније, што је довело до тога да NVIDIA представи свој први GPU уређај, који је унапређена верзија њихових графичких картица и који може да се користи за опште намене, а не само за приказ графике на екрану. То је означило почетак нове технологије GPGPU (енг. General Purpose GPU). NVIDIA је у новембру 2006. представила и CUDA архитектуру (енг. Compute Unified Докторска дисертација Поглавље 6 Тијана Ђукић 58 Device Architecture) [43], која заправо представља софтвер потребан за програмирање хардвера – GPU уређаја. На овај начин је програмирање за GPU постало далеко лакше и то уз помоћ програмских језика вишег нивоа, као што су C и C++. CUDA архитектура проширује програмски језик C, са скупом нових резервисаних речи помоћу којих се GPU уређају задају команде. Велика предност ове архитектуре је и то што се програмирање може вршити у Microsoft Visual Studio окружењу, као да су у питању стандардне апликације. CUDA омогућава покретање кода на две различите платформе, на тзв. host систему (који представља CPU процесор) и на тзв. device систему (који се састоји од једног или више GPU уређаја). Постојећи програм се прилагођава за рад са GPU уређајима тако што се одваја део програма у коме се извршава велики број аритметичких операција на великом скупу података и тај део програма се извршава на GPU. На тај начин GPU ради као копроцесор главног процесора. За то време CPU организује извршавање програма и координира разменом потребних информација. Слика 6.1 - Илустрација мреже блокова и нити [43] Докторска дисертација Поглавље 6 Тијана Ђукић 59 Да би се део програма извршавао на GPU уређају, потребно је испрограмирати специјалне функције, које се називају кернели (енг. kernels). Свака оваква функција се позива из главног дела програма од стране CPU-а и подпрограм написан у оквиру тела кернел функције извршава се паралелно на мрежи нити (енг. grid of threads). Приликом сваког позива кернел функције, потребно је дефинисати број нити и димензије мреже (што дефинише начин на који су те нити распоређене) и наравно проследити одговарајуће аргументе. Број нити и њихов распоред је битан да би GPU уређај могао да одреди колико ресурса треба да ангажује. На слици 6.1 илустрован је изглед мреже блокова и нити у оквиру једног блока, за један позив кернел функције. Променљивом dimGrid дефинише се број блокова у мрежи, а променљивом dimBlock дефинише се број нити у оквиру сваког блока. Наравно, постоје извесна ограничења када су у питању димензије мреже блокова и сваког блока и то зависи од GPU уређаја на коме се кернел функција извршава. Када је у питању Tesla C1060 која је коришћена за паралелизацију LB солвера, блок може имати највише 512 нити, а максимална димензија мреже блокова је 65.535 нити. Више детаља о карактеристикама GPU уређаја дато је Додатку Ц. Да би се у оквиру кернел функције приступало елементима низа или матрице, потребно је знати индекс сваке нити. Свака нит приликом покретања кернел функције добија свој дефинисани идентификациони број (ID) и на основу њега се одређује тачно место те нити у низу или матрици. Сваки блок има свој ID број којим је одређена позиција тог блока у мрежи блокова и ту информацију садржи системски дефинисана променљива blockIdx. Системски дефинисана променљива threadIdx памти позицију нити у оквиру одговарајућег блока. Помоћу ове две променљиве одређује се глобални индекс нити и помоћу њега се приступа одговарајућем елементу матрице (или низа). 6.2 Паралелизација софтвера на бази lattice Boltzmann методе У оквиру нумеричке имплементације lattice Boltzmann методе постоје два основна корака – корак судара и корак пропагације, дефинисани једначинама (2.101) и (2.102). Ова два корака се наизменично понављају у одређеном броју итерација. Сваки од ова два корака се мора применити на све чворове мреже, пре него што почне следећи корак. Корак судара је потпуно локалног карактера, пошто једначина овог корака указује да за израчунавање функције расподеле у посматраном чвору нису потребни подаци о било ком другом чвору. У оквиру корака пропагације потребни су подаци о посматраном чвору и његовим најближим суседним чворовима, али се опет операције могу извршавати за сваки чвор посебно. Проблем симулације применом LB методе се може свести на проблем у коме се велики скуп података може поделити на мање подскупове и над сваким подскупом се може симултано извршавати одређени број углавном аритметичких операција. Због тога је LB метода изузетно погодна за паралелизацију и извршавање на GPU уређајима. На слици 6.2 приказан је алгоритам рада LB солвера који се извршава секвенцијално, са означеним делом који треба паралелизовати. Такође су наведени делови кода помоћу Докторска дисертација Поглавље 6 Тијана Ђукић 60 којих се имплементирају корак судара и корак пропагације. Потребно је испрограмирати две кернел функције које ће се позивати уместо ових еквивалентних секвенцијалних функција. Постоји једна важна ставка о којој треба водити рачуна када се паралелизује постојећи програм, а то је манипулација меморијом. Наиме, CUDA архитектура је тако конципирана да CPU (другачије назван host, у преводу „домаћин“) и GPU (другачије назван device, у преводу „уређај“) имају одвојене DRAM меморије (што у суштини и јесте случај, ако се посматра са хардверске стране) и ове меморије се називају host и device меморија, респективно. CUDA омогућава програму да приступа овим меморијама независно, али при томе постоје одређена ограничења. Део програма који се извршава на главном процесору (CPU) може користити само host меморију, док кернел функције могу приступати само device меморији. Због тога се мора имплементирати експлицитно алоцирање и деалоцирање меморије за оба система, као и трансфер података између ова два типа меморије. Постоје посебне CUDA функције које обезбеђују пренос података са CPU на GPU и обрнуто. Слика 6.2 - Алгоритам рада LB солвера са означеним деловима које треба паралелизовати Докторска дисертација Поглавље 6 Тијана Ђукић 61 Најједноставнији начин паралелизације у случају LB солвера би био да се сви потребни подаци о чворовима (подаци о функцији расподеле, спољашњим силама, граничним условима ...) чувају у host меморији, а да се пре сваког позива кернел функције изврши пренос података у device меморију, па да се ти подаци поново врате назад после завршетка рада кернел функције. Међутим, овај приступ је потпуно погрешан, због тога што је трансфер података између ова два типа меморије изузетно „скупа“ операција, која „троши“ много рачунарског времена. Потребно је избегавати трансфере између ова два типа меморије кад год је то могуће. У овом конкретном случају, то значи да треба применити другачији приступ. Сви потребни подаци треба најпре да се иницијализују на CPU, потом се само једном преносе у device меморију и током извршавања програма се више не преносе. На тај начин је кернел функцијама обезбеђен приступ подацима у device меморији и ови подаци се могу мењати по потреби. За приказивање резултата (поља макроскопске брзине и притиска) назад у host меморију прослеђују се само поменуте макроскопске величине, које се рачунају помоћу посебне кернел функције. Треба имати у виду да се најчешће поље брзине и притиска не приказују за сваку итерацију LB солвера, већ на 100, 1000 или више итерација, у зависности од симулације. На овај начин је драстично смањена количина података која се преноси између меморија, што додатно убрзава рад програма. У оквиру CUDA архитектуре, device меморија је организована у једну хијерархијску структуру, тако да нити током извршавања кернел функција могу да приступају подацима из меморије на више нивоа. Свака нит има своју приватну локалну меморију, којој само та нит има приступ. Сваки блок нити има своју дељиву меморију (енг. shared memory) која је доступна свим нитима у оквиру тог блока. И на крају, све нити имају приступ глобалној меморији уређаја. Наравно, различито време је потребно да би једна нит приступила податку из локалне (регистарске), дељиве или глобалне меморије. Најбржи је приступ локалној, затим дељивој, а најспорији је приступ глобалној меморији [43,79] (у неким случајевима је и до 150 пута спорији од регистарске меморије). Зато треба избегавати приступ глобалној меморији кад год је то могуће. Ако се ове чињенице узму у обзир приликом паралелизације LB солвера, долази се до закључка да је погодније да се уместо две независне кернел функције за корак судара и корак пропагације, имплементира једна функција која спаја ова два корака. На овај начин се ништа не мења у самој теоријској поставци, пошто се та два корака и даље извршавају одвојено, већ је то само питање имплементације. Оваквим приступом се избегава да се вредности функције расподеле између ова два корака записују у глобалну device меморију и потом поново учитавају, већ се памте у регистарској меморији сваке нити појединачно. Без обзира на начин имплементације, корак моделирања судара се мора завршити за све чворове пре него што почне корак пропагације. Да би се то обезбедило користи се системска CUDA функција за синхронизацију свих нити __syncthreads(). Ова функција обезбеђује да све нити заврше са извршавањем свих претходних операција, пре наставка извршавања кернел функције. Докторска дисертација Поглавље 6 Тијана Ђукић 62 Слика 6.3 – Алгоритам извршавања паралелизованог LB солвера – моделирање кретања деформабилног тела у флуиду Докторска дисертација Поглавље 6 Тијана Ђукић 63 Поред паралелизације корака судара и пропагације, било је неопходно паралелизовати и функције које су везане за солид-флуид интеракцију. Уколико би те функције остале непромењене, то би поново узроковало трансфере података између host и device меморије, што није добро решење. Поступак моделирања солид-флуид интеракције у паралелизованој верзији је еквивалентан секвенцијалном алгоритму који је приказан на слици 4.4, само је потребно применити принципе GPU програмирања и написати три кернел функције. Прва кернел функција врши интерполацију брзина граничних тачака солида и ажурира позицију солида и овде свака нит обрађује податке о једном чвору граничне површине солида. Друга кернел функција одређује реакцију солида на деформацију и овде свака нит обрађује податке о једном елементу солида. Трећа врши интерполацију спољашње силе којом солид делује на флуид и овде свака нит обрађује податке о једном чвору флуида. Такође, извршена је и паралелизација дела секвенцијалног програма који је везан за симулацију вишекомпонентног флуида, програмирањем додатних кернел функција. У случају моделирања применом Shan-Chen-овог модела, направљене су измене у кернел функцији која симулира кораке судара и пропагације, док је у случају моделирања флуида различите вискозности имплементирана додатна кернел функција која ажурира карактеристике флуида у сваком чвору мреже после сваке итерације. На слици 6.3 приказан је алгоритам извршавања паралелизованог LB солвера, за случај моделирања кретања деформабилног тела у флуиду. 6.3 Паралелизација на кластерима са више GPU уређаја До сада је објашњен поступак паралелизације софтвера за рад на једном GPU уређају. Наравно, уколико би се уместо једног користило више GPU уређаја, могуће је постићи још веће убрзање паралелизованих програма. Постоје две могућности покретања програма на више GPU уређаја. Први случај је покретање на кластеру код кога су сви GPU уређаји повезани на исти CPU процесор и код овакве архитектуре се користи OpenMP приступ [80]. Други случај је покретање на кластеру код кога су GPU уређаји повезани са више CPU процесора, тако да постоји више тзв. node-ова у архитектури и у оваквим случајевима се користи MPI приступ [81]. Без обзира на архитектуру кластера, процес паралелизације је врло сличан. Једина разлика је у томе што се комуникација између GPU уређаја различито остварује за различите архитектуре. Што се тиче програмирања кернел функција, у том делу такође има мало измена. Користе се исте кернел функције које су испрограмиране за покретање програма на једном GPU уређају. Наравно, пошто се користе ресурси више GPU уређаја у исто време, кернел функције се морају позивати независно за сваки GPU уређај. То аутоматски значи да потребни подаци морају бити алоцирани и иницијализовани у device меморији сваког GPU-а посебно. Зато је неопходно најпре извршити декомпозицију, односно поделу целог домена који се симулира на N делова, где је N број GPU уређаја који се користе у току симулације. У овом раду декомпозиција је Докторска дисертација Поглавље 6 Тијана Ђукић 64 извршена у правцу x осе, као што је илустровано на слици 6.4 на примеру дводимензионалног домена. На овај начин се постиже да сваки GPU буде „задужен“ за само један део домена и да обрађује податке везане за тај део домена. Као што је на почетку програма потребно све податке распоредити на N GPU уређаја, тако је и приликом исписивања резултата (поља макроскопске брзине и притиска) потребно сакупити све податке са свих GPU уређаја назад у host меморију, да би могли да буду приказани или одштампани у одговарајућем формату. Слика 6.4 – Декомпозиција домена симулације дуж x осе Већ је дискутовано да је корак симулације судара локалног карактера, па због тога у току овог корака није потребна никаква синхронизација података узмеђу GPU уређаја. Међутим, у кораку пропагације, користе се и подаци из суседних чворова мреже, па је због тога неопходно у свакој итерацији извршавати пренос података са сваког GPU уређаја на његова два најближа суседа – на први суседни GPU уређај са леве и десне стране од посматраног GPU уређаја, односно на GPU уређај са ID бројем за један мањим и за један већим. У свакој итерацији се дакле најпре извршавају потребна израчунавања, а затим се врши размена података. Подаци које је потребно разменити су одређене компоненте функције расподеле. На слици 6.5 су приказане компоненте које се прослеђују у случају симулације дводимензионалног струјања флуида, док су на слици 6.6 приказане компоненте које се прослеђују у случају симулације тродимензионалног струјања флуида. Докторска дисертација Поглавље 6 Тијана Ђукић 65 Слика 6.5 – Компоненте функције расподеле које се размењују између два суседна GPU уређаја – случај симулације дводимензионалног струјања флуида Слика 6.6 – Компоненте функције расподеле које се размењују између два суседна GPU уређаја – случај симулације тродимензионалног струјања флуида Кад је у питању размена података, односно комуникација између GPU уређаја, битан је тип архитектуре. Најпре ће бити објашњено како се остварује комуникација у случају када су сви GPU уређаји повезани са једним CPU процесором. У овом случају програм се позива само са главног CPU процесора, који контролише све GPU уређаје и покреће кернел функције на свим уређајима. Код овакве архитектуре, два GPU уређаја могу директно да комуницирају и није потребно да CPU посредује у комуникацији. Размена података између GPU уређаја када се примењује OpenMP приступ илустрована је на слици 6.7. Докторска дисертација Поглавље 6 Тијана Ђукић 66 Слика 6.7 – Размена података између GPU уређаја када се примењује OpenMP приступ У случају да су GPU уређаји повезани са више CPU процесора, програм се покреће на више CPU процесора и користи се N MPI нити (енг. threads) и свака MPI нит контролише један GPU уређај. У овом случају, два GPU уређаја не могу директно да комуницирају међусобно, пошто их не контролише исти CPU процесор. Због тога се подаци са једног GPU уређаја (означеног са GPU_1) на други (означен са GPU_2) преносе у три корака: 1. пренос података из device меморије GPU_1 у host меморију CPU_1 2. пренос података са CPU_1 на CPU_2 коришћењем MPI функције MPI_Sendrecv 3. пренос података из host меморије CPU_2 у device меморију GPU_2 Овај процес је илустрован на слици 6.8. Слика 6.8 – Размена података између GPU уређаја када се примењује MPI приступ Докторска дисертација Поглавље 6 Тијана Ђукић 67 6.4 Поређење брзине рада секвенцијалне и паралелизоване верзије LB солвера У наставку ће бити приказано убрзање које се постиже када се користи паралелизована верзија LB солвера. Најпре је поређено време извршавања секвенцијалне верзије са једне стране и паралелизоване верзије која је покретана на једном GPU уређају са друге стране. Потом је поређено време извршавања паралелизоване верзије на више GPU уређаја применом OpenMP приступа. На крају је анализирано убрзање које се постиже на кластеру са великим бројем GPU уређаја, уз примену MPI приступа. 6.4.1 Покретање паралелизоване верзије на једном GPU уређају Најпре ће на једном примеру бити илустровано убрзање које се постиже у случају симулације струјања флуида, како у дводимензионалном, тако и у тродимензионалном домену на једном GPU уређају. Током тестирања извршавана је симулација стационарног струјања флуида у правој цеви. Број итерација је постављен на 20.000. Овај број итерација је одабран зато што се тада програм извршава у разумном временском интервалу, а довољан је да се достигне стационарно стање. Са друге стране, пошто су за сваку итерацију потребни исти рачунарски ресурси и исто време извршавања, резултати овог тест примера су релевантни да се донесе закључак о убрзању паралелизоване верзије и може се предвидети време извршавања комплекснијих симулација. Када се пореди брзина рада, потребно је узети у обзир утицај броја чворова, односно димензија матрица над којима се врше операције. Због тога је укупан број чворова LB мреже вариран и утврђено је да што је број чворова већи, убрзање постигнуто применом GPU уређаја је веће, што је било и очекивано. Одређивано је време извршавања секвенцијалне и паралелизоване верзије LB солвера како за симулацију дводимензионалног, тако и за симулацију тродимензионалног струјања флуида. Промена времена извршавања обе верзије у зависности од броја чворова LB мреже за дводимензионално струјање приказано је на дијаграму 6.9, док је иста промена за тродимензионално струјање приказана на дијаграму 6.10. Докторска дисертација Поглавље 6 Тијана Ђукић 68 Слика 6.9 – Упоредни приказ времена извршавања за различит број чворова LB мреже; случај дводимензионалног струјања Слика 6.10 - Упоредни приказ времена извршавања за различит број чворова LB мреже; случај тродимензионалног струјања На слици 6.11 приказано је како се убрзање паралелизованог LB солвера у односу на секвенцијални солвер мења у зависности од броја чворова LB мреже. Како се број чворова повећава, расте и убрзање. Вредности убрзања су израчунате на основу Докторска дисертација Поглавље 6 Тијана Ђукић 69 измереног времена извршавања само оног дела програма у коме се спроводе операције израчунавања (петља по итерацијама, без приказивања резултата). У том случају максимална вредност убрзања износи 21 и постиже се за 80.000 чворова. Даље повећање броја чворова нема утицаја на убрзање. Слика 6.11 - Дијаграм промене убрзања паралелизованог LB солвера у зависности од броја чворова LB мреже Тестирано је и убрзање паралелизованог програма, када се симулира кретање деформабилног тела у флуиду. Покретан је први пример који је представљен у поглављу 7.4.1 – деформација честице под дејством смичућег профила брзине. Број итерација је поново фиксиран, али овога пута на 2.000, пошто је то довољно да се постигне одређена деформација солида и да се предвиди убрзање програма. Вариран је и број чворова LB мреже којом је дискретизован домен флуида, али се као и претходним случајевима показало да убрзање достиже максималну вредност за 80.000 чворова и после се веома мало мења. Вариран је и број чворова мреже којом је дискретизован домен солида и ту је уочена промена убрзања, што се може видети на слици 6.12. За мањи број чворова, убрзање је слично, али са великим порастом броја чворова, долази до изражаја чињеница да се сви чворови солида у паралелизованој верзији обрађују синхронизовано. Даље повећање броја чворова солида није драстично повећало убрзање паралелизоване верзије. Докторска дисертација Поглавље 6 Тијана Ђукић 70 Слика 6.12 - Упоредни приказ убрзања програма за симулирање кретања деформабилног тела у флуиду за различит број чворова мреже солида Паралелизован је и програм који симулира струјање вишекомпонентног флуида. Број итерација је фиксиран на 2.000, пошто је у овом случају то довољан број итерација да се течности до одређене мере измешају и може да се стекне увид у рад програма. Број чворова је вариран, да би се утврдила зависност убрзања од броја чворова, која је приказана на слици 6.13. Слика 6.13 - Дијаграм промене убрзања паралелизованог LB солвера за симулирање вишекомпонентног струјања у зависности од броја чворова LB мреже Докторска дисертација Поглавље 6 Тијана Ђукић 71 6.4.2 Покретање паралелизоване верзије на више GPU уређаја, уз примену OpenMP приступа Претходно објашњена симулација тродимензионалног стационарног струјања у правој цеви је сада извршавана на конфигурацији која се састоји од 3 GPU уређаја Tesla C1060, који су повезани са истим CPU процесором. Карактеристике GPU уређаја Tesla C1060 дати су у Додатку Ц, у табели Ц.2. Коришћен је OpenMP приступ да би се програм прилагодио овој архитектури. Приликом тестирања поново је број итерација постављен на 20.000, а број чворова је вариран, да би се утврдила зависност убрзања од броја чворова. На дијаграму 6.14 приказано је време извршавања исте симулације на једном, два и три GPU уређаја. Коначно убрзање које је постигнуто када се пореди извршавање LB солвера на једном и два GPU уређаја је 1.9 пута, док се LB солвер на три GPU уређаја извршава 2.9 пута брже од солвера који се извршава на једном GPU уређају. Убрзање не расте линеарно са бројем GPU уређаја, што је последица чињенице да се одређено време потроши на комуникацију између GPU уређаја и пренос података из једне device меморије у другу. Промена убрзања приликом извршавања на два и три GPU уређаја у односу на извршавање на једном уређају је приказана на слици 6.15, при чему је приказана и зависност ових убрзања од броја чворова. Слика 6.14 - Упоредни приказ времена извршавања за различит број чворова LB мреже и различит број GPU уређаја на којима се апликација извршава Докторска дисертација Поглавље 6 Тијана Ђукић 72 Слика 6.15 - Дијаграм промене убрзања паралелизованог LB солвера на више GPU уређаја у зависности од броја чворова LB мреже 6.4.3 Покретање паралелизоване верзије на кластеру са великим бројем GPU уређаја, уз примену MPI приступа За тестирање рада програма на кластеру са великим бројем GPU уређаја који су повезани са више процесора, коришћен је CURIE суперкомпјутер. Он је у власништву GENCI институције, која се налази у Француској и део је PRACE удружења (енг. PRACE – Partnership for Advanced Computing in Europe). Током ангажовања на пројекту Центра за Биоинжењеринг, омогућен је приступ CURIE компјутерском кластеру и коришћење њихових ресурса у оквиру Preparatory Access позива за пројекте PRACE организације. На CURIE GPU уређајима покретане су комплексне симулације струјања флуида кроз артерије, али је и тестиран LB солвер оптимизован за рад са више GPU уређаја који су повезани са више CPU процесора. Детаљи конфигурације CURIE кластера дати су у Додатку Ц. Поново је покретана симулација тродимензионалног стационарног струјања флуида у правој цеви, као и у претходним случајевима тестирања рада LB солвера. На слици 6.16 приказано је убрзање које је постигнуто на већем броју GPU уређаја, у поређењу са временом извршавања паралелизоване верзије на једном GPU уређају. Приказани су тзв. strong scaling резултати. То значи да је исти скуп података обрађиван на променљивом броју GPU уређаја. Са дијаграма 6.16 се види да како се повећава број GPU уређаја, програм је мање ефикасан (одступање од линеарног убрзања је веће). Ово је последица чињенице да приликом извршавања на већем броју GPU уређаја због исте величине скупа података нису искоришћени сви расположиви ресурси и са повећањем броја уређаја, ресурси су све мање искоришћени. Због тога је убрзање мање од идеалног и време које се троши на комуникацију између Докторска дисертација Поглавље 6 Тијана Ђукић 73 GPU-ова има више утицаја на укупно време извршавања. На слици 6.17 такође је приказано убрзање на већем броју GPU уређаја, али овога пута приказани су тзв. weak scaling резулати. То значи да је сваком GPU уређају додељен исти скуп података за обрађивање, односно укупан скуп података је већи са већим бројем GPU-ова. Када се погледа дијаграм 6.17, јасно је да развијени програм користи на адекватан начин све ресурсе GPU уређаја и да време потрошено на комуникацију између GPU-ова нема претерано негативан утицај на укупне перформансе LB солвера. Слика 6.16 - Strong scaling убрзање LB солвера извршаваног на више GPU уређаја, у поређењу са извршавањем на једном GPU уређају Докторска дисертација Поглавље 6 Тијана Ђукић 74 Слика 6.17 - Weak scaling убрзање LB солвера извршаваног на више GPU уређаја, у поређењу са извршавањем на једном GPU уређају Тијана Ђукић 75 7 Резултати симулација Као део истраживачког рада који је претходио овом докторату, развијен је софтвер који представља нумеричку имплементацију LB методе, уз посебну имплементацију солид- флуид интеракције и кретања деформабилног тела у флуиду. У овом поглављу биће приказани резултати добијени применом овог софтвера. Резултати који су добијени су поређени или са аналитичким решењима, или са решењима која су добијена применом других метода у литератури. Сви примери се односе на тродимензионалне (просторне) проблеме, тако да је за прорачун коришћена мрежа са ознаком D3Q27, чија је структура приказана у Додатку Б. С обзиром да је LB метода дискретна и да посматра флуид као скуп честица, у оквиру ове методе се све физичке величине дефинишу и користе у свом бездимезионом облику, који је прилагођен теоријским основама саме методе. Због тога ће најпре бити објашњено како се у LB методи врши израчунавање физичких величина. LB метода се може успешно користити за симулације струјања крви кроз кардиоваскуларни систем човека, што ће бити показано добијеним резултатима симулација тродимензионалног струјања флуида у људским артеријама – у аорти, каротидној и коронарној артерији. Кроз ове примере ће уједно бити приказане и могућности развијеног софтвера да симулира струјања кроз комплексну геометрију, која се може задавати на различите начине - на основу параметарског модела, или директно на основу клиничких података добијених за конкретног пацијента. Потом ће бити приказани резултати симулације мешања две течности у једном суду. На крају ће бити приказане симулације кретања деформабилног тела у флуиду, како честица сферног облика, тако и честица облика црвеног крвног зрнца. Најпре ће бити приказано неколико примера на којима је тестирана тачност целокупног нумеричког модела, а потом и неколико примера симулација конкретних проблема у биоинжењерингу. 7.1 Конверзија бездимензионих величина у физичке величине У оквиру LB методе примењује се посебан принцип рачунања свих физичких величина, тако да се све макроскопске величине које су значајне за симулацију (и оне које се добијају као резултат симулације) дефинишу у тзв. систему јединица мреже (енг. system Докторска дисертација Поглавље 7 Тијана Ђукић 76 of lattice units). Због тога је пре покретања симулације, потребно одредити вредности свих параметара у бездимензионом облику. Такође, приликом интерпретирања резултата, потребно је трансформисати бездимензионе величине у реалне „физичке“ вредности. Постоје три макроскопске величине које се користе за изражавање свих осталих величина. То су време  st , дужина  ml и густина       3m kg  . Нека друга физичка величина x (на пример, брзина или притисак) се може изразити у функцији ове три величине:  ,,ltFx  (7.1) Да би се извршила конверзија бездимензионих величина у физичке величине, потребно је одредити фактор скалирања pF , тако да се може записати релација: pdp Fxx  (7.2) где су са px и dx означене вредности физичке величине у стандардном систему јединица и у бездимензионом облику, респрективно. Једна од карактеристичних величина, која се дефинише пре покретања LB симулације је карактеристична брзина (у бездимензионим јединицама мреже), која се означава са lbu . Овај параметар је јако важан, пошто се одговарајућим избором вредности карактеристичне брзине обезбеђује да струјање остане у области малог Mach-овог броја. На овој претпоставци се базира цела LB метода и то је узето у обзир током извођења једначина, па због тога овај услов мора да буде испоштован. То такође практично значи да се флуид задржава у области мале стишљивости. Вредност карактеристичне брзине мора да задовољи неједнакост: s lb c u 1  (7.3) Кинематска вискозност флуида и Reynolds-ов број повезани су релацијом:  uL Re (7.4) Најчешће се у симулацијама применом LB методе задаје нека карактеристична дужина домена, брзина и вискозност флуида у физичким јединицама. На основу ових података може се одредити Reynolds-ов број. Овај број је бездимензиона величина и он мора бити једнак у еквивалентним симулацијама, без обзира на јединице у којима се Докторска дисертација Поглавље 7 Тијана Ђукић 77 изражавају преостале величине. Поред поменутих параметара, у LB симулацијама се дефинише и резолуција, односно број чворова у LB мрежи у правцу сваке осе. На основу ових почетних података могуће је извршити конверзију физичких величина у бездимензионе и обрнуто. Другим речима, могу се одредити фактори скалирања за поменуте три макроскопске величине – дужину, време и густину. Фактор скалирања за дужину је једнак: lb p L L x  (7.5) где је pL физичка дужина, lbL је дефинисана резолуција lbN , односно број чворова којим се дискретизује посматрана дужина. Да би се испунио услов нестишљивости и малог Mach-овог броја, најчешће се бира вредност карактеристичне брзине 01,0lbu . На основу фактора скалирања за дужину и карактеристичне брзине, одређује се фактор скалирања за време: xut lb  (7.6) Фактор скалирања за густину се дефинише у зависности од јединица које се користе приликом дефинисања проблема. Уколико су макроскопске величине дефинисане у килограмима, метрима и секундама, онда је фактор скалирања за густину 1000 , а уколико су макроскопске величине дефинисане у грамима, центиметрима и секундама, онда је 1 . Пре покретања симулације, потребно је израчунати и време релаксације, као један од најбитнијих параметара симулације. Током извођења једначина LB методе, једначином (2.103) дефинисана је веза између времена релаксације и динамичке вискозности, па се може написати веза између кинематске вискозности и времена релаксације:        2 12  slb c (7.7) Кинематска вискозност у физичким величинама се може конвертовати у бездимензиони облик применом претходно дефинисаних фактора скалирања: 2x t lbp    (7.8) Докторска дисертација Поглавље 7 Тијана Ђукић 78 На основу једначина (7.7) и (7.8), ако је позната кинематска вискозност у физичким величинама и ако су претходно израчунати фактори скалирања, може се одредити време релаксације: 2 11 2 2     t x c p s  (7.9) Да би се обезбедила добра стабилност решења и LB симулације, потребно је да се води рачуна да добијена вредност времена релаксације задовољава неједнакост: 15,0  (7.10) и да ова вредност буде што ближа вредности 1. Уколико током дефинисања параметара не буде задовољен овај услов, треба кориговати неки други параметар, односно или карактеристичну брзину lbu или број чворова (резолуцију) N , да би се добило задовољавајуће време релаксације. 7.2 Симулације струјања крви кроз кардиоваскуларни систем Кардиоваскуларни систем је један од најважнијих система у људском организму. Због тога ће у наставку овог поглавља најпре бити дат кратак увод у кардиоваскуларни систем, са описом главних крвних судова, а затим ће бити приказани резултати симулација струјања крви кроз поменуте крвне судове. Циркулаторни систем се састоји од крвних судова чија је главна улога транспорт хранљивих материја и кисеоника до периферних ткива. Може се поделити на два дела: системску циркулацију која служи за транспорт хранљивих материја и плућну циркулацију у току које се врши оксидација хемоглобина у крви. Ова два дела су повезана преко срца које делује као пумпа крвног система. Системска циркулација се назива још и велики крвоток или периферна циркулација, зато што снабдева крвљу сва ткива изузев плућа. Циркулаторни систем се дели на артерије, артериоле, капиларе, вене и венуле. Артерије преносе крв из срца у ткива. Артериоле су последње мале гране артеријског система и имају улогу контролних валвула помоћу којих се омогућава или спречава пролаз крви у капиларе, у зависности од тренутне потребе ткива. Капилари имају улогу у размени течности, хранљивих материја, електролита, хормона и других супстанци. Венуле прикупљају крв из капилара и уливају се у веће вене. Вене транспортују крв из ткива у срце. Докторска дисертација Поглавље 7 Тијана Ђукић 79 Срце пумпа крв континуално у аорту при средњем притиску од 100 mmHg (12.3 kPa). Ово пумпање је пулзаторно и због тога крвни притисак у артеријама осцилује између 120 mmHg у систоли (периоду контракције срца) до 80 mmHg у дијастоли (периоду релаксације срца). Са смањењем величине крвног суда смањује се и крвни притисак, па на крају циркулаторног система у капиларима средњи притисак износи око 17 mmHg. Просечни проток крви у циркулацији одраслог човека износи 5 литара у минуту. Како кроз сваки крвни суд пролази иста количина крви у јединици времена, брзина крви се мења у зависности од пречника крвног суда. У стању мировања брзина кроз аорту износи око 30 cm/s, а у капиларима је брзина хиљаду пута мања и износи око 0.30 mm/s. У хемодинамици се најчешће проучава струјање крви кроз три артерије, а то су: абдоминална аорта, каротидна бифуркација и лева коронарна артерија. Ови крвни судови су најзначајнији, са једне стране зато што најчешће подлежу болестима и са друге стране зато што најбоље репрезентују струјање у целом артеријском стаблу. У наставку овог поглавља биће приказани резултати симулација струјања крви кроз поменуте крвне судове. Уједно ће кроз приказ резултата ових симулација бити представљене и могућности развијеног LB солвера када је у питању задавање геометрије различитим пре-процесирањем података. 7.2.1 Моделирање струјања крви кроз људску аорту Аорта је велики крвни суд који полази из срца, тачније из леве срчане коморе и пролази уз задњи зид абдомена у правцу средње линије. Аорта носи и дистрибуира крв богату кисеоником до свих артерија системског крвотока. Из проширења почетног дела аорте (који се назива bulbus) полазе две коронарне артерије, а из лука се гранају артерије које снабдевају крвљу главу, врат и горње екстремитете. Аорта касније прелази у трбушну аорту, која се рачва у две илијачне артерије које снабдевају крвљу доње екстремитете. На слици 7.1 је приказан положај аорте и поменутих срчаних артерија у људском организму. Људска аорта има компликовану геометрију која укључује закривљења која нису у једној равни. Због тога је тешко направити компјутерски генерисан модел аорте. Модел аорте потребан за симулацију направљен је директно на основу клиничких података о једном конкретном пацијенту. Наиме, на основу DICOM снимака добијених са 64- MSCT (multi-slice Computed Tomography) скенера извршена је тродимензионална реконструкција. За анализу DICOM снимака коришћен је посебно развијен софтвер за пре-процесирање података потребних за LB солвер. Овај софтвер одређује пикселе снимка који припадају аорти. На тај начин се креира тзв. вокселизовани модел, који се користи при дефинисању геометрије у оквиру LB солвера. Сам модел се ограничава и прави се „оквир“ око аорте, тако да цео домен који се разматра буде у облику квадра и они чворови који одговарају обележеним пикселима са снимка су чворови у којима постоји флуид, док се остали чворови овако ограниченог домена третирају као зидови. Докторска дисертација Поглавље 7 Тијана Ђукић 80 Слика 7.1 – Положај аорте и срчаних артерија у људском организму На слици 7.2 приказан је изглед поменутог софтвера за анализу DICOM снимака. Са леве стране је приказан један пресек (енг. slice), а са десне је приказана аорта која је реконструисана са свих снимака. Црвени правоугаоник означава где се тренутно посматрани пресек налази у оквиру креираног модела аорте. Слика 7.2 – Софтвер за генерисање геометрије на основу DICOM снимака Што се тиче физичких параметара потребних за симулацију, сматра се да је крв нестишљив Њутновски флуид, који има густину 305,1 cmg и кинематску вискозност scm 2035,0 . Гранични услови су задати тако да је дефинисан параболични профил брзине на улазу и на излазима је дефинисано ограничење да нема Докторска дисертација Поглавље 7 Тијана Ђукић 81 површинског напона на излазима. Такође, на зидовима аорте брзина је једнака нули. Карактеристична брзина флуида је постављена на 01,0lbu . Укупан број чворова у LB мрежи износи 3.781.701. На слици 7.3 приказана је расподела притиска добијена симулацијом, док је на слици 7.4 приказана расподела брзина. Како је брзина на границама домена једнака нули, на слици су приказане струјне линије, уз одговарајућом бојом означен интензитет брзине дуж струјних линија. Слика 7.3 – Расподела притиска у људској аорти Слика 7.4 – Струјне линије и расподела брзина у људској аорти Докторска дисертација Поглавље 7 Тијана Ђукић 82 7.2.2 Моделирање струјања крви кроз каротидну бифуркацију Каротидне артеријске бифуркације се налазе са обе стране у пределу врата и снабдевају мозак и лице крвљу. Бифуркацију чини заједничка каротидна артерија која се рачва на унутрашњу каротидну артерију која снабдева мозак и спољашњу која снабдева лице крвљу. Артерије се после бифуркације рачвају под углом од око 25° у односу на осу улазне артерије. На слици 7.5 је шематски приказана каротидна артерија и њен положај у људском организму. Каротидна артерија је од посебног значаја јер се код пацијената у највећем проценту овде јавља стеноза. Највећи пречник се налази одмах иза рачвања на унутрашњој каротидној артерији и тај део се назива синус. Управо у том делу, на спољашњем зиду, јавља се атеросклеротични плак. Слика 7.5 – Положај каротидне артерије у људском организму Каротидна артерија се може параметарски дефинисати. Пречник заједничке каротиде износи cmD 61.0 и ово је основни параметар, на основу кога се одређују сви остали параметри којима се описује каротидна артерија. На слици 7.6 приказани су сви остали геометријски подаци [82,83]. Овакав модел каротиде се незнатно разликује од реалне анатомије код човека, али је у складу са литературом, симулација струјања спроведена на оваквом моделу [83,84]. Докторска дисертација Поглавље 7 Тијана Ђукић 83 Слика 7.6 – Геометријски подаци за модел каротидне артерије (равански пресек) Што се тиче физичких параметара потребних за симулацију, поново се (као и у претходном примеру) сматра се да је крв нестишљив Њутновски флуид, који има густину 305,1 cmg и кинематску вискозност scm2035,0 . Гранични услови за овај проблем су задати тако да је брзина на зидовима крвног суда једнака нули, на улазу је задат параболични профил брзине, а на излазу обе артерије задат је услов да су површинске силе једнаке нули. Ова симулација је покретана коришћењем улазног фајла који се иначе користи у софтверском пакету PakF, који представља имплементацију методе коначних елемената (МКЕ) [2,35]. Наиме, током креирања софтвера на бази LB методе, омогућено је да се за дефинисање параметара симулације користи и улазни фајл који је креиран за PakF. Оно што је потребно додатно урадити, а што је ван стандардне и до сада описане процедуре када је у питању покретање LB симулације, јесте вокселизација дефинисаног домена, пошто се LB симулације спроводе на мрежи фиксираних чворова, а не на мрежи произвољно распоређених елемената у простору, као што је случај са PakF-ом. Вокселизација је спроведена применом процедуре описане у литератури [85] при чему је са LB солвером повезан додатни open-source софтвер Binvox који је развио Patrick Min [86]. Вокселизовани модел се даље користи да би се дефинисали чворови LB мреже. Наиме, они чворови који су приликом вокселизације означени са 1 (значи да су унутар модела), заправо су чворови у којима постоји флуид. Са друге стране, чворови који су приликом вокселизације означени са 0 (значи да су ван модела) су означени као Bounce-back чворови, у овим чворовима се не решавају класичне једначине LB методе и они представљају зидове од којих се флуид одбија. За сваки коначни елемент у оквиру геометрије модела каротидне артерије може се одредити вектор нормале тог елемента. Овај вектор нормале се у LB симулацији користи да би се на основу једначине (2.115) одредиле вредности смичућег напона у свакој тачки зида. Резултати прорачуна LB методом се штампају коришћењем мреже коначних елемената, односно израчунате вредности се интерполирају и приказују у чворним тачкама мреже коначних елемената. На тај начин резултати се приказују на Докторска дисертација Поглавље 7 Тијана Ђукић 84 глатком моделу, а не на грубљем вокселизованом моделу који је коришћен за прорачун и постиже се прегледнија визуелизација резултата. Резултати добијени применом LB методе и методе коначних елемената приказани су на сликама 7.7, 7.8 и 7.9. На слици 7.7 приказана је расподела притиска, на слици 7.8 приказане су струјне линије, при чему је одговарајућом бојом означен интензитет брзине дуж струјних линија, а на слици 7.9 је приказана расподела смичућих напона. Као што се може видети, добијено је добро поклапање решења, а грешке између израчунатих вредности у оба софтвера су мале. Разлика у расподели смичућег напона у једној и другој симулацији је последица чињенице да је симулација у PakF-у спровођена за случај да су зидови каротидне артерије деформабилни, што није случај са LB методом. Слика 7.7 – Расподела притиска у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) Слика 7.8 – Струјне линије и расподела брзина у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) Докторска дисертација Поглавље 7 Тијана Ђукић 85 Слика 7.9 – Расподела смичућих напона у каротидној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) 7.2.3 Моделирање струјања крви кроз коронарну артерију Коронарне артерије се још називају и срчане артерије и то су крвни судови који снабдевају крвљу срчани мишић. Постоје две коронарне артерије – лева и десна и оне полазе из почетног дела аорте. На слици 7.10 су приказане коронарне артерије и њихов положај у људском организму. Слика 7.10 – Положај коронарних артерија у људском организму Циљ ове симулације је да се утврди како се струјање флуида одиграва у левој коронарној артерији (енг. left coronary artery – LCA). Комплексна тродимензионална геометрија модела потребна за симулацију креирана је за једног конкретног пацијента. Подаци о овом пацијенту су прикупљени у оквиру FP7 ICT пројекта ARTreat, у Институту за клиничку физиологију Националног истраживачког центра у Пизи у Италији (Institute of Clinical Physiology, National Research Council, Pisa, Italy). Докторска дисертација Поглавље 7 Тијана Ђукић 86 Симулације струјања крви су у оквиру овог пројекта вршене применом методе коначних елемената. Циљ симулација је био да се утврде критична места на којима је дошло или може да дође до настанка атеросклеротичног плака. У овом раду ће симулација струјања бити спроведена применом LB методе и резултати за једног пацијента ће бити упоређени са резултатима добијеним применом МКЕ. Наравно, како циљ овог рада нису првенствено симулације струјања у великим крвним судовима, овај пример само треба да илуструје могућност развијеног софтвера да симулира и струјања крви у већим крвним судовима комплексне геометрије. На основу DICOM снимака добијених са 64- MSCT (multi-slice Computed Tomography) скенера извршена је тродимензионална реконструкција геометрије употребом комерцијалног софтвера Mimics. Да би се добио бољи квалитет геометрије генерисаног модела, по потреби је вршена мануелна контрола и корекција геометрије. За потребе симулације применом LB методе поново је вршена вокселизација геометрије, као што је описано за случај каротидне артерије у поглављу 7.2.2. Једина разлика је што је у овом случају уместо улазног фајла за PakF геометрија задата коришћењем STL фајла, који је учитаван директно у LB солвер, па је коришћењем мреже троуглова вршена вокселизација. У оквиру STL фајла је за сваки троугао већ дефинисан вектор нормале, који је у LB солверу коришћен за израчунавање смичућег напона на зидовима. Физички параметри и гранични услови су исти као и у претходним симулацијама - сматра се да је крв нестишљив Њутновски флуид, који има густину 305,1 cmg и кинематску вискозност scm 2035,0 . Гранични услови су задати тако да је брзина на зидовима крвног суда једнака нули, на улазу је задат параболични профил брзине, а на свим излазима је задат услов да су површинске силе једнаке нули. Слика 7.11 - Расподела притиска у коронарној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) Докторска дисертација Поглавље 7 Тијана Ђукић 87 Резултати добијени применом LB методе могу се упоредити са резултатима добијеним применом методе коначних елемената. На слици 7.11 приказана је расподела притиска, а на слици 7.12 приказане су струјне линије, при чему је одговарајућом бојом означен интензитет брзине дуж струјних линија. На овој слици се уједно може видети како изгледа вокселизована геометрија. На слици 7.13 приказана је расподела смичућих напона. Као што се може видети, добијено је добро поклапање решења. Слика 7.12 - Струјне линије и расподела брзина у коронарној артерији – симулација применом LB методе Слика 7.13 - Расподела смичућих напона у коронарној артерији – симулација применом LB методе (лево) и применом МКЕ (десно) Докторска дисертација Поглавље 7 Тијана Ђукић 88 7.3 Симулација мешања две течности Овде ће бити приказани резултати симулације мешања две течности. Течности које су разматране су вода (густине 3 1 1000 mkg и кинематске вискозности sm271 1094.8  ) и уље (густине 32 920 mkg и кинематске вискозности sm251 10804.8  ). Оба флуида се налазе у округлом суду и на почетку симулације горњу половину посуде заузима вода, док је доле уље. На оба флуида делује гравитациона сила. Мешање течности је симулирано применом Shan-Chen-овог модела, који је већ описан у поглављу 5.1. На слици 7.14 приказан је симулирани процес мешања после одређеног броја итерација, при чему је приликом визуелизације урађено неколико пресека суда, да би се прегледније пратио ток симулације. почетни тренутак – после 1 итерације после 1000 итерација после 2000 итерација после 3000 итерација Докторска дисертација Поглавље 7 Тијана Ђукић 89 после 5000 итерација после 7000 итерација после 8000 итерација после 10000 итерација после 12000 итерација после 15000 итерација Слика 7.14 – Моделирање процеса мешања две течности, у више временских тренутака Докторска дисертација Поглавље 7 Тијана Ђукић 90 7.4 Симулације струјања крви кроз капиларне крвне судове Крв представља двокомпонентни систем који се састоји од ћелија (формираних ћелијских елемената) и крвне плазме. Постоји неколико типова крвних ћелија, а то су еритроцити, леукоцити и тромбоцити. Црвена крвна зрнца (енг. red blood cell – RBC) или еритроцити су зреле, високо диференциране ћелије које се састоје само од ћелијске мембране коју чини двоструки слој липида и цитоплазме, унутрашње течности са издвојеним протеином хемоглобином. Хемоглобин боји еритроците у црвено. Концентрација црвених крвних зрнаца у крви варира од 3,7 до 5,8 1210 по литру. Пречник еритроцита износи приближно m8 у недеформисаном стању. Пошто црвена крвна зрнца чине највећи проценат од укупног броја ћелија у крви, боја крви потиче управо од ових ћелија. Црвена крвна зрнца током свог кретања кроз кардиоваскуларни систем пролазе кроз потпуно различите услове струјања крви, од високог Reynolds-овог броја у аорти, до веома малог Reynolds-овог броја и великог смичућег струјања у капиларним крвним судовима. Када се разматра струјање крви кроз велике крвне судове, чији је пречник већи од m200 , хидродинамичке силе које делују на ћелије су велике и зато ћелијска мембрана треба да буде довољно крута и да не подлегне њиховом утицају. Са друге стране, у малим крвним судовима пречника од 4 до m100 , долази до изражаја деформабилност ћелија. Особине мембране као што су флексибилност и велика способност деформације посебно су потребне када ћелија треба да прође кроз уске капиларне крвне судове чији је пречник истог реда величине као сама ћелија или чак и до два пута мањи, а притом могу имати и додатна сужења. Просечно време за које једно црвено крвно зрнце прође кроз цео организам је један минут, што значи да се све ове промене облика и услова у којима се ћелија креће дешавају веома брзо. Због свега тога, црвена крвна зрнца морају да имају способност да се брзо и лако деформишу са једне стране, а да са друге стране могу да задрже своју структуру и да се после деформисања врате у почетно стање. Црвено крвно зрнце је тзв. биконкавног облика. Тај облик изгледа као да је једна велика сфера са две стране у правцу исте осе притиснута са две сфере мањег пречника, услед чега је површина мембране конкавна. Због оваког облика црвено крвно зрнце има око 50% већу површину од сферне честице еквивалентног пречника. Више детаља о облику еритроцита дато је у поглављу 7.4.2. Услед саме структуре ових ћелија и услова којима су оне изложене током кретања, јављају се четири различита отпора деформацијама, који су детаљно анализирани са теоријског аспекта у поглављу 3. У литератури су измерене вредности параметара којима се дефинишу ови отпори [87,88]. У овом раду су за моделирање понашања црвених крвних зрнаца током кретања усвојене вредности у складу са онима које су предложене у поменутој литератури, а уједно и са вредностима које су користили Dupin и други [25]. Коефицијент отпора запреминској деформацији износи 150  mpNKV  . Коефицијент отпора локалној промени површине износи Докторска дисертација Поглавље 7 Тијана Ђукић 91 111067,1  mpNK Al  , док коефицијент отпора промени укупне површине износи 167,1  mpNK At  . Коефицијенти који су коришћени приликом дефинисања Skalak-овог материјалног модела ћелијске мембране износе 131075.0  mpNks  и 131075  mpNk  , према литератури [37]. Коефицијент отпора савијању мембране износи 1110  mpNK B  . У наставку овог поглавља ће бити приказано неколико примера моделирања понашања честице сферног облика и облика црвеног крвног зрнца. Наведени коефијенти отпора ће бити коришћени у оба случаја, тако да се сферна честица заправо понаша исто као крвно зрнце у датим условима. Добијени резултати ће бити поређени са аналитичким и резултатима презентованим у литератури. 7.4.1 Моделирање понашања сферне честице под дејством смичућег профила брзине У овом примеру, домен флуида је моделиран тако да су задате брзине горњег и доњег зида и ове брзине су супротног смера, као што је приказано на слици 7.15. Димензије домена у сва три правца су једнаке - zyx LLL  . На границама домена у правцу x и z осе задат је периодични гранични услов. Полупречник сферне честице једнак је једној десетини димензије домена флуида - xLr 10 1  . Честица је постављена тако да на почетку симулације буде у центру домена флуида, па су координате центра једнаке: 2 x C L x  ; 2 x C L y  ; 2 x C L z  Ако су вредности брзина горњег и доњег зида означене са wu и ако је висина домена означена са zL као на слици 7.15, онда се може дефинисати и коефицијент смицања: z w L u2  (7.11) Најпре је претпостављено да су густине флуида унутар и ван сферне честице једнаке    outin , што важи и за кинематске вискозности. У првом случају је дакле међусобни однос вискозности једнак јединици: 1 out in    (7.12) Докторска дисертација Поглавље 7 Тијана Ђукић 92 док ће у наставку бити разматрани и случајеви са различитим коефицијентом  . Reynolds-ов број се може дефинисати као:   2 Re r  (7.13) где је r полупречник сферне честице. Reynolds-ов број је у овом примеру веома мали, много мањи од један и у LB симулацијама је коришћена вредност 0.025, у складу са истим примером који је анализиран у литератури [36,42]. Слика 7.15 – Пример 1 – Шема геометрије домена флуида Постоји још један параметар који је потребно дефинисати, а то је бездимензиони коефицијент смицања: sk r G   (7.14) где је sk коефицијент отпора на смичућу деформацију у равни мембране, који је већ дефинисан у поглављу 3.1. У зависности од вредности овог параметра, мења се и деформација честице под дејством флуида. Докторска дисертација Поглавље 7 Тијана Ђукић 93 Да би се најбоље квантификовала деформација честице, користи се Taylor-ов параметар облика, који се дефинише на следећи начин: SS SS xy BL BL D    (7.15) где су са SL и SB означене дужа и краћа полуоса елипсе, која се добија као пресек сферне честице и равни која је паралелна са координатном равни yx0 и која пролази кроз тежиште честице. Полуосе елипсе приказане су на слици 7.16. Слика 7.16 – Пресек сферне честице у току деформисања и полуосе добијене елипсе Параметри који су подешени у оквиру LB симулације су дати у табели 7.1. У складу са примером из литературе са којим су поређена решења [36,42], параметри су одабрани тако да вредност времена релаксације  буде приближно једнака јединици. Сви остали потребни параметри су из физичких конвертовани у своје бездимензионе вредности, применом поступка који је описан у поглављу 7.1. Табела 7.1 – Подаци за LB солвер за пример 1 lbN 50 wu 0.004 Re 0.025 Током симулације се честица постепено деформише, док не достигне коначан степен деформације, када почиње ротационо кретање око тежишта. Поље брзине флуида и положај солида у том тренутку, када више нема додатне деформације честице приказани су на слици 7.17, док су на слици 7.18 приказане струјне линије, одакле се јасно види како флуид опструјава око честице која ротира око своје осе. Докторска дисертација Поглавље 7 Тијана Ђукић 94 Слика 7.17 – Пример 1 – Поље брзине флуида и тренутни положај сферне честице Слика 7.18 – Пример 1 – Струјне линије Варирањем параметра G , добијају се различите вредности Taylor-овог параметра облика xyD . На слици 7.19 је приказан коначни облик сферне честице за различито G . Облик који је приказан испрекиданом линијом је добијен применом LB солвера, док је Докторска дисертација Поглавље 7 Тијана Ђукић 95 плавом линијом приказан облик који је добијен у литератури [42] и као што се види, добијено је добро поклапање. Слика 7.19 – Пример 1 – Поређење коначног облика сферне честице у LB симулацији (црвена линија) са резултатима из литературе [42] (плава линија) На слици 7.20 приказана је промена Taylor-овог параметра облика у току времена за 4 различите вредности параметра G . Тачке представљају вредности добијене симулацијом применом LB солвера, а линије су добијене применом других метода у литератури [42]. Време је на дијаграму приказано као бездимензиона величина и израчунато је применом формуле:  itertd (7.16) где је iter тренутна итерација LB солвера. Слика 7.20 – Пример 1 – Промена Taylor-овог параметра облика током времена, за различите вредности параметра G (пуне линије представљају решења из литературе [42], а тачке решења LB симулације) Докторска дисертација Поглавље 7 Тијана Ђукић 96 Такође, угао који већа полуоса елипсе заклапа са осом x (видети слику 7.16) мења се током времена и током деформације честице. Вредности овог угла такође зависе од вредности параметра G и њихову зависност су испитивали Ramanujan и Pozrikidis [18]. На слици 7.21 су поређена њихова решења са резултатима добијеним применом LB солвера. Слика 7.21 - Пример 1 – Промена угла закошења честице током времена, за различите вредности параметра G (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) У истом раду су Ramanujan и Pozrikidis [18] разматрали утицај флуида различите вискозности који се налази унутар честице. Наиме, однос вискозности флуида унутар и ван честице који је дефинисан једначином (7.12) варира између 5 и 10, када је у питању црвено крвно зрнце уроњено у крвну плазму, док у случају када је црвено крвно зрнце уроњено у стандардни медијум у експериментима овај однос варира између 0.1 и 1. За различите вредности параметра  , деформација честице се разликује. На слици 7.22 је приказана промена Taylor-овог параметра облика у току времена за 4 различите вредности параметра G , за случај када је 2.0 , односно када је вискозност унутрашњег флуида 5 пута мања од вискозности спољашњег флуида. Иста промена је приказана на слици 7.23, али је сада 5 , односно вискозност унутрашњег флуида је 5 пута већа од вискозности спољашњег флуида. Треба приметити да се у овом случају решење добијено применом LB методе мало разликује од решења добијеног у литератури. То је последица чињенице да је у овом раду узет у обзир утицај промене 4 параметра деформације, као што је објашњено у поглављу 3, док су у литератури узета у обзир само површинска напрезања мембране. Докторска дисертација Поглавље 7 Тијана Ђукић 97 Слика 7.22 – Пример 1 - Промена Taylor-овог параметра облика током времена, за различите вредности параметра G, 2.0 (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) Слика 7.23 - Пример 1 - Промена Taylor-овог параметра облика током времена, за различите вредности параметра G, 5 (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) Да би се боље илустровао утицај вискозности унутрашњег флуида на деформацију сферне честице, на слици 7.24 је приказан тренутни облик честице за случај када је Докторска дисертација Поглавље 7 Тијана Ђукић 98 1.0G , за све три разматране вредности параметра  . Види се да је највећа деформација остварена када је 5 (линија црвене боје), а најмања када је 2.0 (линија плаве боје), док је зеленом бојом обележена деформација за 1 . Слика 7.24 – Пример 1 – Деформација честице за G=0.1, за различите вредности коефицијента  : црвена линија - 5 ; зелена линија - 1 ; плава линија - 2.0 Слика 7.25 - Пример 1 - Промена Taylor-овог параметра облика током времена, за случај растерећења честице Оно што је такође битно када се разматрају деформације честица које су хипереластичне и које се понашају као црвена крвна зрнца јесте способност повратка у почетно стање након престанка спољашњих утицаја који доводе до деформисања. У литератури [42] је овакво понашање тестирано управо на сферној честици, која је изложена утицају смичућег профила брзине. После неког времена (у бездимензионим јединицама, то време износи 1dt ), горњи и доњи зид домена престају да се крећу и анализира се понашање честице. Поново се прати промена Taylor-овог параметра Докторска дисертација Поглавље 7 Тијана Ђукић 99 облика у току времена. У овој симулацији вредност параметра G износи 0.075, а Reynlods-ов број износи 0.25. На слици 7.25 приказана је ова промена, при чему су резултати из литературе приказани пуном линијом, а резултати добијени применом LB солвера су приказани црвеним тачкама. Као што се може видети, добијено је добро поклапање решења, што показује да развијени софтвер може добро да моделира и овај вид понашања деформабилних честица. После извесног времена (у бездимензионим јединицама, то време износи 5dt ), честица се враћа у потпуно исто стање као на почетку симулације, што је илустровано на слици 7.26. Слика 7.26 - Пример 1 – Деформација честице током растерећења за G=0.075; плава линија - 0dt ; зелена линија - 1dt ; црвена линија - 5dt 7.4.2 Моделирање понашања црвеног крвног зрнца под дејством смичућег профила брзине Црвено крвно зрнце је биконкавног облика, пречника од 6,7 до 7,7 m , просечне дебљине од 1,7 до 2,5 m и запремине од 76 до 96 fL . Облик ових ћелија је прилагођен тако да оне могу лако да прођу и кроз најмање крвне судове (капиларе) без промене своје првобитне запремине или површине. Једначину којом се описује облик црвеног крвног зрнца емпиријски су дефинисали Evans и Fung [89]: 222 ryx                           2 2 0 22 22 0 22 10 5.0 2 0 22 0 1 2 1 R yx C R yx CC R yx Rz (7.17) где су вредности константи: mR 91.30  ; 207.00 C ; 0.21 C ; 123.12 C (7.18) а са r је означен полупречник црвеног крвног зрнца. Докторска дисертација Поглавље 7 Тијана Ђукић 100 Модел црвеног крвног зрнца који је дискретизован на троуглове и који је коришћен у симулацијама у овом раду је приказан на слици 7.27. На истој слици је приказан и попречни пресек еритроцита у равни која је паралелна zx0 равни и која пролази кроз тежиште крвног зрнца. Слика 7.27 – Дискретизован модел и попречни пресек црвеног крвног зрнца У овом примеру је домен флуида дефинисан на потпуно исти начин као у примеру са сферном честицом (у поглављу 7.4.1). Све вредности свих параметара LB симулације су исте. На исти начин као у претходном примеру дефинисан је и параметар G и у овој симулацији његова вредност је износила 0.2. И у овом случају праћена је промена Taylor-овог параметра облика током времена и резултати добијени применом LB солвера су поређени са резултатима које су публиковали Ramanujan и Pozrikidis [18]. Добро поклапање решења је добијено за однос вискозности унутрашњег и спољашњег флуида 1 (што је приказано на слици 7.28), као и за однос вискозности 5 (што је приказано на слици 7.29). Слика 7.28 - Пример 2 - Промена Taylor-овог параметра облика током времена, 1 Докторска дисертација Поглавље 7 Тијана Ђукић 101 Слика 7.29 - Пример 2 - Промена Taylor-овог параметра облика током времена, 5 Слика 7.30 - Пример 2 – Промена угла закошења црвеног крвног зрнца током времена, за различите вредности почетног угла (пуне линије представљају решења из литературе [18], а тачке решења LB симулације) Поред промене облика и деформације, црвено крвно зрнце у овом случају током времена обавља и ротационо кретање око своје осе. У зависности од почетног угла који већа „полуоса“ црвеног крвног зрнца заклапа са x осом, промена угла ротације се разликује. У симулацијама је однос вискозности унутрашњег и спољашњег флуида Докторска дисертација Поглавље 7 Тијана Ђукић 102 постављен на 5 , као што је реални случај код црвених крвних зрнаца. Резултати добијени применом LB солвера су поређени са резултатима из литературе [18] и добијено је добро поклапање решења, што се може видети на слици 7.30. У току кретања црвених крвних зрнаца у домену флуида под утицајем смичућег профила брзине уочава се још једна појава. Црвено крвно зрнце може или да ротира око своје осе без већих деформација (енг. tumbling behavior) или да се црвено крвно зрнце налази под одређених углом, а да притом мембрана ротира око тежишта (енг. tank- threading behavior), у зависности од вредности тзв. капиларног броја [90]. Капиларни број Ca се дефинише као однос вискозног напона и напона који се јавља на ћелијској мембрани, а напон на ћелијској мембрани се повезује са коефицијентом еластичности мембране. s out k r Ca    (7.19) Капиларни број је заправо једнак бездимензионом коефицијенту смицања G који је дефинисан једначином (7.14). Krüger и други [90] су закључили да ако је вредност овог параметра испод критичне, крвно зрнце ће ротирати око своје осе, без већих деформација, док ако је вредност овог параметра изнад критичне, онда ће крвно зрнце бити фиксирано, а мембрана ће ротирати око тежишта. Такође су израчунали да је критична вредност између 0.2 и 0.3. Промена облика крвног зрнца током времена када је 1.0Ca приказана је на слици 7.31, док је на слици 7.32 приказана промена облика крвног зрнца током времена када је 5.0Ca . На обе слике плава линија означава резултате које су добили Krüger и други [90], док црвена испрекидана линија представља резултате које су добијени применом LB солвера. Слика 7.31 - Пример 2 – Промена облика црвеног крвног зрнца током времена, када је 1.0Ca ; плава линија – резултати из литературе [90]; црвена линија – LB симулација Слика 7.32 - Пример 2 – Промена облика црвеног крвног зрнца током времена, када је 5.0Ca ; плава линија – резултати из литературе [90]; црвена линија – LB симулација Докторска дисертација Поглавље 7 Тијана Ђукић 103 Струјне линије за оба разматрана случаја приказане су на слици 7.33 (када је 1.0Ca ) и на слици 7.34 (када је 5.0Ca ). Слика 7.33 - Пример 2 – Струјне линије за случај када је 1.0Ca Слика 7.34 – Пример 2 – Струјне линије за случај када је 5.0Ca Докторска дисертација Поглавље 7 Тијана Ђукић 104 7.4.3 Моделирање кретања сферне честице у правом каналу У овом примеру је разматрано кретање сферне честице кроз прави канал, квадратног попречног пресека. Гранични услови су постављени тако да су брзине на зидовима једнаке нули, а да је на улазу и излазу задата разлика притисака, која доводи до кретања флуида. Геометрија домена флуида је приказана на слици 7.35. Дужине домена у правцу све три осе су дефинисане тако да важе следећи услови: zy LL  и yx LL 2 . Полупречник сферне честице једнак је једној четвртини ширине канала yLr 4 1  . Честица је постављена тако да на почетку симулације буде у центру правог канала, па су координате центра једнаке: 2 y C L x  ; 2 y C L y  ; 2 y C L z  . Слика 7.35 - Пример 3 – Шема геометрије домена флуида Параметри који су подешени у оквиру LB симулације су дати у табели 7.2 и усклађени су са примером из литературе са којим су поређена решења [42]. Сви остали потребни параметри су из физичких конвертовани у своје бездимензионе вредности, применом поступка који је описан у поглављу 7.1. Постављено је да су вискозности унутрашњег и спољашњег флуида у овом примеру једнаке. На слици 7.36 приказано је поље брзине флуида и положај солида у једном тренутку симулације применом LB солвера, док је на слици 7.37 приказана промена облика сферне честице током времена, која је упоређена са резултатима из литературе [42]. Докторска дисертација Поглавље 7 Тијана Ђукић 105 Табела 7.2 – Подаци за LB солвер за пример 3 lbN 50 wu 0.02 Re 0.71 Слика 7.36 – Пример 3 - Поље брзине флуида и тренутни положај сферне честице Слика 7.37 – Пример 3 – Промена облика сферне честице током времена, која је у центру правог канала у почетном тренутку; плава линија – резултати из литературе [42]; црвена линија – LB симулација У литератури [42] је анализирана и промена облика сферне честице, уколико се она помери за одређену малу вредност у односу на центар правог канала. Релативна позиција центра честице је дефинисана као:              22 yy Cr LL yy (7.20) Докторска дисертација Поглавље 7 Тијана Ђукић 106 Наравно, у случају који је прво разматран, сферна честица је постављена у центар правог канала, па је 0ry . Сферна честица је у литератури померана за 063.0ry и 13.0ry . Симулације са помереном честицом су урађене и применом LB солвера и добијено је добро поклапање решења, као што се може видети на сликама 7.38 и 7.39. Слика 7.38 - Пример 3 – Промена облика сферне честице током времена, која је у почетном тренутку померена за 063.0ry ; плава линија – резултати из литературе [42]; црвена линија – LB симулација Слика 7.39 - Пример 3 – Промена облика сферне честице током времена, која је у почетном тренутку померена за 13.0ry ; плава линија – резултати из литературе [42]; црвена линија – LB симулација 7.4.4 Моделирање кретања црвеног крвног зрнца у правом каналу Кретање црвеног крвног зрнца у правом каналу је моделирано у литератури [91]. У овом примеру су задати идентични гранични услови што се флуида тиче као и у примеру 7.4.3, с тим што је овде задата пет пута већа вискозност флуида унутар крвног зрнца, као што је случај и у реалним условима. С обзиром да је у примеру који је обрађен у литератури задат пречник канала m10 и максимална брзина флуида на улазу smm5.0 , исти услови су задати и у LB симулацији. Симулације су спроведене за два различита почетна положаја крвног зрнца. Први положај је када је крвно зрнце постављено вертикално, односно тако да је угао између најдуже осе крвног зрнца и x осе једнак 90 степени, а други положај је када је овај угао једнак 45 степени. На слици 7.40 је приказано поље флуида и положај крвног зрнца у току симулације, за први поменути положај, док је на слици 7.41 приказано поље флуида у току Докторска дисертација Поглавље 7 Тијана Ђукић 107 симулације, за други положај. Промена облика црвеног крвног зрнца током времена која је добијена применом LB солвера у оба случаја упоређена је са резултатима из литературе и добијена су добра поклапања решења, као што се види на сликама 7.42 и 7.43. Слика 7.40 - Пример 4 - Поље брзине флуида и тренутни положај крвног зрнца – Први случај Слика 7.41 - Пример 4 - Поље брзине флуида и тренутни положај крвног зрнца – Други случај Докторска дисертација Поглавље 7 Тијана Ђукић 108 Слика 7.42 - Пример 4 – Промена облика крвног зрнца током времена – Први случај; плава линија – резултати из литературе [91]; црвена линија – LB симулација Слика 7.43 - Пример 4 – Промена облика крвног зрнца током времена – Други случај; плава линија – резултати из литературе [91]; црвена линија – LB симулација 7.4.5 Моделирање кретања сферне честице кроз артерију са сужењем Атеросклероза је прогресивна болест артерија коју карактерише скупљање липида, холестерола, калцијума и ћелијских елемената унутар зида крвног суда. Овај процес узрокује појаву плака, промену попречног пресека крвног суда, што доводи до поремећаја протока крви. Разлика између изгледа нормалне артерије и артерије са стенозом приказана је на слици 7.44. Управо због појаве плака, неки крвни судови могу бити знатно сужени у односу на свој почетни пречник, па је приликом проласка честица кроз овакав крвни суд потребно да се честице додатно деформишу. Управо због тога је један овакав пример разматран и моделирано је понашање честице када пролази кроз сужење. Докторска дисертација Поглавље 7 Тијана Ђукић 109 Слика 7.44 – Шематски приказ нормалне и артерије са стенозом У овом примеру ће бити моделирано кретање сферне честице кроз једну артерију са стенозом. Шематски приказ геометрије домена флуида дат је на слици 7.45, при чему је крвни суд кружног попречног пресека. Дужине домена у правцу све три осе су дефинисане тако да важе следећи услови: zy LL  и yx LL 3 . Сферна честица је постављена тако да су на почетку симулације координате центра једнаке: 5 2 y C L x  ; 2 y C L y  ; 2 y C L z  . Пречник честице је 1.2 пута већи од пречника сужења крвног суда. Гранични услов је дефинисан тако да су брзине на зидовима једнаке нули, а почетни услов је задат као параболични профил брзине на улазу. Параметри који су подешени у оквиру LB симулације су дати у табели 7.3. Постављено је да су вискозности унутрашњег и спољашњег флуида у овом примеру једнаке. Табела 7.3 – Подаци за LB солвер за пример 5 lbN 50 wu 0.01 Re 1 Докторска дисертација Поглавље 7 Тијана Ђукић 110 Слика 7.45 – Пример 5 – Шема геометрије домена флуида На слици 7.46 приказано је поље брзине флуида и положај солида у почетном тренутку симулације применом LB солвера, док је на слици 7.47 приказано поље притиска такође у почетном тренутку симулације. Слика 7.46 – Пример 5 - Поље брзине флуида и почетни положај сферне честице Слика 7.47 - Пример 5 - Поље притиска флуида и почетни положај сферне честице Докторска дисертација Поглавље 7 Тијана Ђукић 111 Liu и други [26] су симулирали овај пример применом методе коначних елемената и одређена је промена пада притиска између улаза и излаза током трајања симулације. Исти пад притиска је одређен и у LB симулацијама и добијено је добро поклапање решења, као што се може видети на слици 7.48. Слика 7.48 – Пример 5 – Поређење промене пада притиска током времена добијене применом LB солвера са резултатима из литературе [26] Процес деформације честице током времена приказан је на слици 7.49. Приликом проласка кроз сужење, честица мора значајно да промени свој облик, с обзиром да је пречник суженог крвног суда мањи од пречника честице. Након проласка кроз сужење, честица се враћа у облик сличан почетном, с тим што непосредно после изласка има облик падобрана, са увученим делом ближе средини, који се касније постепено враћа у почетни облик. t = 0.1s Докторска дисертација Поглавље 7 Тијана Ђукић 112 t = 0.3s t = 0.4s t = 0.45s t = 0.5s Докторска дисертација Поглавље 7 Тијана Ђукић 113 t = 0.6s t = 0.7s t = 1.25s Слика 7.49 – Пример 5 - Симулација проласка сферне честице кроз артерију са стенозом Докторска дисертација Поглавље 7 Тијана Ђукић 114 7.4.6 Моделирање кретања црвеног крвног зрнца и канцер ћелије кроз микрофлуидни чип за сепарацију Проучавање динамике микрофлуида почело је још осамдесетих година прошлог века. Од тада, примена уређаја који проучавају понашање флуида на микроскопском нивоу је значајно проширена [92]. Последњих година почели су да се израђују тзв. микрофлуидни чипови (енг. microfluidic chip). Овакви уређаји се примењују у многим областима, пре свега молекуларне биологије за анализу различитих особина ћелија [93]. Компјутерске симулације могу да буду од велике помоћи у току пројектовања микрофлуидних чипова, јер могу да смање време пројектовања и цену, пошто се све анализе могућег понашања ћелија могу компјутерски предвидети и анализирати пре реалних експеримената. Применом нумеричких модела, као што је модел који је представљен у овом раду, могуће је моделирати кретање различитих објеката чије се понашање прати у оквиру микрофлуидног чипа, као што је на пример праћење кретања крвних ћелија кроз крв, кретања бактерија кроз воду, других микроорганизама у течности и слично. Како је тема овог рада примена нумеричких симулација у биомедицинском инжењерингу, у наставку ће бити приказано како се развијени нумерички модел може применити на анализу микрофлуидног чипа за сепарацију канцер ћелија. Канцер је тип болести која се карактерише неконтролисаним повећањем броја ћелија, које временом почињу да се шире и у околна ткива. Такође, одређени број канцер ћелија, тзв. циркулишуће ћелије тумора (енг. circulating tumour cells) се одвајају од примарног тумора, доспевају у крв и на тај начин даље настављају да се шире по организму и да метастазирају у другим органима. Чак и ако се примарни тумор одстрани, ове ћелије остају у организму. Циркулишуће ћелије тумора се могу идентификовати у крви пацијента узимањем узорка. Додатна биолошка анализа ових ћелија може значајно да допринесе лакшем дијагностиковању болести, али и одређивању адекватне терапије специфично прилагођене и пацијенту и типу тумора. Међутим, ове ћелије се јако ретко јављају – у просеку се јавља једна циркулишућа ћелија тумора на 5-10 милиона црвених крвних зрнаца. Због тога је потребно применити неку савременију методу за изолацију ових ћелија из крвног узорка, имајући притом у виду карактеристике ћелија. У ову сврху се примењују претходно поменути микрофлуидни чипови, у којима се обавља сепарација ћелија. На слици 7.50 је приказан микроскопски снимак једног оваквог микрофлуидног чипа. Микрофлуидни чип заправо представља једну минијатурну комору, која је испресецана преградама, између којих струји крв. Једна могућност је да се у оквиру микрофлуидног чипа поставе преграде које су прекривене антителима, које ће да привлаче канцер ћелије и да их задржавају [94]. Друга могућност је да се у оквиру микрофлуидног чипа поставе фиксиране преграде на тачно утврђеном растојању, које ће да обезбеде филтрацију на основу величине ћелија [95,96]. У наставку ће бити даље анализиран овај други приступ. Докторска дисертација Поглавље 7 Тијана Ђукић 115 Слика 7.50 – Микроскопски снимак микрофлуидног чипа [94] Варирањем размака између преграда у чипу, постиже се да одређене ћелије прођу између преграда, а да се друге задрже. Наиме, циркулишуће ћелије тумора су веће и чвршће од црвених крвних зрнаца, а мање од белих крвних зрнаца. Применом два различита филтера, могуће је најпре издвојити највеће ћелије (бела крвна зрнца), а затим оне најмање и веома деформабилне (црвена крвна зрнца). Да би се одредио најповољнији размак између фиксираних преграда, најлакше је приликом пројектовања микрофлуидног чипа применити нумеричке симулације. Битно је тачно и прецизно одредити ово растојање да би се постигао што већи проценат задржаних канцер ћелија. Нумерички модел који је тема овог рада разматра струјање флуида на микроскопском нивоу и моделира понашање појединачних ћелија у оквиру крвне плазме, па се може веома успешно применити у овом случају. Карактеристике циркулишућих ћелија тумора се наравно разликују у зависности од типа тумора. У овом случају је моделирано понашање једне ћелије канцера дојке. У литератури су спроведени експерименти за утврђивање карактеристика овог типа ћелије [97,98]. По угледу на рад објављен у литератури [99], једне истраживачке групе која се бави управо истраживањима у овој области и анализом микрофлуидних чипова за сепарацију канцер ћелија, у овом раду су коришћене исте карактеристике канцер ћелија. У табели 7.4 приказане су вредности параметара које су коришћене у овом раду. С обзиром да карактеристике канцер ћелија нису толико испитиване као карактеристике еритроцита, није познато да ли се Skalak-ов модел мембране може применити и на овај тип ћелија. Због тога је у овом примеру за моделирање површинског напрезања канцер ћелије коришћен једноставнији модел, представљен у поглављу 3.1, уместо моделирања применом методе коначних елемената. Такође, за разлику од црвеног крвног зрнца, које је биконкавног облика, канцер ћелија је сферног облика. Пречник канцер ћелије је два пута већи од пречника црвеног крвног зрнца. Докторска дисертација Поглавље 7 Тијана Ђукић 116 Табела 7.4 - Вредности параметара за канцер ћелију дојке Параметар Вредност VK 15.14 mpN AtK 148.0 mpN AlK 111048.0  mpN BK 111029.0  mpN SK 11104375.1  mpN Шематски приказ геометрије домена флуида дат је на слици 7.51. Дужине домена у правцу све три осе су дефинисане тако да важе следећи услови: zy LL  и yx LL 3 . И канцер ћелија (моделирана као сферна честица), као и црвено крвно зрнце су на почетку симулација постављени тако да су координате центра једнаке: 5 2 y C L x  ; 2 y C L y  ; 2 y C L z  . Пречник црвеног крвног зрнца је 1.35 пута већи од растојања између фиксираних препрека, док је пречник канцер ћелије, као што је већ поменуто, 2 пута већи од пречника црвеног крвног зрнца. Гранични услов је дефинисан тако да су брзине на зидовима једнаке нули, а почетни услов је задат као разлика притисака између улаза и излаза. Параметри који су подешени у оквиру LB симулације су дати у табели 7.5. Постављено је да су вискозности унутрашњег и спољашњег флуида за канцер ћелију у овом примеру једнаке, док је за случај црвеног крвног зрнца постављено да је вискозност унутрашњег флуида пет пута већа од вискозности спољашњег флуида. Табела 7.5 – Подаци за LB солвер за пример 6 lbN 50 wu 0.01 Re 1 Докторска дисертација Поглавље 7 Тијана Ђукић 117 Слика 7.51 - Пример 6 – Шема геометрије домена флуида На слици 7.52 приказано је поље брзине флуида и положај канцер ћелије у почетном тренутку симулације применом LB солвера, док је на слици 7.53 приказано поље притиска и положај црвеног крвног зрнца такође у почетном тренутку симулације. Слика 7.52 - Пример 6 - Поље брзине флуида и почетни положај канцер ћелије Докторска дисертација Поглавље 7 Тијана Ђукић 118 Слика 7.53 - Пример 6 - Поље притиска флуида и почетни положај црвеног крвног зрнца Процес деформације сферне честице, односно канцер ћелије, после одређеног броја итерација приказан је на слици 7.54. Као што се види, ова ћелија нема довољно велику способност деформације, тако да остаје заглављена у пролазу и не може да прође између две преграде. За разлику од ње, црвено крвно зрнце је веома деформабилно и лако пролази између преграда, као што се види на слици 7.55, где је приказано његово кретање кроз микрофлуидни чип. Управо ово и јесте циљ, да црвена крвна зрнца прођу између преграда, а да канцер ћелије остану заглављене и да се на тај начин изврши жељена сепарација. после 1000 итерација Докторска дисертација Поглавље 7 Тијана Ђукић 119 после 3000 итерација после 5000 итерација после 7000 итерација Слика 7.54 - Пример 6 - Симулација проласка канцер ћелије кроз микрофлуидни чип Докторска дисертација Поглавље 7 Тијана Ђукић 120 после 3000 итерација после 5000 итерација после 6000 итерација Докторска дисертација Поглавље 7 Тијана Ђукић 121 после 7000 итерација после 8000 итерација Слика 7.55 - Пример 6 - Симулација проласка црвеног крвног зрнца кроз микрофлуидни чип Тијана Ђукић 122 8 Поређење решења нумеричких симулација са експерименталним резултатима У овом поглављу ће бити приказано поређење резултата који су добијени нумеричким симулацијама применом LB солвера са експерименталним резултатима. Сви параметри симулација ће бити подешени тако да у потпуности одговарају условима који су дефинисани у експериментима. Кроз ове примере ће бити доказано да се развијена нумеричка метода и развијени софтвер могу успешно користити за моделирање реалних проблема у биомедицинском инжењерингу. Најпре ће бити анализирано кретање црвених крвних зрнаца кроз стаклену цев, што представља један пример експеримента који је спроведен in-vitro, у потпуно лабораторијским условима. Потом ће бити анализирано кретање црвених крвних зрнаца кроз бифуркацију пацова, која је извађена из тела пацова и држана је у одговарајућим лабораторијским условима. Ово представља један пример експеримента који је спроведен ex-vivo, у условима који су врло слични условима у реалном организму. На крају ће бити спроведена анализа кретања црвених крвних зрнаца кроз сплет репних вена у живој зебра риби, што представља један експеримент in-vivo, где је посматран проток крви кроз реалан организам у реалном времену. 8.1 Кретање црвених крвних зрнаца кроз стаклену цев in-vitro Овај експеримент су спровели Secomb и други [100] и у оквиру овог експеримента посматрали су понашање људских црвених крвних зрнаца унутар стаклене цеви. Пречник стаклене цеви је m7 . Цео процес је наравно посматран под микроскопом. Крв је узета од волонтера, спроведене су одговарајуће процедуре почетне припреме крвних ћелија, као што је објашњено у литератури [100] и потом је крв пуштана из шире цеви у поменуту цев пречника m7 , тако да је на крајевима цеви постојао фиксиран пад притиска и тако да је брзина црвених крвних зрнаца била приближно једнака smm1 . На слици 8.1 приказан је микроскопски снимак проласка црвених крвних зрнаца кроз стаклену цев, при чему се проток дешава са лева на десно. Ови резултати су меродавни и зато што је добијен изглед који је врло сличан реалном, што се може видети на слици 8.2. На овој слици је приказан микроскопски снимак протока крви кроз крвни суд Докторска дисертација Поглавље 8 Тијана Ђукић 123 пацова. Проток се и овде дешава са лева на десно. Ова слика је добијена у оквиру експеримента који су спровели Pries и Secomb [9]. Слика 8.1 – Микроскопски снимак проласка црвених крвних зрнаца кроз стаклену цев [9] Слика 8.2 – Микроскопски снимак струјања крви кроз крвни суд пацова [9] Три различите симулације су спроведене применом LB солвера, са различитим почетним положајем црвеног крвног зрнца. Сви услови су задати на идентичан начин као у експерименту – вискозност крвне плазме је постављена на scm 2018,0 , у складу са [101]; задата је максимална брзина на улазу од smm1 ; шира стаклена цев је на почетку пречника m14 , а потом се сужава до експериментом фиксираних m7 . Пречник црвеног крвног зрнца је m8 . У првој симулацији је крвно зрнце постављено на средину цеви, у другој је померено за m5.0 на горе, док је у трећој померено за m1 на горе. Како се крвно зрнце креће кроз цев, оно постепено мења облик, тако да постаје конвексно на почетном делу, а конкавно на задњем делу. Такође, оне ћелије које у почетном тренутку нису биле постављене у центру цеви, временом су добиле асиметричан облик, издужен на оној страни која је ближе зиду. На слици 8.3 су упоредо приказани облици крвног зрнца који су добијени нумеричком симулацијом применом LB солвера и облици који су добијени експериментом, у сва три случаја. Облици добијени симулацијом се поклапају са облицима који су добијени експериментом. Докторска дисертација Поглавље 8 Тијана Ђукић 124 Слика 8.3 – Упоредни приказ облика црвеног крвног зрнца добијених нумеричком симулацијом применом LB солвера (слике лево) и експериментом (слике десно); релативна почетна позиција крвног зрнца у односу на центар стаклене цеви одозго на доле - 00 y ; my 5.00  ; my 10  На слици 8.4 приказано је поље притиска добијено применом LB солвера, као и крајњи положај црвеног крвног зрнца које је у почетку постављено на средину цеви. На сликама 8.5 и 8.6 приказано је поље брзине добијено применом LB солвера, као и крајњи положај црвеног крвног зрнца које је у почетку померено у односу на средину цеви за m5.0 и m1 , респективно. Слика 8.4 – Поље притиска флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку постављено на средину стаклене цеви Докторска дисертација Поглавље 8 Тијана Ђукић 125 Слика 8.5 - Поље брзине флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку померено за m5.0 у односу на средину стаклене цеви Слика 8.6 - Поље брзине флуида и крајњи положај црвеног крвног зрнца које је у почетном тренутку померено за m1 у односу на средину стаклене цеви 8.2 Кретање црвених крвних зрнаца кроз каротидну бифуркацију пацова ex-vivo Овај експеримент су, као и претходни, спровели Secomb и други [9], с тим што је овај експеримент извођен ex-vivo на пацовима. Из тела пацова је извађен део абдомена, који је потом држан у одговарајућим условима и крв је циркулисала кроз артерије у оквиру тог сегмента, а проток крви је сниман тако да су добијени видео снимци и фотографије. Детаљи целокупног експерименталног поступка описани су у литератури [9,100]. Посебна пажња је посвећена артерији са бифуркацијом. На слици 8.7 приказана је једна артерија са гранањем и крвна зрнца која пролазе кроз њу. Докторска дисертација Поглавље 8 Тијана Ђукић 126 Слика 8.7 – Кретање црвених крвних зрнаца кроз артерију са бифуркацијом Применом LB солвера спроведена је симулација кретања црвеног крвног зрнца кроз геометријски сличну артерију са бифуркацијом. Шематски приказ геометрије домена флуида дат је на слици 8.8. Пречник улазне артерије је постављен на m12 , пречник излазних артерија је m10 . Улазна брзина износи smm1 , а брзина на зидовима артерије једнака је нули. Приликом моделирања црвеног крвног зрнца, симулирана је различита вискозност флуида, тако да је вискозност унутрашњег флуида била пет пута већа од вискозности спољашњег флуида. Слика 8.8 - Шема геометрије артерије са бифуркацијом Докторска дисертација Поглавље 8 Тијана Ђукић 127 На слици 8.9 приказан је положај крвног зрнца у LB симулацији, као и поље притиска флуида у почетном тренутку, када крвно зрнце почиње да се креће кроз артерију. Моделирани процес кретања црвеног крвног зрнца кроз артерију са бифуркацијом приказан је на слици 8.10. Слика 8.9 – Поље притиска флуида и почетни положај црвеног крвног зрнца почетни тренутак после 5 ms Докторска дисертација Поглавље 8 Тијана Ђукић 128 после 10 ms после 15 ms после 20 ms после 25 ms после 30 ms после 40 ms Слика 8.10 - Симулација проласка црвеног крвног зрнца кроз артерију са бифуркацијом Докторска дисертација Поглавље 8 Тијана Ђукић 129 Промена облика крвног зрнца илустрована је на слици 8.11. Добијено је добро поклапање облика са онима који су добијени у екперименту (што се може видети са слике 8.7). Слика 8.11 – Промена облика црвеног крвног зрнца током времена 8.3 Моделирање протока крви кроз кардиоваскуларни систем зебра рибе Процес ембрионалног стварања и каснијег функционисања кардиоваскуларног система је тешко проучавати in-vivo, применом класичне ембриолошке и генетске анализе [102]. Код кичмењака постоји развијен кардиоваскуларни систем, са срцем које пумпа крв и крвним судовима које имају своје јасно дефинисане зидове. Међутим, посматрање ових крвних судова у живим ембрионима је тешко или зато што се ембриони развијају унутар мајчине утробе или зато што нису довољно транспарентни. Врста риба праве кошљорибе (teleostei) и једна подврста која се назива зебра риба (zebrafish) имају велике предности наспрам других врста када је проучавање у питању. Пре свега, оне су малих димензија. Одрасла зебра риба најчешће може да порасте до 4 cm. Затим, веома су плодне и полажу велики број јаја. Размножавање је спољашње, ван тела мајке. Такође, зебра рибе се веома добро могу оптички анализирати. Због малих димензија зебра рибе, ембриони могу да преживе и ако им се екстерно додаје јако мала количина кисеоника. Иако постоје извесне варијације у детаљима анатомије кардиоваскуларног система код ових риба, основе овог система су сличне као и код других кичмењака. Због свега поменутог, у литератури се могу наћи бројни радови у којима су анализирани многи аспекти развоја кардиоваскуларног система код зебра риба, као и бројне генетске анализе мутација различитих гена [103,104]. Постоје и бројне студије које су проучавале настанак еритроцита како код кичмењака, тако и код кошљориба Докторска дисертација Поглавље 8 Тијана Ђукић 130 [105,106,107]. Коришћењем технике конфокалне микроангиографије [104] или применом зелених флуоресцентних протеина (енг. green fluorescent protein – GFP) који се генеришу посебним поступцима [108,109], при томе се ослањајући на чињеницу да су зебра рибе малих димензија и веома прозрачне, могуће је сликати крвне судове кроз целу дебљину рибе. На тај начин се посматра анатомија крвних судова и снима се проток крвних зрнаца кроз крвни суд. За разлику од неких других техника анализе васкулатуре, поменути приступ не угрожава ни на који начин саму рибу и слике које се добијају односе се на потпуно активну циркулацију крви. Слика 8.12 – Фазе ембрионалног развоја зебра рибе [110] Зебра риба (латински назив Danio rerio) је тропска врста која се најчешће може наћи у југоисточном делу Хималаја, у деловима Индије, Пакистана, Бангладеша, Непала и Бурме. Процес ембрионалног развоја из положеног јајета је веома брз и приказан је на слици 8.12. С обзиром да су јаја провидна, цео процес се лако може пратити. Овај процес је заправо врло сличан процесу развоја људског ембриона и управо зато и јесте интересантна ова врста риба, јер се анализом процеса код ове врсте могу извући важни закључци који се могу касније проширити и на људски организам. Термини на слици 8.12 нису превођени и слика је у оригиналу преузета са сајта [110], и ово је био део истраживања у оквиру Универзитета у Орегону у САД, као и једне истраживачке групе Универзитета Santa Barbara у Калифорнији. Сви унутрашњи органи, као што су очи, мозак, срце, унутрашње уво се развијају већ у прва 3 дана. Кардиоваскуларни систем је један од првих система који се развија током ембрионалног развоја зебра рибе. Већ Докторска дисертација Поглавље 8 Тијана Ђукић 131 после једног дана од тренутка када је ларва положена, зебра риба је на нивоу развоја када почиње да циркулише крв. Детаљи самог процеса развоја кардиоваскуларног система нису тема овог рада, а детаљно објашњење процеса настанка свих већих крвних судова у зебра риби може се наћи у литератури [102,111]. На слици 8.13 приказана је једна одрасла зебра риба. Ова слика је добијена применом зелених флуоресцентних протеина [108,109], тако да се лако могу видети главни крвни судови. Крвни притисак током кретања крви кроз кардиоваскуларни систем одрасле зебра рибе се креће од 0.8 mmHg у срчаним преткоморама, око 1-1.6 mmHg у главним артеријама, па до приближно 2 mmHg у срчаним коморама [112]. Брзина протока крви је измерена у рибама старим 3-5 дана и износила је око smm1 [113]. Слика 8.13 – Одрасла зебра риба и њен кардиоваскуларни систем Слика 8.14 – Кардиоваскуларни систем у 3 дана старој зебра риби [111] На слици 8.14 је мало детаљније објашњено где се који крвни суд налази у оквиру зебра рибе. Скраћенице су преузете из литературе [102] - AA означава лукове аорте (енг. Докторска дисертација Поглавље 8 Тијана Ђукић 132 aortic arches); CV је репна вена (енг. caudal vein); CCV је заједничка главна вена (енг. common cardinal vein); DA је леђна аорта (енг. dorsal aorta); PCV је задња главна вена (енг. posterior cardinal vein); PHS је примарни синусни чвор главе (енг. primary head sinus); SA је попречна (сегментна) артерија (енг. segmental artery); SV је попречна (сегментна) вена (енг. segmental vein); DLAV је леђни подужни крвни суд (енг. dorsal longitudinal anastromotic vessel). У овом раду посматрано је кретање крвних зрнаца кроз сплет репних вена (енг. caudal vein plexus), које су означене црвеном стрелицом на слици 8.15, која је преузета из литературе [114]. Слика 8.15 – Сплет репних вена зебра рибе [114] Симулације применом LB солвера су урађене за два случаја кретања крвног зрнца у венама, у риби чија је циркулација посматрана 32 сата после фертилизације (енг. hours post fertilization – hpf). Геометрија домена флуида је креирана директно на основу одговарајућих микроскопских снимака добијених током експеримената са зебра рибама. Ови експерименти су спроведени на Институту за анатомију, Универзитета у Берну, у Швајцарској, за потребе истраживања у оквиру пројекта на коме је ангажован и аутор овог рада. „Прорези“ који се могу видети на слици 8.15 су тзв. пилари (енг. pillars), који се формирају, током времена расту и у току каснијег развоја омогућавају даље гранање вена и стварање нових крвних судова. Ширина сплета вена износи m100 , а брзина на улазу је измерена на основу видео снимака добијених експериментално (заправо је рађено микроскопско сликање рибе сваке 2 ms и спајањем тих слика добијени су поменути снимци) и износи sm60 . На слици 8.16 приказани су резултати за први посматрани почетни положај крвног зрнца, за неколико временских тренутака. Лево су приказане микроскопске слике сплета вена зебра рибе, при чему је конкретно крвно зрнце означено жутим кругом, десно су приказани резултати добијени применом LB солвера, док су у средини издвојени само облици конретног крвног зрнца, ради лакшег поређења (плавом бојом је обележен облик добијен симулацијом, а црвеном бојом облик извучен са микроскопског снимка). На слици 8.17 приказани су исти резултати, само за други посматрани почетни положај крвног зрнца, такође за неколико временских тренутака. Као што се види са издвојених облика посматраног крвног зрнца, резултати симулација се веома добро слажу са експерименталним резултатима. Докторска дисертација Поглавље 8 Тијана Ђукић 133 почетни тренутак после 8 ms после 18 ms после 28 ms Докторска дисертација Поглавље 8 Тијана Ђукић 134 после 40 ms Слика 8.16 – Поређење резултата добијених из експеримента и LB симулације за први почетни положај крвног зрнца; лево – микроскопски снимак зебра рибе, са означеним конкретним крвним зрнцем; у средини – издвојени облици конкретног крвног зрнца (црвено – експеримент, плаво – симулација); десно резултати добијени применом LB симулације почетни тренутак после 5 ms Докторска дисертација Поглавље 8 Тијана Ђукић 135 после 15 ms после 21 ms Слика 8.17 - Поређење резултата добијених из експеримента и LB симулације за други почетни положај крвног зрнца; лево – микроскопски снимак зебра рибе, са означеним конкретним крвним зрнцем; у средини – издвојени облици конкретног крвног зрнца (црвено – експеримент, плаво – симулација); десно резултати добијени применом LB симулације Тијана Ђукић 136 9 Закључак Кардиоваскуларни систем је један од најважнијих система у људском организму. Начин на који крв струји кроз многобројне мале крвне судове у микроциркулацији људског организма је од великог значаја за исправно функционисање органа. Анализа процеса који се овде дешавају је веома битна због великог утицаја микроскопских појава на цео организам. Ови процеси се тешко експериментално испитују и због тога у овој области нумеричке симулације могу да буду јако корисне. Постоји велика примена оваквих нумеричких симулација у ћелијској биологији, хемијском инжењерингу, медицини и биомедицинском инжењерингу. Моделирање кретања деформабилних тела уроњених у флуид је најважнији корак ка моделирању много комплекснијих проблема, као што је филтрација канцер ћелија и испорука лекова кроз крвоток. Управо због тога је ова тема обрађивана у оквиру овог рада. У оквиру овог рада представљен је комплетан нумерички модел који успешно симулира тродимензионално струјање крви кроз крвне судове у људском организму. Ово је показано на неколико примера, где је моделирано струјање крви кроз људску аорту, каротидну и коронарну артерију, при чему је геометрија самих крвних судова била или параметарски задата или креирана на основу реалних клиничких података о пацијенту, што је додатна вредност развијеног софтвера. Такође, представљен је нумерички модел који симулира комплексно кретање деформабилних честица које су изложене значајним променама облика у капиларним крвним судовима. Поређењем резултата спроведених симулација са резултатима који су објављени у литератури, као и са експерименталним резултатима, показано је да се развијени софтвер може успешно примењивати у симулацијама реалних проблема у биомедицинском инжењерингу на микроскопском нивоу. Симулирано је укупно 9 различитих проблема, од којих су неки веома значајни у области биомедицине, као што је кретање црвених крвних зрнаца кроз бифуркацију или артерију са стенозом, као и кретање црвених крвних зрнаца и канцер ћелија кроз микрофлуидни чип за сепарацију канцер ћелија. Такође, у оквиру нумеричке методе омогућено је да се моделира кретање деформабилне честице, односно црвеног крвног зрнца, које је испуњено флуидом који је различитих карактеристика од околне крвне плазме, што је апсолутно идентично реалним условима у људском организму. С обзиром да један овакав софтвер за научне симулације захтева много рачунарских ресурса и да извршавање једне симулације траје доста дуго, приликом имплементације нумеричке методе и током развоја софтвера коришћена је посебна технологија коју је развила NVIDIA, тзв. CUDA архитектура. На тај начин је омогућено да се развијени софтвер може покретати на компјутерским кластерима са GPU уређајима, што омогућава драстично убрзање постојећих програма. У овом раду је применом CUDA архитектуре комплетан софтвер паралелизован и поређена је брзина рада секвенцијалне Докторска дисертација Поглавље 9 Тијана Ђукић 137 и паралелизоване верзије. На специјализованом Tesla уређају постигнуто је убрзање од чак 20 пута, док је у оквиру међународних PRACE пројеката приступано и кластерима са великим бројем GPU уређаја и постигнуто је и много пута веће убрзање. Симулације струјања флуида и солид-флуид интеракције, за које је раније требало и по неколико сати, сада се извршавају за неколико минута. То практично значи да се резултати симулација и кретање деформабилних честица сада могу интерактивно пратити у реалном времену. Крајњи циљ овог рада био је да се креира нумерички модел који симулира микроциркулацију у људском организму, до крвних судова величине микрометра, укључујући при томе комплетно симулирање динамике кретања ћелија. Већина модела који су до сада представљени у литератури су ограничени на симулације у геометријски једноставним крвним судовима. Метода која је представљена у овом раду и развијени софтвер могу да се користе за моделирање сложеног кретања вискоеластичних тела, као што су црвена крвна зрнца, у капиларним крвним судовима чија је геометрија комплексна. У овом раду је разматрано кретање једне честице кроз флуид. Уз врло мале адаптације, додавањем одговарајуће силе интеракције између крвних зрнаца, могуће је моделирати и кретање више честица истовремено. Такође, на сличан начин би могло да се моделира и истовремено кретање свих крвних конституената, што би допринело да нумеричке симулације буду још ближе реалним процесима у људском организму. Ово је свакако један од даљих праваца истраживања и усавршавања развијене методе и софтвера. Представљени нумерички модел и софтвер могу значајно да помогну у бољем разумевању и решавању бројних реалних проблема у биомедицинском инжењерингу, као и у унапређивању медицинских третмана, као што су на пример процеси стварања плака и тромба, понашање, кретање и филтрација канцер ћелија, као и моделирање свих других процеса који се дешавају у оквиру лабораторијских микрофлуидних чипова, затим анализа и унапређење процеса интраваскуларне дистрибуције и испоруке матичних ћелија, лекова и генетског материјала до циљаних органа. У свим овим процесима није тако једноставно и јефтино спроводити експерименте, а овакве нумеричке симулације нуде јако добру и ефикасну алтернативу. Због нумеричке методе која доста реално осликава услове у људском организму, велике тачности која је показана кроз поређење са експерименталним резултатима и бројним другим решењима из литературе, као и због брзог извршавања на компјутерским кластерима, које је омогућено паралелизацијом, нумеричка метода представљена у овом раду и развијени софтвер имају велики потенцијал када је реч о практичној примени у симулацијама реалних проблема у биоинжењерингу. Тијана Ђукић 138 Додатак А Boltzmann-ова једначина се трансформише у Навије-Стоксове једначине тако што се најпре уведе оператор I као интеграл:          ξξxξ x x dtfh t thI ,, , 1 ,   (А.1) где је  ξh интеграбилна функција, а  густина флуида. Једначина која се добија за 1h одговара маси, а једначина која се добија за ξh одговара количини кретања. Могу се израчунати вредности дефинисаног интеграла за одговарајуће вредности функције h [55]:   11 I ;   uξ I ;   uu P ξξ   I (А.2) где су са u и P означени вектор брзине и тензор напона, респективно. Ако се обе стране Boltzmann-ове једначине помноже са  ξh и интеграле по пољу брзине, добија се:    ξξ ξ g x ξ dhffdh ff t f                01  (А.3) Десна страна једнакости (А.2) је једнака нули, за наведене вредности функције h , пошто морају да важе закони одржања масе и количине кретања за идеалне сударе, а такви се овде разматрају. Докторска дисертација Додатак А Тијана Ђукић 139 Лева страна се даље може трансформисати (при томе се има у виду да функција h зависи само од брзине ξ и примењују се правила за извод производа):                                                                                        ξ ξ xgξξ x ξ ξ ξ gngξξ x ξ ξg ξ ξg ξ ξξξ x ξξ x ξ ξ ξ gξ x ξξ htIhIhI t dhfdfhdhfdhf t dhfdhfdhfdhfdhf t dh f dh f dh t f , 0 0           (А.4) Заменом   1ξh у горњој једначини, узимајући у обзир вредности дефинисане у (А.2), добија се једначина одржања масе (једначина континуитета):   0      u x   t (А.5) Заменом   ξξ h у горњој једначини, поново узимајући у обзир вредности дефинисане у (А.2), добија се једначина одржања количине кретања (Навије-Стоксова једначина):     guuP x u        t (А.6) Треба напоменути да је тензор напона непознат и потребно га је накнадно дефинисати, што је урађено у поглављу 2.2. Тијана Ђукић 140 Додатак Б Структура мреже са q различитих праваца брзина у d -димензионом простору се означава као DdQq мрежа (енг. lattice). За симулације струјања изотермалног нестишљивог флуида користи се неколико 2D и 3D структура мреже и они су овде наведени. За сваки тип је дата слика која илуструје правце брзина. Такође, дати су и тежински коефицијенти i и вектори који се користе за дефинисање апсцисних оса iξ . Треба напоменути да је за све наведене структуре коефицијент скалирања означен као „брзина звука мреже“ sc исти и износи: 3 12 sc 92QD          8,7,6,5 4,3,2,1 0 361 91 94 i i i i                 8,7,6,5 4,3,2,1 0 3,3 3,0;0,3 0,0 i i i iξ 153QD          14,,7 6,,1 0 721 91 92   i i i i                   14,,7 6,,1 0 3,3,3 3,0,0;0,3,0;0,0,3 0,0,0   i i i iξ Докторска дисертација Додатак Б Тијана Ђукић 141 193QD          18,,7 6,,1 0 361 181 31   i i i i                       18,,7 6,,1 0 3,0,3;3,3,0;0,3,3 3,0,0;0,3,0;0,0,3 0,0,0   i i i iξ 273QD             26,,15 14,,7 6,,1 0 541 2161 272 278    i i i i i                               26,,15 14,,7 6,,1 0 3,0,3;3,3,0;0,3,3 3,3,3 3,0,0;0,3,0;0,0,3 0,0,0    i i i i iξ Тијана Ђукић 142 Додатак Ц Већина савремених графичких картица, који су део конфигурација данашњих персоналних рачунара, могу се користити као GPU уређај. Карактеристике једне такве графичке картице NVIDIA GeForce GT 630 датe су у табели Ц.1. Табела Ц.1 – Карактеристике NVIDIA GeForce GT 630 CUDA верзија уређаја (енг. compute capability) 2.1 Укупна расположива глобална меморија 1024 MB Број мултипроцесора 2 Број CUDA језгара (8 по мултипроцесору) 96 Радни такт GPU уређаја 1.62 GHz Укупан број регистара на располагању једног блока 32768 Максималан број нити у једном блоку 1024 Максималне димензије блока 1024 x 1024 x 64 Максималне димензије мреже блокова 65535 x 65535 x 65535 Ипак, NVIDIA у својој понуди има и специјализоване GPU уређаје, такозвану Tesla серију уређаја. Они су примарно намењени за извршавање комплексних научних симулација. Основна разлика између класичне графичке картице и Tesla уређаја је у томе што Tesla нема могућност приказа слике на екрану. Конфигурација која је коришћена приликом паралелизације LB солвера и тестирања програма које је описано у поглављу 6.4.2 садржи 3 GPU уређаја Tesla C1060, који су повезани на CPU процесор Intel Xeon E5520, радног такта 2.27 GHz. Карактеристике овог GPU уређаја приказане су у табели Ц.2. Докторска дисертација Додатак Ц Тијана Ђукић 143 Табела Ц.2 – Карактеристике NVIDIA Tesla C1060 CUDA верзија уређаја (енг. compute capability) 1.3 Укупна расположива глобална меморија 4096 MB Број мултипроцесора 30 Број CUDA језгара (8 по мултипроцесору) 240 Радни такт GPU уређаја 1.30 GHz Укупан број регистара на располагању једног блока 16384 Максималан број нити у једном блоку 512 Максималне димензије блока 512 x 512 x 64 Максималне димензије мреже блокова 65535 x 65535 x 1 У оквиру овог рада, за извршавање комплексних прорачуна коришћен је CURIE суперкомпјутер, који је у власништву GENCI институције, која се налази у Француској и део је PRACE удружења (енг. PRACE – Partnership for Advanced Computing in Europe). Током ангажовања на пројекту Центра за Биоинжењеринг, омогућен је приступ овом компјутерском кластеру и коришћење њихових ресурса у оквиру Preparatory Access позива за пројекте PRACE организације. CURIE кластер има више сегмената, а у оквиру овог рада и поменутог пројекта коришћени су тзв. хибридни чворови (енг. hybrid nodes) кластера. Овај део кластера садржи 144 хибридна чвора, а сваки чвор има по 2 CPU процесора Intel Westmere X5650, радног такта 2.66 GHz и 2 GPU уређаја Tesla M2090. Укупно у целом кластеру има 288 GPU уређаја, а карактеристике ових GPU уређаја су дате у табели Ц.3. Табела Ц.3 – Карактеристике NVIDIA Tesla M2090 CUDA верзија уређаја (енг. compute capability) 2.0 Укупна расположива глобална меморија 5375 MB Број мултипроцесора 16 Број CUDA језгара (8 по мултипроцесору) 512 Радни такт GPU уређаја 1.30 GHz Укупан број регистара на располагању једног блока 32768 Максималан број нити у једном блоку 1024 Максималне димензије блока 1024 x 1024 x 64 Максималне димензије мреже блокова 65535 x 65535 x 65535 Тијана Ђукић 144 Литература [1] G. Rakocevic, T. Djukic, N. Filipovic, and V. Milutinovic, Computational Medicine in Data Mining and Modeling. New York, USA: Springer, 2013. [2] M. Kojic, N. Filipovic, B. Stojanovic, and N. Kojic, Computer modeling in bioengineering: Thеoretical Background, Examples and Software. Chichester, England: John Wiley and Sons, 2008. [3] E.A. Evans and R. Skalak, Mechanics and thermodynamics of biomembranes. Boca Raton FL, USA: CRC Press, 1980. [4] R. Skalak and S. Chien, Handbook of Bioengineering. New York, USA: McGraw- Hill, 1987. [5] R. E. Waugh and R. M. Hochmuth, "Mechanics and deformability of hematocytes," in Biomedical Engineering Fundamentals. Boca Raton, Florida, USA: CRC Press, 2006. [6] K.S. Chang and W.L. Olbricht, "Experimental studies of the deformation and breakup of a synthetic capsule in steady and unsteady simple shear flow," J. Fluid Mech., vol. 250, pp. 609–633, 1993. [7] A. Walter, H. Rehage, and H. Leonhard, "Shear-induced deformations of polyamide microcapsules," Colloid Polym. Sci., vol. 278, no. 2, pp. 169–175, 2000. [8] P. Gaehtgens, C. Duhrssen, and K.H. Albrecht, "Motions, Deformation and Interaction of Blood Cells and Plasma During Flow Through Narrow Capillary Tubes," Blood cells, vol. 6, pp. 799-812, 1980. [9] A. R. Pries and T. W. Secomb, "Blood Flow in Microvascular Networks," Comprehensive Physiology, pp. 3–36, 2011. [10] R. Skalak, "Rheology of red blood cell membrane," Microcirculation, vol. I, pp. 53- 70, 1976. [11] R.M. Hochmuth and R.E. Waugh, "Erythrocyte membrane elasticity and viscosity," Annu. Rev. Physiol., vol. 49, pp. 209-219, 1987. [12] J. Li, M. Dao, C.T. Lim, and S. Suresh, "Spectrin-level modeling of the cytoskeleton and optical tweezers stretching of the erythrocyte," Biophys J., vol. 88, no. 5, pp. Докторска дисертација Литература Тијана Ђукић 145 3707-19, 2005. [13] D. Kuzman, S. Svetina, R.E. Waugh, and B. Zeks, "Elastic properties of the red blood cell membrane that determine echinocyte deformability," Eur Biophys J., vol. 33, no. 1, pp. 1-15, 2004. [14] R. Mukhopadhyay, H. W. G. Lim, and M. Wortis, "Echinocyte shapes: bending, stretching, and shear determine spicule shape and spacing," Biophys J., vol. 82, no. 4, pp. 1756-72, 2002. [15] C. Pozrikidis, "Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow," J. Fluid Mech., vol. 297, pp. 123–152, 1995. [16] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky, "Fluid vesicles in shear flow," Phys. Rev. Lett., vol. 77, no. 17, p. 3685, 1996. [17] E. Lac, D. Barthès-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., vol. 516, pp. 303–334, 2004. [18] 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., vol. 361, pp. 117–143, 1998. [19] A. Diaz, N. Pelekasis, and D. Barthes-Biesel, "Transient response of a capsule subjected to varying flow conditions: effect of internal fluid viscosity and membrane elasticity," Phys. Fluids, vol. 12, no. 5, pp. 948–957, 2000. [20] C.D. Eggleton and A.S. Popel, "Large deformation of red blood cell ghosts in a simple shear flow," Phys. Fluids, vol. 10, no. 8, pp. 1834–1845, 1998. [21] Y. Sui, Y.T. Chew, P. Roy, and H.T. Low, "A hybrid method to study flow-induced deformation of three-dimensional capsules," J. Comput. Phys., vol. 227, no. 12, pp. 6351–6371, 2008. [22] K. Boryczko, W. Dzwinel, and D. A. Yuen, "Dynamical clustering of red blood cells in capillary vessels," J. Mol. Model., vol. 9, pp. 16–33, 2003. [23] S. Koshizuka, H. Tamako, and Y. Oka, "A particle method for incompressible viscous flow with fluid fragmentation," Comput. Fluid. Dyn. J., vol. 4, pp. 29–46, 1995. [24] K. Tsubota, S. Wada, and T. Yamaguchi, "Particle method for computer simulation of red blood cell motion in blood flow," Comput. Methods. Prog. Biomed., vol. 83, pp. 139–146, 2006. [25] M. Dupin, I. Halliday, C. Care, L. Alboul, and L. Munn, "Modeling the flow of dence Докторска дисертација Литература Тијана Ђукић 146 suspensions of deformable particles in three dimensions," Phys Rev E Stat Nonlin Soft Matter Phys, vol. 6, no. 75, p. 066707, 2007. [26] Y. Liu, L. Zhang, X. Wang, and W. K. Liu, "Coupling of Navier–Stokes equations with protein molecular dynamics and its application to hemodynamics," International journal for numerical methods in fluids, vol. 46, pp. 1237–1252, 2004. [27] S.K. Doddi and P. Bagchi, "Three-dimensional computational modeling of multiple deformable cells flowing in microvessels," Phys. Rev. E, vol. 79, no. 4, pp. 046318– 046331, 2009. [28] X., He and L.-S. Luo, "Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation," Phys. Rev. E, vol. 56, no. 6, pp. 6811– 6817 , 1997. [29] G. McNamara and G. Zanetti, "Use of the Boltzmann equation to simulate lattice-gas automata," Phys. Rev. Let., vol. 61, no. 20, pp. 2332-2335, 1988. [30] C.Y. Lim, C. Shu, X. D. Niu, and Y. T. Chew, "Application of the lattice Boltzmann Method to simulate microchannel flows," Physics of Fluids, vol. 14, no. 7, pp. 2299– 2308., 2002. [31] S. Chen and Z. Tian, "Simulation of microchannel flow by lattice Boltzmann method," Physica A, vol. 388, no. 23, pp. 4803-4810, 2009. [32] N. Filipovic et al., "Numerical and experimental LDL transport through arterial wall," Microfluidics and Nanofluidics, vol. 16, no. 3, pp. 455-464, 2014. [33] N. Filipović, V. Isailović, T. Đukić, M. Ferrari, and M. Kojić, "Multi-scale modeling of circular and elliptical particles in laminar shear flow," IEEE Trans Biomed Eng, vol. 59, no. 1, pp. 50-53, 2012. [34] M. Kojić, R. Slavković, M. Živković, and N. Grujović, Metod konačnih elemenata 1. Kragujevac: Mašinski fakultet Kragujevac, 1998. [35] N. Filipović, Numeričko rešavanje spregnutih problema deformabilnog tela i strujanja fluida. Mašinski fakultet Kragujevac: PhD dissertation, 1999. [36] T. Krüger, F. Varnik, and D. Raabe, "Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method," Computers and Mathematics with Applications, vol. 61, pp. 3485–3505, 2011. [37] R. Skalak, A. Tozeren, R.P. Zarda, and S. Chien, "Strain energy function of red blood cell membranes," Biophys. J., vol. 3, no. 13, pp. 245-264, 1973. Докторска дисертација Литература Тијана Ђукић 147 [38] C. S. Peskin, "Numerical analysis of blood flow in the heart," Journal of Computational Physics, vol. 25, no. 3, pp. 220-252, 1977. [39] Z. Feng and E. Michaelides, "The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problem," Journal of Computational Physics, vol. 195, no. 2, pp. 602-628, 2004. [40] J. Wu and C. Shu, "Particulate flow simulation via a boundary condition-enforced immersed boundary-lattice Boltzmann scheme," Commun. Comput. Phys., vol. 7, no. 4, pp. 793-812, 2010. [41] T. Đukić, Modeliranje solid-fluid interakcije primenom LB metode, Master rad. Kragujevac, 2012. [42] T. Murayama, M. Yoshino, and T. Hirata, "Three-Dimensional Lattice Boltzmann Simulation of Two-Phase Flow Containing a Deformable Body with a Viscoelastic Membrane," Commun. Comput. Phys., vol. 9, no. 5, pp. 1397-1413, 2011. [43] NVIDIA. (2011, May) CUDA Programming guide, Version 4.0. [Online]. http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_ Programming_Guide.pdf [44] C. Sun, C. Migliorini, and L.L. Munn, "Red blood cells initiate leukocyte rolling in postcapillary expansions: A lattice Boltzmann analysis," Biophysical Journal, vol. 85, no. 1, pp. 208-222, 2003. [45] S. Wolfram, Cellular Automaton Fluids 1: Basic Theory.: J. Stat. Phys., 3/4:471–526, 1986. [46] P. L. Bhatnagar, E. P. Gross, and M. Krook, "A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems," Phys. Rev. E, vol. 77, no. 5, pp. 511-525, 1954. [47] A. H. Carter, Classical and Statistical Thermodynamics. New Jersey: Prentice–Hall, Inc., 2001. [48] E. H. Kennard, Kinetic Theory of Gases. New York: McGraw-Hill, 1938. [49] H. Grad, "Note on the N-dimensional Hermite polynomials," Communications on Pure and Applied Mathematics, vol. 2, no. 4, pp. 325–330, 1949. [50] H. Grad, "On the kinetic theory of rarefied gases," Communications on Pure and Applied Mathematics, vol. 2, no. 4, pp. 331–407, 1949. [51] S. Chapman and T. G. Cowling, The mathematical theory of nonuniform gases. Cambridge: Cambridge University Press, 1970. Докторска дисертација Литература Тијана Ђукић 148 [52] К. Huang, Statistical Mechanics. New York: J. Wiley, 1987. [53] O. P. Malaspinas, Lattice Boltzmann Method for the Simulation of Viscoelastic Fluid Flows. Switzerland: PhD dissertation, 2009. [54] T. Đukić, Lattice Boltzmann metod, Završni rad. Kragujevac: Mašinski fakultet, 2010. [55] C. Peng, The Lattice Boltzmann Method for Fluid Dynamics: Theory and Applications. Switzerland: MSc thesis, 2009. [56] V. I. Krylov, Approximate Calculation of Integrals. New York: Macmillan, 1962. [57] P.J. Davis and P. Rabinowitz, Mеthods of Numerical Integration. New York, 1984. [58] X. He, X. Shan, and G. D. Doolean, "Discrete Boltzmann equation model for non- ideal gases," Phys. Rev. E, vol. 57, no. 1, pp. R13 - R16, 1998. [59] P. J. Dellar, "Bulk and shear viscosities in lattice Boltzmann equations," Phys. Rev. E, vol. 64, no. 3, p. 031203, 2001. [60] A. J. Wagner, A Practical Introduction to the Lattice Boltzmann Method. North Dakota State University: PhD dissertation, 2008. [61] S. Chapman and T. G. Cowling, The Mathematical Theory of Non-Uniform Gases.: Cambridge University Press, 1995. [62] T. Shiga, N. Maeda, and K. Kon, "Erythrocyte rheology," Crit. Rev. Oncol. Hematol., vol. 10, pp. 9–48, 1990. [63] N. Maeda and T. Shiga, "Red cell aggregation, due to interactions with plasma proteins," J.Blood. Rheol., vol. 7, pp. 3–12, 1993. [64] D. Barthès-Biesel, A. Diaz, and E. Dhenin, "Effect of constitutive laws for two- dimensional membranes on flow-induced capsule deformation," J. Fluid Mech., vol. 460, pp. 211–222, 2002. [65] M. Živković, Nelinearna analiza konstrukcija. Kragujevac, Srbija: Mašinski fakultet, Kragujevac, 2006. [66] M. Živkovic, Mehanika kontinuuma za analizu metodom konačnih elemenata. Kragujevac, Srbija: WUS Austria, 2011. [67] K. J. Bathe, Finite Element Procedures. Englewood Cliffs, New Jersey, USA: Prentice-Hall, Inc., 1996. [68] D. Mateescu, A. Mekanik, and M.P. Paidoussis, "Analysis of 2-D and 3-D unsteady annular flows with oscillating," Journal of Fluids and Structures, vol. 10, pp. 57–77, Докторска дисертација Литература Тијана Ђукић 149 1996. [69] R. Li, T. Tang, and P.W. Zhang, "A moving mesh finite element algorithm for singular problems in two and three space," Journal of Computational Physics, vol. 177, pp. 365–393, 2002. [70] H.H. Hu, N.A. Patankar, and M.Y. Zhu, "Direct numerical simulations of fluid–solid systems using arbitrary Lagrangian–," Journal of Computational Physics, vol. 169, pp. 427–462, 2001. [71] Z. Guo, C. Zheng, and B. Shi, "Discrete lattice effects on the forcing term in the lattice Boltzmann method," Phys. Rev. E, vol. 65, no. 4, p. 046308, 2002. [72] X. Shan and H. Chen, "Lattice Boltzmann model for simulating flows with multiple phases and components," Phys. Rev. E, vol. 47, no. 3, pp. 1815–1819, 1993. [73] A. Kuzmin, Multiphase simulations with lattice Boltzmann scheme, PhD dissertation ed. Calgary, Canada: University of Calgary, 2009. [74] M.M. Dupin, I. Halliday, and C.M. Care, "Multi-component lattice Boltzmann equation for mesoscale blood flow," J. Phys. A: Math. Gen., vol. 36, pp. 8517–8534, 2003. [75] R. Zhang and H. Chen, "Lattice Boltzmann method for simulations of liquid-vapor thermal flows," Phys. Rev. E, vol. 67, pp. 1–6, 2003. [76] Z. Yu, O. Hemminger, and L.-S. Fan, "Experiment and lattice Boltzmann simulation of two-phase gas-liquid flows in microchannels," Chem. Eng. Sci., vol. 62, pp. 7172– 7183, 2007. [77] G. Tryggvason et al., "A front-tracking method for the computations of multiphase flow," J. Comput. Phys., vol. 169, no. 2, pp. 708–759, 2001. [78] M. Francois, E. Uzgoren, J. Jackson, and W. Shyy, "Multigrid computations with the immersed boundary technique for multiphase flows," Int. J. Numer. Methods Heat Fluid Flow, vol. 14, no. 1, pp. 98–115, 2003. [79] NVIDIA. (2011, May) CUDA C Best practices guide, Version 4.0. [Online]. http://developer.download.nvidia.com/compute/DevZone/docs/html/C/doc/CUDA_C_ Best_Practices_Guide.pdf [80] B. Chapman, G. Jost, and R. van der Pas, Using OpenMP, Portable Shared Memory Parallel Programming. Cambidge, Ma: The MIT Press, 2007. [81] M. Snir, S. Otto, S. Huss-Lederman, D. Walker, and J. Dongarra, MPI-The Complete Reference, Volume 1: The MPI Core, 2nd ed. Cambridge, MA, USA: MIT Press, Докторска дисертација Литература Тијана Ђукић 150 1998. [82] K. Perktold, R. O. K. Peter, M. Resch, and G. Langs, "Pulsatile non-Newtonian blood flow in three-dimensional carotid bifurcation models: a numerical study of flow phenomena under different bifurcation angles," Journal of Biomedical Engineering, vol. 13, pp. 507-515, 1991. [83] N. Filipović, Osnovi bioinženjeringa. Kragujevac, Srbija: Fakultet inženjerskih nauka, 2012. [84] D. N. Ku and D. P. Giddens, "Laser Doppler anemometer measurements of pulsatile flow in a model carotid bifurcation," J. Biomechanics, vol. 20, pp. 407-421, 1987. [85] F. S. Nooruddin and G. Turk, "Simplification and Repair of Polygonal Models Using Volumetric Techniques," IEEE Transactions on Visualization and Computer Graphics, vol. 9, no. 2, pp. 191-205, 2003. [86] Patrick Min. (2013, October) Binvox 3D mesh voxelizer. [Online]. http://www.cs.princeton.edu/~min/binvox/ [87] P. Bagchi, P.C. Johnson, and A.S. Popel, "Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow," Journal of biomechanical engineering, vol. 127, no. 7, pp. 1070-1080, 2005. [88] M. Dao, C.T. Lim, and S. Suresh, "Mechanics of human red blood cell deformed by optical tweezers," J Mech Phys Solids, vol. 51, pp. 2259-2280, 2003. [89] E.A. Evans and Y.C. Fung, "Improved measurements of the erythrocyte geometry," Microvasc Res., vol. 4, no. 4, pp. 335-47, 1972. [90] T. Krüger, M. Gross, D. Raabe, and F. Varnik, "Crossover from tumbling to tank- treading-like motion in dense simulated suspensions of red blood cells," Soft Matter, vol. 9, no. 37, 2013. [91] Y.-Q. Xu and F.-B. Tian, "An efficient red blood cell model in the frame of IB-LBM and its application," International Journal of Biomathematics, vol. 6, no. 1, pp. 1250061-1250083, 2013. [92] S. Haeberle and R. Zengerle, "Microfluidic platforms for lab-on-a-chip applications," Lab on a Chip, vol. 7, pp. 1094–1110, 2007. [93] J. Berthier and P. Silberzan, Microfluidics for Biotechnology. Boston, USA: Artech House, 2006. [94] S. Nagrath et al., "Isolation of rare circulating tumour cells in cancer patients by microchip technology," Nature, vol. 450, pp. 1235-1239, 2007. Докторска дисертација Литература Тијана Ђукић 151 [95] S. Zheng et al., "Membrane microfilter device for selective capture, electrolysis and genomic analysis of human circulating tumor cells," J. Chromatogr. A, vol. 1162, pp. 154–161, 2007. [96] B. Lu, T. Xu, S. Zheng, A. Goldkorn, and T. Yu-Chong, "Parylene membrane slot filter for the capture, analysis and culture of viable circulating tumor cells," in Micro Electro Mechanical Systems (MEMS), Wanchai, Hong Kong , 2010, pp. 935-938. [97] J. Guck et al., "The optical stretcher: a novel laser tool to micromanipulate cells," Biophysical Journal, vol. 81, no. 2, pp. 767-784, 2001. [98] J. Guck et al., "Optical Deformability as an Inherent Cell Marker for Testing Malignant Transformation and Metastatic Competence," Biophysical Journal, vol. 88, no. 5, pp. 3689-3698, 2005. [99] M. Gusenbauer et al., "A tunable cancer cell filter using magnetic beads: cellular and fluid dynamic simulations," arXiv:1110.0995v1 [physics.flu-dyn], 2011. [100] T. Secomb, B. Styp-Rekowska, and A. R. Pries, "Two-Dimensional Simulation of Red Blood Cell Deformation and Lateral Migration in Microvessels," Annals of Biomedical Engineering, vol. 35, no. 5, pp. 755–765, 2007. [101] R. E. Klabunde, Cardiovascular Physiology Concepts. Philadelphia, Pennsylvania, USA: Lippincott Williams & Wilkins, 2011. [102] S. Isogai, M. Horiguchi, and B. M. Weinstein, "The Vascular Anatomy of the Developing Zebrafish: An Atlas of Embryonic and Early Larval Development," Developmental Biology, vol. 230, pp. 278–301, 2001. [103] D. Y. Stainier, B. M. Weinstein, H. W., 3rd Detrich, L. I. Zon, and M. C. Fishman, "Cloche, an early acting zebrafish gene, is required by both the endothelial and hematopoietic lineages," Development, vol. 121, pp. 3141–3150, 1995. [104] B. M. Weinstein, D. L. Stemple, W. D. Driever, and M. C. Fishman, "Gridlock, a localized heritable vascular patterning defect in the zebrafish," Nat. Med., vol. 11, pp. 1143–1147, 1995. [105] A. Swaen and A. Brachet, "Etude sur les premières phases du diveloppement des organes dérivés du mésoblaste chez les poissons téléostéens," Arch. Biol., vol. 16, pp. 173–311, 1899. [106] M. S. Strawinski, "The development of the liver vessels of the sea-trout - Salmo trutta," Bull. Acad. Sci. Cracov, vol. 2, pp. 435–446, 1949. [107] J.-M. Vernier, "Table chronologique du développement embryonnaire de la truite arc- en-ciel Salmo gairdneri Rich 1863," Ann. Embryol. Morphol., vol. 2, pp. 495–520, Докторска дисертација Литература Тијана Ђукић 152 1969. [108] N. D. Lawson and B. M. Weinstein, "In vivo imaging of embryonic vascular development using transgenic zebrafish," Dev. Biol., vol. 248, pp. 307-318, 2002. [109] T. Motoike et al., "Universal GFP reporter for the study of vascular development," Genesis, vol. 28, pp. 75-81, 2000. [110] Jonathan Knight. (2014, April) Facts and information about Zebrafish. [Online]. http://uoneuro.uoregon.edu/k12/zfk12.html [111] E. Ellertsdóttir et al., "Vascular morphogenesis in the zebrafish embryo," Developmental Biology, vol. 341, pp. 56–65, 2010. [112] N. Hu, H. J. Yost, and E. B. Clark, "Cardiac Morphology and Blood Pressure in the Adult Zebrafish," The anatomical record, vol. 264, pp. 1–12, 2001. [113] S. C. Watkins et al., "High Resolution Imaging of Vascular Function in Zebrafish," PLoS ONE, vol. 7, no. 8, p. e44018, 2012. [114] C.C. Huang, C.W. Huang, Y.S Cheng, and J. Yu, "Histamine metabolism influences blood vessel branching in zebrafish reg6 mutants," BMC Dev Biol., vol. 8, p. 31, 2008.