基于贝塞尔函数的森林函数化简方法

基于贝塞尔函数的森林函数化简方法

0快速衰减中点电流源的特征分析对能源系统接地网络特征的分析是电气系统兼容性（ec）研究的重要课题。近年来，高压直接供电系统（hvdc）中的电流对交流系统的影响引起了极大的关注。分层土壤中点电流源产生的电流场计算公式是进行地网仿真计算的基本公式之一。文推导了4层分层土壤中点电流源的格林函数公式,但这些公式很复杂。且因被积函数中含有贝塞尔函数,数值积分难以收敛。本文对文导出的格林函数公式进行了化简,得到了一个较简明的表达式。然后分析了公式中的被积函数,将其中衰减较慢的成分提取出来,使用解析法计算其积分;对剩余的快速衰减部分采用数值求积,以贝塞尔函数的零点划分积分区间,逐段积分,并使用外推技术加速积分的收敛。最后使用数值方法计算了点电流源在分层土壤中引起的电流分布,以便建立一个直观的地电流流动图像。1不同电流源的地表电位4层水平分层土壤模型见图1,典型的参数见表1。ρk、hk分别是各层土壤的电阻率和深度,k(k=1,2,3,4)为土壤层数。建立坐标系见图1,r、z分别为场点位置;z1是点电流源所在的位置;设入地电流为I。文导出,位于第1层土壤中的点电流源在该层土壤中产生的电位U=ρ1Ι4π(1√r2+(z-z1)2+∫∞0(θ1e-λz+φ1eλz)J0(λr)dλ)。(1)θ1=e-λz1L4((ρ2+ρ1)(ρ3+ρ2)(ρ4+ρ3)+(ρ2-ρ1)(ρ3+ρ2)(ρ4+ρ3)e-2λh1+(ρ2+ρ1)(ρ3-ρ2)(ρ4+ρ3)e-2λh2+(ρ2+ρ1)(ρ3+ρ2)(ρ4-ρ3)e-2λh3+(ρ2-ρ1)(ρ3-ρ2)(ρ4+ρ3)e-2λ(h2-h1)+(ρ2-ρ1)(ρ3+ρ2)(ρ4-ρ3)e-2λ(h3-h1)+(ρ2+ρ1)(ρ3-ρ2)(ρ4-ρ3)e-2λ(h3-h2)+(ρ2-ρ1)(ρ3-ρ2)(ρ4-ρ3)e-2λ(h3-h2+h1));φ1=1L4((ρ2-ρ1)(ρ3+ρ2)(ρ4+ρ3)e-λ(2h1+z1)+(ρ2-ρ1)(ρ3+ρ2)(ρ4+ρ3)e-λ(2h1-z1)+(ρ2+ρ1)(ρ3-ρ2)(ρ4+ρ3)e-λ(2h2+z1)+(ρ2+ρ1)(ρ3-ρ2)(ρ4+ρ3)e-λ(2h2-z1)+(ρ2+ρ1)(ρ3+ρ2)(ρ4-ρ3)e-λ(2h3+z1)+(ρ2+ρ1)(ρ3+ρ2)(ρ4-ρ3)e-λ(2h3-z1)+(ρ2-ρ1)(ρ3-ρ2)(ρ4-ρ3)e-λ(2h3-2h2+2h1+z1)+(ρ2-ρ1)(ρ3-ρ2)(ρ4-ρ3)e-λ(2h3-2h2+2h1-z1));L4=(ρ2+ρ1)(ρ3+ρ2)(ρ4+ρ3)-(ρ2-ρ1)(ρ3+ρ2)(ρ4+ρ3)e-2λh1-(ρ2+ρ1)(ρ3-ρ2)(ρ4+ρ3)e-2λh2-(ρ2+ρ1)(ρ3+ρ2)(ρ4-ρ3)e-2λh3+(ρ2-ρ1)(ρ3-ρ2)(ρ4+ρ3)e-2λ(h2-h1)+(ρ2-ρ1)(ρ3+ρ2)(ρ4-ρ3)e-2λ(h3-h1)+(ρ2+ρ1)(ρ3-ρ2)(ρ4-ρ3)e-2λ(h3-h2)-(ρ2-ρ1)(ρ3-ρ2)(ρ4-ρ3)e-2λ(h3-h2+h1)。上式形式非常复杂。定义:μk=ρk+1-ρkρk+1+ρk,k=1,2,3。经过化简,U可表示为U=ρ1Ι4π(1/√r2+(z-z1)2+1/√r2+(z+z1)2+∫∞04ΜLcosh(λz1)cosh(λz)J0(λr)dλ)。(2)其中,M=μ1e-2λh1+μ2e-2λh2+μ3e-2λh3+μ1μ2μ3e-2λ(h1-h2+h3);L=1-M+μ1μ2e-2λ(h2-h1)+μ1μ3e-2λ(h3-h1)+μ2μ3e-2λ(h3-h2)。工程上一般关心远离接地点处的地表电位(ESP)。表2给出了I=1kA时按表1参数计算的不同z1电流源引起的地表电位,可见距离接地点>200m时,电流源深度z1对地表电位的影响<0.1%。因z和z1可以换,z的影响可忽略,它意味着在远离接地点处表层土壤中电流是平流的,这可在后文计算结果中得到印证。故在计算远离接地点地表电位时可令式(2)中z1=z=0,得式(3)。至于靠近接地点区域的计算,可认为h2=∞,当作2层土壤处理,即U=ρ1Ι2πr+ρ1Ιπ∫∞0ΜLJ0(λr)dλ。(3)2算法的基本思想式(3)中当r很大时被积函数(M/L)J0(λr)是一个高速振荡函数。对这样的函数求积很困难。观察当λ数值较大时,1-μ1e-2λh1成为L中的主要部分,因此把M/L改写为ΜL=Μ1-μ1e-2λh1+(ΜL-Μ1-μ1e-2λh1)。显然,第1项M/(1-μ1e-2λh1)集中了函数中衰减缓慢的主要部分;第2项M/L-M/(1-μe-2λh1)则随着λ的增大迅速衰减为0,见图2。相应的地表电位计算公式分解为U=ρ1Ι2πr+U0+ˉU。(4)其中,U0=ρ1Ιπ∫∞0Μ1-μ1e-2λh1J0(λr)dλ;ˉU=ρ1Ιπ∫∞0(ΜL-Μ1-μ1e-2λh1)J0(λr)dλ。先讨论U0,注意到分式M/(1-Q)在Q<1时可在形式上展开为Μ/(1-Q)=∞∑n=0QnΜ,故:U0=ρ1Ιπ∞∑n=0∫∞0μn1e-2nλh1ΜJ0(λr)dλ=ρ1Ιπ∞∑n=0∫∞0(μn+11e-2λ(n+1)h1+μn1μ2e-2λ(nh1+h2)+μn1μ3e-2λ(nh1+h3)+μn+11μ2μ3e-2λ(nh1+h1-h2+h3))J0(λr)dλ。由积分公式∫∞0J0(λr)e-λzdλ=1√(r2+z2),得:U0=ρ1Ιπ∞∑n=0(μn+11√r2+4(nh1+h1)2+μn1μ2√r2+4(nh1+h2)2+μn1μ3√r2+4(nh1+h3)2+μn+11μ2μ3√r2+4(nh1+h1-h2+h3)2)。(5)式(5)为一收敛很快的级数。再讨论ˉU,考虑积分S=∫∞0f(x)J0(x)dx。以贝塞尔函数的第偶数个零点x2k(k=1,2,…,并记x0=0)为界划分积分区间,构造积分序列:{ΔSk=∫x2kx2k-2f(x)J0(x)dx,k=0,1,2,⋯}。在每个子区间(x2k-2,x2k)内可以使用高斯—洛巴托求积公式计算积分。J0(x)的第k个零点xk计算公式为xk=β8+1β-1243β3+12092815β5-401743168105β7+867712403β9。式中,β=(8k-2)π。式中最后一项作者添加,主要为改善k较小时xk的计算精度。令Sk=k∑n=1ΔSk,那么S=limk→∞Sk。为加速积分序列的收敛,假定{ΔSk}为等比序列,得到外推式(6);也可以采用多项式外推式(7);或有理式外推(8)。S∞=Sk+ΔSk-1ΔSkΔSk-1-ΔSk；(6)S∞=Sk+x2k-2ΔSkx2k-x2k-2；(7)S∞=Sk+x2k-2ΔS(x2k-x2k-2)-x2kΔS/Sk。(8)外推式(6)～(8)在级数收敛很慢时具有显著的加速效果。但在计算ˉU时,因被积函数中已剔除了衰减很慢的成分,加速效果并不明显。3地表径流分布为得到一个直观的地电流流动图像,使用有限元方法计算了点电流源在分层土壤中产生的电流分布,见图3。计算采用的参数为:h1～h4分别为5、10、15、50m,其它参数同表1。计算结果表明,在表层土壤中,当远离接地点时地电流以水平