Диссертация (1097617), страница 10

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 10)

Плотность падающего потокаэнергии на поперечную поверхность хрящевой ткани имеет пространственноераспределение, соответствующее распределению Гаусса с эффективным радиусом пучкаr0(x), учитывающим расхождение пучка в поперечном направлении по мере прохождениячерез среды в направлении x. y 2 + z 2  κ exp ( −κ x )G ( x, y, z,τ ) = P (τ ) exp  − 2cρ r0 ( x ) (2.2)здесь P(τ) — зависящая от времени мощность лазерного излучения, c — удельнаятеплоемкость, ρ — плотность.На границах раздела сред записывались граничное условие постоянства плотноститеплового потока:λm∂TSm (τ )∂T n (τ )= λn S∂i∂i(2.3)или конвективного теплообмена между окружающей средой и поверхностью хрящаносовой перегородки, покрытой слизистой оболочкой:qSm (τ ) = α (TS (τ ) − TS 0 (τ ) )(2.4)где m и n — индексы, нумерующие соприкасающиеся среды, индекс i — x, y или z, λ— коэффициент теплопроводности, S0 — ближайшие к поверхности раздела точкиокружающей среды, S — точки поверхности раздела.Для решения дифференциального уравнения теплопроводности (2.1) с учетомграничных условий (2.3) и (2.4), удовлетворяющих конкретной геометрии лазерноиндуцированного нагрева, применялся метод конечных разностей (метод сеток), основанныйна замене производных их приближенным значениями, выраженными через разностизначений функций в узлах сетки.С учетом цилиндрической симметрии рассматриваемой задачи уравнениетеплопроводности можно записать в цилиндрической системе координат (ЦСК) в виде:cvmγ m ∂ 2T (r, θ , z, τ ) 1 ∂T (r, θ , z, τ ) 1 ∂ 2T (r, θ , z, τ ) ∂ 2T (r, θ , z, τ ) ∂T (r, θ , z, τ ) + G (r, θ , z, τ )= λm ++ 2+∂τ∂r 2r∂rr∂θ 2∂z 2(2.5)43Так как в рассматриваемом случае ось нагреваемого пространства совпадает с осью z,а начальные и граничные условия и симметрия пространства не зависят от угловойкоординаты, то имеем: ∂ 2T (r, z, τ ) 1 ∂T (r, z, τ ) ∂ 2T (r, z, τ ) ∂T (r, z, τ ) + G (r, z, τ )cvmγ m= λm ++∂τ∂r 2r∂r∂z 2(2.6) м2 здесь m - номер среды, am  — коэффициент тепловой диффузии (по Кельвину) сек или коэффициент температуропроводности при постоянном объеме, выражающийся через кг плотность среды γm  3  , коэффициент теплопроводности λmм  Вт  м град  и удельную Дж λизохорную теплоемкость cv : am = m .γ m cvm кг град Функция Gm (r, z, τ ) в ЦСК, описывает воздействие лазерного излучения, какобъемного источника тепла, и включает в себя коэффициент1cvmγ m, зависящий отпараметров каждой среды.

Плотность падающего потока энергии на поперечнуюповерхность среды m имеет пространственное распределение, соответствующее вr2рассматриваемой задаче распределению Гаусса: exp[−] , где величина r0 ( z ) r0 ( z ) 2эффективный радиус пучка, который, в общем случае, учитывает возможное расхождениепучка по мере прохождения через среды и зависит от координаты z распространения вглубину.Предполагается, что при прохождении через поглощающую и рассеивающую средуослабление интенсивности плотности падающего потока в среде m описывается закономБугера-Ламберта-Бера с эффективным показателем поглощения κ m :r2q(r, z,τ ) = (1 - K )q 0 (τ ) exp[− 2 ] f (t ) exp(−κ m z )r0rmGm (r, z, τ ) =dq(r, z,τ )1r2=(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mdzcvmγ mr0 ( z ) 21(2.7)(2.8)44 ∂ 2T (r, z, τ ) 1 ∂T (r, z, τ ) ∂ 2T (r, z, τ ) ∂T (r, z, τ ) += am ++22∂τ∂rr∂r∂z(2.9)r2+(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 ( z ) 21При рассмотрении данной задачи положим K mr = 0 , то есть пренебрегаем отражениемот границ оптически непрозрачных сред, так как эта величина сравнима с погрешностьюизмерения мощности.

Положим только на первой границе оптоволокно-хрящ отражениесоответствующее 4% затухания мощности [Борн и др., 1973].Мощность излучения Pm [Вт], проходящая через поперечный срез m пространства,равная поверхностному интегралу теплового потока с амплитудой q Sm (τ ) или объемномуинтегралу от всего источника тепла с амплитудой q Vm (τ ) , в зависимости от объемного илиплоскостного представления источника излучения:Pm = ∫ q mS (τ ) exp[−Sy2 + z2y2 + z2V]dydz=q(τ)exp[−]κ m exp(−κ m x) dxdydz∫ mr0 ( x) 2r0 ( x) 2V(2.10)Вычисляя (численно) эти интегралы и приравнивая их измеренной величинемощности, можно получить амплитуды q mS (τ ) и q Vm (τ ) . Рассмотрим амплитуду объемногоисточника тепла. Введя обозначение: qVm (τ ) = PowerLaserm (τ ) , получим: ∂ 2T (r, z, τ ) 1 ∂T (r, z, τ ) ∂ 2T (r, z, τ ) ∂T (r, z, τ ) += am ++∂τ∂r 2r∂r∂z 2r2+(1 - K mr ) PowerLaserm exp[−]κ m exp(−κ m z )cvmγ mr0 ( z ) 2(2.11)12.2.

Численное моделирование лазерного нагрева биологической ткани при внешнеммеханическом воздействии на ее поверхностьС учетом геометрии, дифференциальное уравнение теплопроводности невозможнобыло решить аналитически, и для решения задачи использовался метод конечных разностей(метод сеток), основанный на замене производных их приближенными значениями,выраженными через разности значений функций в узлах сетки. При этом дифференциальныеуравнения заменялись эквивалентными соотношениями в конечных разностях, итемпература на следующем шаге временной сетки выражалась через значения температурына двух предыдущих временных слоях [Бахвалов и др., 1977].45Введем индексацию по осям: ir = 1, NR + 1 , iz = 1, NZ + 1 , и индексацию повремени j = 1, NJ , где черта сверху обозначает изменение индексов по осям в пределах,указанных под чертой с единичным шагом.T (r, z, τ ) = Tir , iz , jj +1jj +1j −1∂T (r, z, τ ) Tir , iz − Tir , iz Tir , iz − Tir , iz==∂τdt2dt∂T (r, z, τ ) Tir +1, iz , j − Tir , iz , j=∂rdr∂T (r, z, τ ) Tir , iz +1, j − Tir , iz , j=∂zdzjjjjj +1j −1j∂ 2T j (r, z, τ ) Tir +1, iz − 2Tir , iz + Tir −1, iz Tir +1, iz − Tir , iz − Tir , iz + Tir −1, iz===∂r 2(dr )2(dr )2==Tirj+1, iz + Tirj−1, iz(dr )2−Tirj,+iz1 + Tirj,−iz1(dr )2Tirj+1, iz − 2Tirj, iz + Tirj−1, iz(dr )2± 2Tjir , iz=Tirj+1, iz − 2Tirj, iz + Tirj−1, iz(dr )2−Tirj,+iy1 , iz − 2Tirj, iz + Tirj,−iz1(dr )2j +1jj −12(dt ) Tir , iz − 2Tir , iz + Tir , iz−,(dr )2(dt )2=jjjjj +1j −1j∂ 2T j (r, z, τ ) Tir , iz +1 − 2Tir , iz + Tir , iz −1 Tir , iz +1 − Tir , iz − Tir , iz + Tir , iz −1===∂z 2(dz )2(dz )2==Tirj, iz +1 + Tirj, iz −1(dz )2−Tirj,+iz1 − Tirj,−iz1(dz )2Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1(dz )2−± 2Tjir , iz=Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1(dz )2(dt )2 Tirj,+iz1 − 2Tirj, iz − Tirj,−iz1 ,(dz )2(dt )2−Tirj,+iz1 − 2Tirj, iz − Tirj,−iz1(dz )2тогда уравнение теплопроводности (2.11) имеет вид:Tirj,+iz1 − Tirj,−iz12dt ∂ 2T (r, z, τ ) 1 ∂T (r, z, τ ) ∂ 2T (r, z, τ ) 1r2r= am +++(1K)q()exp[−]κ m ×τmm c γ∂r 2r∂r∂z 2r0 (r ) 2 vm m Tirj+1, iz − 2Tirj, iz + Tirj−1, iz (dt )2 Tirj,+iz1 − 2Tirj, iz + Tirj,−iz1 1 Tir +1, iz , j − Tir , iz , j× exp(−κ m z ) = am −++222rdr()()()drdrdt+Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1(dz )2j +1jj −12(dt ) Tir , iz − 2Tir , iz − Tir , iz 1r2r−+(1 - K m )q m (τ ) exp[−]κ m exp(−κ m z ) cvmγ mr0 ( r ) 2(dt )2(dz )2(2.12)Произведем перегруппировку:46Tirj,+iz12dt= (dt )2 Tirj,+iz1 − 2Tirj, iz + Tirj,−iz1 (dt )2 Tirj,+iz1 − 2Tirj, iz − Tirj,−iz1 +− am + (dr )22dt(dt )2(dz )2(dt )2Tirj,−iz1 T j − 2Tirj, iz + Tirj−1, iz 1 Tir +1, iz , j − Tir , iz , j Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1 ++ am  ir +1, iz++22rdr()()drdz21r+(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 (r ) 2члены видаTirj,+iz1 − 2Tirj, iz + Tirj,−iz1(dt )2по двум осям, являются второй производной повремени и дают устойчивость сетки.

Это трехслойная схема «ромб» или схема ФранкелаДюфорта [Калиткин, 1978], в которой значения функции на втором временном слоерассчитывались по явной центральной четырехточечной схеме, а значение сеточнойфункции на третьем временном слое рассчитывались по ее значениям на двух предыдущихвременных слоях.Tirj,+iz1Tirj,−iz1 Tirj,+iz1 − 2Tirj, iz + Tirj,−iz1 Tirj,+iz1 − 2Tirj, iz − Tirj,−iz1 +=− am +222dt2dt()()drdz Tirj+1, iz − 2Tirj, iz + Tirj−1, iz 1 Tirj+1, iz − Tirj, iz Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1 ++ am ++22rdr(dr)(dz)21r+(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 (r ) 2Tirj,+iz12dt=Tirj,−iz1 11 +− am Tirj,+iz1 − 2Tirj, iz + Tirj,−iz1 +22 2dt(dr)(dz)() T j − 2Tirj, iz + Tirj−1, iz 1 Tirj+1, iz − Tirj, iz Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1 ++ am  ir +1, iz++22rdr(dr)(dz)21r+(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 (r ) 2Tirj,+iz1j −1 11  j +1 Tir , iz+ am +T=− am − 2Tirj, iz + Tirj,−iz122  ir , iz2dt2dt (dr ) (dz ) ()11 ++22  (dr ) (dz )  Tirj+1, iz − 2Tirj, iz + Tirj−1, iz 1 Tirj+1, iz − Tirj, iz Tirj, iz +1 − 2Tirj, iz + Tirj, iz −1 ++ am ++rdr(dr )2(dz )221r+(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 (r ) 247 1 1 1 1 11  1  1  j = Tirj,−iz1  + am 2Tir , iz −Tirj,+iz1 + am +− am ++22 22 22  (dr ) (dz )   (dr ) (dz )   (dr ) (dz )  2dt 2dt Tirj+1, iz + Tirj−1, iz 1 Tirj+1, iz Tirj, iz +1 + Tirj, iz −1  11 11  j+2Tir , iz + am − am ++++22 r dr(dr )2(dz )2  (dr ) 2r dr (dz ) r2+(1 - K )q m (τ ) exp[−]κ m exp(−κ m z )cvmγ mr0 ( r ) 21rm 1 1 1 11  1  1 j = Tirj,−iz1  − amTirj,+iz1 + am +− am +Tir , iz +22 22 rdr (dr ) (dz )   (dr ) (dz )   2dt 2dtjT j +T jT j +T j 1T1r2+ am  ir +1, iz 2ir −1, iz + ir +1, iz + ir , iz +1 2ir , iz −1  +(1 - K mr )q m (τ ) exp[−]κ m exp(−κ m z )2 c γrdrr(r)()()drdz0vmm 1 11  1  1 j = Tirj,−iz1 1 − 2dtam  − 2dtamTirj,+iz1 1 + 2dtam ++Tir , iz +22 22rdr(dr)(dz)(dr)(dz)jT j +T jT j +T j 1T(1 - K mr )r2+ 2dtam  ir +1, iz 2ir −1, iz + ir +1, iz + ir , iz +1 2ir , iz −1  + 2dtq m (τ ) exp[−]κ m exp(−κ m z )2γrdrcr(r)()()drdzvmm0(2.13)Отсюда можно выразить температуру на следующем шаге временной сетки череззначения температуры на двух предыдущих слоях.Для дальнейшего численного моделирования введем некоторые обозначения:arrm =dtdtdta , arm =am , azz m =a , Coef m = dt (1 − K mr )κ m , тогда (2.13) можно2 m2 mrdr(dr )(dz )записать в виде:Tirj,+iz1 = Tirj,−iz1+(1 − 2(arrm + azzm )) +2arm2arrm((Tirj+1, iz − Tirj, iz ) +Tirj+1, iz + Tirj−1, iz ) +(1 + 2(arrm + azzm )) (1 + 2(arrm + azzm ))(1 + 2(arrm + azzm ))()2azzm2Coef mr2Tirj, iz +1 + Tirj, iz −1 +q m (τ ) exp[−] exp(−κ m z )(1 + 2(arrm + azzm ))(1 + 2(arrm + azzm )) cvmγ mr0 (r ) 2(2.14)Граничные условия:На внешних границах рассматриваемого пространства имеем граничное условие 2-огорода или 3-его, соответствующее заданию плотности теплового потока для каждой точкиповерхности в любой момент времени или оттоку тепла с поверхности в окружающую среду:Граничное условие2-ого рода на стыке двух сред при r = rm48∂T m (rm , z, τ )∂T m+1 (rm , z, τ )= λm+1∂r∂rmm +1T (rm , z, τ ) = T (rm , z, τ )λmАналогичные уравнения записываются для всех остальных внутренних границсоприкосновения, перпендикулярных оси z при z = z n∂T n (r, z m , τ )∂T n+1 (r, z n , τ )= λn+1∂z∂znn +1T (r, z n , τ ) = T (r, z n , τ )λn(2.15)Граничные условия 3-ого рода:qпов = -λgradT = -λen Вт∂T= β(Tсреды-Tпов), где β — коэффициент теплообмена  2 град∂nмТаким образом, для стыков среды с воздухом, например, на стороне биологическойткани, противоположной стороне лазерного облучения, можно записать:∂T (r, z β , τ )∂z=β(T (r, z β , τ ) − Tair ),λT ir , NZ , j − T ir , NZ −1, jdz=β(T (r, z β τ ) − Tair )λ(2.16)2.3.

Характеристики

Список файлов диссертации