生物組織中任意傾斜多柱狀電極的電場計算

生物組織中任意傾斜多柱狀電極的電場計算

針灸和針灸廣泛應(yīng)用于醫(yī)學(xué)。這是因為針灸和針灸可以臨床上的傷口小，使用靈活，可以治療淺表、深部腫瘤和內(nèi)臟器官病變。它是一種有效的補救措施，逐漸成為癌癥和其他疾病綜合治療的有效補充手段。由于單針電極的時間分布是不均勻的[10.12]，治療組織非常大，需要擴大治療范圍，通常需要增加針灸和脊柱針的數(shù)量。同時，為了殺死被破壞的組織的薄細(xì)胞，在殺死組織中需要強烈的橫向分布，以最小的損傷。為了達(dá)到不同的治療效果和目的，對多針傾斜電極的治療要求非常嚴(yán)格。獲得的文獻(xiàn)主要是單元網(wǎng)格阻抗法和有限元法。在現(xiàn)有的限元法中，只使用了兩個針電極的計算示例。多針電極的計算示例是二維平面模型，垂直針灸和莫哈布西翁系統(tǒng)[9.13]。近年來，人們對針灸和莫哈布西翁系統(tǒng)空間分布的定量研究較少。在文獻(xiàn)中，主要有單元網(wǎng)格阻抗法和有限元法。在有限元法中，只使用了兩個針電極。多針電極的計算示例是二維平面模型和垂直針灸電極[9.13]。近年來，人們主要使用積分法原理對多針電極的垂直分布進(jìn)行了計算，但目前還不夠。此外，還有許多計算針灸中產(chǎn)生的電流的計算方法，這些算法主要基于簡化模型而使用電路理論。這是接近計算方法或?qū)ΨQ電離作用的問題。然而，如果這種設(shè)計的電離作用是不規(guī)則的，或多針傾斜電極的電離作用是困難的，它是無法計算復(fù)雜電離作用的?，F(xiàn)在，醫(yī)學(xué)上有許多關(guān)于應(yīng)用的報道，但許多關(guān)于臨床實踐報告的定性描述。在臨床上，為了改變電壓和電極距離來改變治療效果，人們需要對多針傾斜電極的場分布進(jìn)行定量分析和計算。這對醫(yī)學(xué)的繼續(xù)應(yīng)用非常重要。此外，對于多針傾斜電極的動態(tài)分布特征進(jìn)行了進(jìn)一步的探索。在數(shù)值計算中所使用的有限元或有限差分法,其缺點是方程涉及的未知量太多,使計算的速度和精度受到限制.若用于計算傾斜多針狀電極系的場分布,把針狀電極當(dāng)柱狀處理,則網(wǎng)格需劃分的非常小;如果把針狀電極當(dāng)作線狀處理,則會出現(xiàn)奇點;傾斜針狀電極的存在還會出現(xiàn)不規(guī)則立體網(wǎng)格的劃分問題,利用有限元或有限差分法計算傾斜多針狀電極系的場,在邊界的處理上存在困難,因此目前尚無傾斜多針狀電極系所產(chǎn)生場的算例.積分方程法所涉及的未知量只分布在不同介質(zhì)的邊界上,可使計算量減小,從而較大幅度地提高計算速度[23~32].同時積分方程法將電流分布按線狀處理,可以很好地解決用有限元素或有限差分法計算傾斜多針狀電極系邊界處理所遇到的困難.但是用積分方程法計算靜電類問題時,有時離散方程會存在病態(tài)問題.我們利用算例對邊界元積分方程系數(shù)矩陣進(jìn)行了考查,所選算例屬于良態(tài)方程.本文實現(xiàn)了利用積分方程法計算傾斜多針狀電極系的電場分布,并獲得了較快的計算速度.1計算電流分布的性能1.1含多個含電極的x-u-t-ln點電流源I在生物導(dǎo)電組織中任意一點的電勢可以寫為其中r為場點到源點的距離,σ是生物組織的電導(dǎo)率.圖1為任意多針狀電極的不規(guī)則空間分布圖,若有L1,L2,…,Lm,…,Ln個針狀電極(編號),并用L1,L2,…,Lm,…,Ln代表針狀電極的長度.設(shè)由單位長度針狀電極流入生物組織的電流為I1,I2,…,Im,…In,且I1,I2,…,Im,…,In在同一個電極上的分布是不均勻的、隨空間坐標(biāo)而變化,則有其中U=[U1,U2,L,Um,L,Un]T,U1,U2,…,Um,…,Un分別為針狀電極上的電勢,上標(biāo)T表示矩陣的轉(zhuǎn)置,r1,r2,…,rm,…,rn分別為場點到源點的距離.1.2轉(zhuǎn)化單元i,j設(shè)場點和源點都在針狀電極上,將L1,L2,…Lm,…,Ln個電極分別分為M1,M2,…,Mm,…,Mn個網(wǎng)格,采用塊中心網(wǎng)格系統(tǒng),總網(wǎng)格數(shù)為M=M1+M2+L+Mm+L+Mn,則第j個網(wǎng)格在第i個網(wǎng)格所產(chǎn)生的電勢為其中ir是網(wǎng)格編號為i的中心點的矢徑,r′為第j個網(wǎng)格上任意一微分元dl的矢徑.若i和j分別為多針狀電極系中任意兩針狀電極Lm1和Lm2上的任意兩個網(wǎng)格(i和j也可以在同一根針狀電極上,如圖2),設(shè)任意兩針狀電極Lm1和Lm2下端點的坐標(biāo)為(x01,y01,z01)和(x02,y02,z02),第i個網(wǎng)格中心點的坐標(biāo)為(x1,y1,z1),針狀電極Lm1和Lm2的方向余弦為(cosα1,cosβ1,cosγ1)和(cosα2,cosβ2,cosγ2),α11β,γ1和α2,β2,γ2分別為針狀電極Lm1和Lm2與坐標(biāo)軸x,y,z的夾角.在第j個網(wǎng)格中認(rèn)為I(j)為常數(shù),若I(j)系數(shù)用a(i,j)表示,則有其中在這里il和lj表示由針狀電極Lm1和Lm2的下端點到第i和j個網(wǎng)格中心的距離.同樣有則令則(3)式積分可表示為其中l(wèi)j+12和lj-12是第j個網(wǎng)格上下邊界點到電極下端點的距離.若i=j,則x01=x02,y01=y02,z01=z02,α1=α2,β1=β2,γ1=γ2,有此時(3)式的積分為含有奇點的積分,可采取如下處理辦法:在i點附近挖去長度為2?的一小段?→0,則在實際計算中根據(jù)精度要求適當(dāng)選取?的大小.當(dāng)i和j均在同一個電極上時,不論j點在i點的上方還是下方,可以證明矩陣元a(i,j)可用統(tǒng)一的(4)式計算.電極上第i點的電勢為將(6)式寫成矩陣形式為2計算空間勢和電池電壓的值2.1第j個網(wǎng)格設(shè)電極mL下端點的坐標(biāo)為(x0m,y0m,z0m),第j個網(wǎng)格是電極mL上的點,r′為第j個網(wǎng)格上任意微分元dl的位矢,r為空間任意一點P的位矢,則P點的電勢為2.2eyj-mezj模型點電流元在電導(dǎo)率為σ的生物組織中的電場可寫為其中r和r′為場點和源點的位置矢徑.與計算電勢方法類似,第j個網(wǎng)格電流元在生物組織中任意一點所產(chǎn)生場強的x分量為其中αm為向量(r-r′)與x軸的夾角.Eyj,Ezj與Exj類似有總電場的各分量為其中βm和γm為向量(r-r′)與y和z軸的夾角.總電場強度為3計算值的例和計算結(jié)果的驗證3.1實驗結(jié)果和討論在傾斜多針狀電極無解析解和沒有其他數(shù)值計算結(jié)果的條件下,為考察和檢驗計算方法的可行性和計算結(jié)果的可靠性,今用文獻(xiàn)[15,16]在水中的實驗測量結(jié)果對傾斜多針狀電極系計算方法進(jìn)行了檢驗,將傾斜多針狀電極系程序簡化(通過改變程序中的電極系參數(shù))為豎直2針狀電極系.得到了如表1的對比結(jié)果(計算及實驗參數(shù),兩電極長度L1=L2=0.07m,陰極1L電勢1U=0V,電源陽極2L電勢U2=3V,1L平行于2L,1L距2L0.04m,兩電極下端點在垂直于兩電極的同一個平面內(nèi),水的電導(dǎo)率σ=1.2×10-4(??m)-1,實驗時整個裝置放在水中,電勢參考0點為陰極),在表1中UC表示由邊界元方法計算的結(jié)果,Um為由實驗測量所測得的電勢分布.可以看出計算結(jié)果與實驗值符合得非常好,而在離陽極較近的點電勢測量值低于計算值,其相對誤差最大值小于1%,在電極上電勢測量值等于計算值,說明計算方法可靠.3.2傾斜多軸電極系統(tǒng)的場勢值計算3.2.1電極形貌和位置分布利用MATLAB軟件平臺,我們編程實現(xiàn)了任意傾斜多針狀電極系的場勢計算,下面給出了一個7針狀傾斜電極的算例,圖3中A,B,C,D,E,F,G為7針狀傾斜電極下端在xoy坐標(biāo)平面上的位置,呈正六邊形排列,各點之間距離相同均為0.02m,各點的坐標(biāo)(x,y)分別為:A(0.1m,0.135m),B(0.12m,0.135m),C(0.13m,0.117m),D(0.12m,0.1m),E(0.1m,0.1m),F(0.09m,0.117m),G(0.11m,0.117m),各電極與x,y,z軸的夾角α,β,γ分別為LA(70°,70°,29°),LB(50°,50°,66°),LC(60°,60°,45°),LD(70°,70°,29°),LE(70°,70°,29°),LF(60°,60°,45°),LG(90°,90°,0°).在這7針狀電極中只有電極G是豎直放置的,每個電極的電勢均為U=200V,電極長度均為0.05m,每根針狀電極被離散為50個網(wǎng)格,生物組織的電導(dǎo)率σ=0.286(??m)-1.計算結(jié)果見圖4~6.圖4表示在坐標(biāo)平面z=0上的電勢分布,在靠近電極的地方電勢存在大的梯度變化,在電極上利用邊界元法計算的結(jié)果等于電極的電勢設(shè)定值200V由于所有電極的下端都在該平面上,所以電勢出現(xiàn)7個銳鋒.盡管電極下端點在xoy平面(z=0)上是均勻分布的,但因電極是傾斜的,所以其電勢等值線圖4(b)為非對稱分布.圖4(b)明顯顯示出電極所在的位置.圖5表示在平面z=0.04上的電勢分布,該圖只顯示出4根針狀電極的存在,這是因為只有這4根針狀傾斜電極被該平面所切割,電勢等值線圖5(b)顯示電極在該平面的位置間距存在較大差別,在電極區(qū)內(nèi)電勢分布更不均勻.圖6表示在直線1L:y=0.117m,z=0.0m;2L:y=0.117m,z=0.01m;3L:y=0.117m,z=0.05m上的電場強度分布.1L,2L,3L互相平行,直線1L穿過三個針狀電極F、G、C的下端,三條直線均穿過豎直電極G,但2L,3L不與F、C電極相交,圖中的曲線1,2,3分別對應(yīng)直線1L,2L,3L上的電場分布,圖6(a)~(c)為這三條直線上的電場x,y,z分量,圖6(d)為總電場強度的大小,曲線2顯示,在x=0.14m附近直線2L與其中一個電極比較靠近.電場的Ex分量在電極的左右兩側(cè)存在方向的改變,顯示出傾斜多針狀電極電場分布的不均勻性.3.2.2算例3:美國法上的最佳計算方法積分方程法的優(yōu)勢是所涉及的未知量只分布在不同介質(zhì)的邊界上,但計算發(fā)現(xiàn)離散方程有時解的穩(wěn)定性差.對上述電極結(jié)構(gòu)我們建立了以傾斜針狀電極系的電流密度為未知量的離散方程,分別計算了離散方程系數(shù)矩陣的條件數(shù)(直接調(diào)用MATLAB計算條件數(shù)的指令可實現(xiàn)系數(shù)矩陣條件數(shù)的計算),結(jié)果表明隨著針狀電極劃分網(wǎng)格數(shù)目的增加,系數(shù)矩陣的條件數(shù)在增大,但增加的非常緩慢.電極網(wǎng)格總數(shù)由70變到2100的過程中,方程系數(shù)矩陣條件數(shù)由1.709變?yōu)?.405,將方程系數(shù)矩陣條件數(shù)與機器精度(10-16)相乘則可粗略估計計算的誤差,若使用雙精度,對于傾斜多針狀電極系的計算機截斷誤差最大在10-15數(shù)量級范圍內(nèi),方程屬于良態(tài)方程,計算方法的穩(wěn)定性較好.表2還給出了該方法的計算程序運行時間與針狀電極網(wǎng)格數(shù)目的關(guān)系,計算時將空間場點劃分為100×100×100=1000000個網(wǎng)格.在不改變空間網(wǎng)格總數(shù)條件下,研究了程序運行時間與電極系網(wǎng)格劃分個數(shù)的關(guān)系(見表2)(軟件:MATLAB軟件平臺,程序:傾斜多針狀電極程序(自己編寫的),機器:PentiumⅣ650微型計算機).表2顯示,隨針狀電極上網(wǎng)格總數(shù)的增加,程序計算運行時間在非線性增加,但當(dāng)針狀電極網(wǎng)格數(shù)為2100時計算空間100萬個場點的運行時間為116.2s,這在僅有0.05m長的針狀電極上就劃分了300個網(wǎng)格(每個網(wǎng)格的長度不到0.2mm),對這種排列的傾斜7針狀電極系,有限元或有限差分方法是難以精確到這種網(wǎng)格尺度的,在這么短的時間內(nèi)無法實現(xiàn)這么多場點的電勢計算,本算例顯示了邊界元積分方程法在計算速度上的優(yōu)勢在實際應(yīng)用中在0.05m長的針狀電極上劃分20~30個網(wǎng)格足以滿足精度的要求,那么空間100萬個場點電勢的運行時間僅為3~5s,這為射頻治療儀的職能控制提供了快速計算方法,該表還反映出電極網(wǎng)格數(shù)的多少對程序運行時間影響較大,而計算空間場點數(shù)目的多少對程序運行時間影響較小.4結(jié)果的評價和應(yīng)用本文利用邊界元法實現(xiàn)了生物組織中傾斜多針狀電極所產(chǎn)生場的數(shù)值計算,克服了利用有限元素或有限差分法計算生物組織中傾斜多針狀電極系問題所遇到的邊界處理上的困難.并用實驗測量結(jié)果對計算方法進(jìn)行了檢驗,結(jié)果吻合得較好.還給出了7針狀傾斜電極算例,利用算例對方法的計算速度和穩(wěn)定性進(jìn)行了考查,獲得了較快的計算速度和穩(wěn)定的計算結(jié)果,表明邊界元積分方程法對計算傾斜多針狀電極系的電場分布是合適的,較其他方法有計算速度快的明顯優(yōu)勢.這對解決醫(yī)學(xué)實際問題具有重要意義,只要給