留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

传输线方程的高精度龙格−库塔数值求解方法

王旭桐 周辉 马良 程引会 李进玺 刘逸飞 赵墨 郭景海 王文兵

引用本文:
Citation:

传输线方程的高精度龙格−库塔数值求解方法

    作者简介: 王旭桐(1996—),男,硕士,研究实习员,从事电磁场数值仿真研究;wangxutong@nint.ac.cn.
  • 基金项目: 强脉冲辐射环境模拟与效应国家重点实验室基金(SKLIPR1702)
  • 中图分类号: O441.4

High-precision Runge-Kutta method for transmission line equation

  • CLC number: O441.4

图(8)
计量
  • 文章访问数:  139
  • HTML全文浏览量:  73
  • PDF下载量:  3
出版历程
  • 收稿日期:  2019-10-14
  • 录用日期:  2019-11-25
  • 网络出版日期:  2019-12-26

传输线方程的高精度龙格−库塔数值求解方法

    作者简介: 王旭桐(1996—),男,硕士,研究实习员,从事电磁场数值仿真研究;wangxutong@nint.ac.cn
  • 西北核技术研究院 强脉冲辐射环境模拟与效应国家重点实验室,西安 710024

摘要: 提出了一种求解传输线方程的高精度龙格-库塔(RK)方法。此方法在空间上采取高阶泰勒展开,提高了对空间微分的近似精度,减少了数值色散所带来的误差。与传统的时域有限差分法(FDTD)方法相比,在每波长采样数相同时,RK方法的计算精度更高。同时,根据Taylor模型,对外界平面波激励源进行离散,成功利用RK方法对外部场激励传输线进行求解,扩大了龙格−库塔方法在求解传输线方程时的应用范围。通过编程对平面波辐照下无限大地平面上的单导体与双导体的算例分别应用FDTD方法与RK方法进行了计算,验证了RK方法的正确性。结果表明同等计算条件下RK方法的计算精度更高。

English Abstract

  • 传输线方程是一阶双曲型偏微分方程,求解此类型方程没有现成的差分格式,因此需要把传输线方程转换成一阶拟线性偏微分方程组,然后采用Lax差分格式,偏心格式或者Lax-Wendroff格式对方程组进行求解,此种方法计算量比较大,且转换为一阶拟线性偏微分方程组的推导过程比较复杂[1-2]。本文采用常微分方程数值解法求解传输线的电磁瞬态过程,求解分为两个步骤,一是将传输线方程转换为常微分方程组,二是在初始条件下求解常微分方程组。求解常微分方程组的数值解法有很多,常用的有Euler方法、显示单步法、龙格−库塔(Runge-Kutta,RK)方法、线性多步法等[1-6]。为了提高计算精度与计算效率,本文将RK方法与空间高阶展开结合应用至传输线方程的求解当中,并扩展了此种算法的应用范围。

    • 无源的传输线方程如下

      $ \left\{ \begin{array}{l} - \dfrac{{\partial u}}{{\partial z}} = {R_0}i + {L_0}\dfrac{{\partial i}}{{\partial t}} \\ - \dfrac{{\partial i}}{{\partial z}} = {G_0}u + {C_0}\dfrac{{\partial u}}{{\partial t}} \end{array} \right. $

      (1)

      两端的边界条件如下

      $ {u_0} = - {i_0}{R_s} $

      (2)

      $ {u_l} = {i_l}{R_l} $

      (3)

      式中:$u$为线上电压;${i}$为线上电流;${R_0}$为单位长度电阻;${L_0}$为单位长度电感;${G_0}$为单位长度电导;${C_0}$为单位长度电容。

      将传输线分为N段,每一个节点上都有包含电压值和电流值,整个方程组共有2N个方程,含有$2N + 2$个变量,添加2个边界条件后,只需求解2N个独立变量。

      为了得到空间四阶中心差分近似,采用泰勒级数展开进行处理。

      $ \begin{split} & f(x + \dfrac{1}{2}\Delta x) = f\left( x \right) + \left( { \pm \dfrac{1}{2}\Delta x} \right)\dfrac{\partial }{{\partial x}}f\left( x \right) + \dfrac{1}{{2!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^2}\dfrac{{{\partial ^2}}}{{\partial {x^2}}}f\left( x \right) + \\ & \dfrac{1}{{3!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^3}\dfrac{{{\partial ^3}}}{{\partial {x^3}}}f\left( x \right) + \dfrac{1}{{4!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^4}\dfrac{{{\partial ^4}}}{{\partial {x^4}}}f\left( x \right) + \dfrac{1}{{5!}}{\left( { \pm \dfrac{1}{2}\Delta x} \right)^5}\dfrac{{{\partial ^5}}}{{\partial {x^5}}}f\left( x \right) + \cdots \end{split} $

      (4)

      $ \begin{split} & f(x + \dfrac{3}{2}\Delta x) = f\left( x \right) + \left( { \pm \dfrac{3}{2}\Delta x} \right)\dfrac{\partial }{{\partial x}}f\left( x \right) + \dfrac{1}{{2!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^2}\dfrac{{{\partial ^2}}}{{\partial {x^2}}}f\left( x \right) + \\ & \dfrac{1}{{3!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^3}\dfrac{{{\partial ^3}}}{{\partial {x^3}}}f\left( x \right) + \dfrac{1}{{4!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^4}\dfrac{{{\partial ^4}}}{{\partial {x^4}}}f\left( x \right) + \dfrac{1}{{5!}}{\left( { \pm \dfrac{3}{2}\Delta x} \right)^5}\dfrac{{{\partial ^5}}}{{\partial {x^5}}}f\left( x \right) + \cdots \end{split} $

      (5)

      根据式(4)和式(5)可以得到多项式如下

      $ 27\left[ {f(x + \dfrac{1}{2}\Delta x) - f(x - \dfrac{1}{2}\Delta x)} \right] - \left[ {f(x + \dfrac{3}{2}\Delta x) - f(x - \dfrac{3}{2}\Delta x)} \right] $

      (6)

      采用式(6)代替式(1)中的空间微分,得到空间四阶中心差分后的方程

      $ \left\{ \begin{array}{l} \dfrac{{\partial {u_k}}}{{\partial t}} \approx \dfrac{{\left( {27{i_{k - 1}} + {i_{k + 1}} - 27{i_k} - {i_{k - 2}}} \right)}}{{24\Delta z{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_k} \\ \dfrac{{\partial {i_k}}}{{\partial t}} \approx \dfrac{{\left( {27{u_k} + {u_{k + 2}} - 27{u_{k + 1}} - {u_{k - 1}}} \right)}}{{24\Delta z{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_k} \end{array} \right. $

      (7)

      式中:${u_k} = {u_1},{u_2} \ldots {u_N}$${i_k} = {i_0},{i_1} \ldots {i_{N - 1}}$$N = l/\Delta z$

      由于在空间高阶展开时,求一个点的电压值和电流值需要邻近多个网格点的值,因此在边界以及边界附近采用二阶的展开,对式(1)在$\left( {k + 1} \right)\Delta z$处采取后向欧拉公式进行离散,对式(1)在$k\Delta z$处采取前向欧拉公式进行离散得到

      $ \left\{ \begin{array}{l} \dfrac{{\partial {u_{k + 1}}}}{{\partial t}} \approx \dfrac{{{i_k} - {i_{k + 1}}}}{{\Delta z{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_{k + 1}} \\ \dfrac{{\partial {i_k}}}{{\partial t}} \approx \dfrac{{{u_k} - {u_{k + 1}}}}{{\Delta z{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_k} \end{array} \right. $

      (8)

      将边界条件代入式(8)可以得到边界上的迭代公式

      $ \left\{ \begin{array}{l} \dfrac{{\partial {u_N}}}{{\partial t}} \approx \dfrac{{{i_{N{\rm{ - 1}}}} - {R_l}^{ - 1}{u_N}}}{{\Delta z{C_0}}} - \dfrac{{{G_0}}}{{{C_0}}}{u_N} \\ \dfrac{{\partial {i_0}}}{{\partial t}} \approx \dfrac{{ - {R_{\rm{s}}}{i_0} - {u_N}}}{{\Delta z{L_0}}} - \dfrac{{{R_0}}}{{{L_0}}}{i_0} \end{array} \right. $

      (9)

      式中:${R_{\rm{s}}}$${R_l}$分别为源端和终端电阻。

      将电报方程的空间坐标离散后总共得到2N个相互独立的状态变量,并得到2N阶的线性常态状态方程组:

      $ \frac{{{\rm{d}}{{X}}}}{{{\rm{d}}t}} = {{FX}} + {{s}}(t) $

      (10)

      其中${{s}}\left( t \right)$是离散后的源项,${{X}} = \left[ {array}{l} u \\ i \\ {array} \right]$F为常数矩阵。

      $ {{F}} = \left[\!\!\!\!\!\! {\begin{array}{*{20}{c}} {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\! \ddots \!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24C\Delta {\textit{z}}}}} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{C\Delta {\textit{z}}}}} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{{R_l}C\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{C\Delta {\textit{z}}}}} \\ {\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - {R_s}}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {\dfrac{1}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {\dfrac{{ - 1}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\! \ddots \!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 27}}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{24L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \\ {}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{\dfrac{1}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{\dfrac{{ - 1}}{{L\Delta {\textit{z}}}}}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{}\!\!\!\!\!&\!\!\!\!\!{} \end{array}} \!\!\!\!\!\!\right] $

      (11)

      4级4阶的RK公式如下

      $ \left\{ \begin{array}{l} {y_{n + 1}} = {y_n} + \dfrac{1}{6}h\left( {{k_1} + 2{k_2} + 2{k_3} + {k_4}} \right) \\ {k_1} = f\left( {{x_n},{y_n}} \right) \\ {k_2} = f\left( {{x_n} + \dfrac{1}{2}h,{y_n} + \dfrac{1}{2}h{k_1}} \right) \\ {k_3} = f\left( {{x_n} + \dfrac{1}{2}h,{y_n} + \dfrac{1}{2}h{k_2}} \right) \\ {k_4} = f\left( {{x_n} + h,{y_n} + h{k_3}} \right) \end{array} \right. $

      (12)

      其中,${k_1},{k_2},{k_3},{k_4}$为RK方法的系数,h为步长,$x \in \left[ {{x_0},b} \right]$b可以是有限的,也可以是$ + \infty $y$\left[ {{x_0},b} \right]$R的连续可微函数。

      考虑外界平面波入射,由于RK方法中电流点电压点在同一空间位置,可以得到源项如下

      $ \left\{ \begin{array}{l} - \dfrac{\partial }{{\partial z}}{E_{\rm{T}}}\left( {z,t} \right) + {E_{\rm{L}}}\left( {z,t} \right) = (\dfrac{{{A_{\rm{T}}}}}{{{v_z}}} - \dfrac{{{A_{\rm{L}}}}}{{{v_x}}})\left[ {{\xi _0}\left( {{t^{n + 3/2}} - \dfrac{{k\Delta z}}{{{v_z}}}} \right) - {\xi _0}\left( {{t^{n + 1/2}} - \dfrac{{k\Delta z}}{{{v_z}}}} \right)} \right]/\Delta t \\ - C\dfrac{\partial }{{\partial z}}{E_{\rm{T}}}\left( {z,t} \right) = - C{A_{\rm{T}}}\left[ {{\xi _0}\left( {{t^{n + 3/2}} - \dfrac{{k\Delta z}}{{{v_z}}}} \right) - {\xi _0}\left( {{t^{n + 1/2}} - \dfrac{{k\Delta z}}{{{v_z}}}} \right)} \right]/\Delta t \end{array} \right. $

      (13)

      其中,${E_{\rm{T}}}\left( {z,t} \right)$${E_{\rm{L}}}\left( {z,t} \right)$分别为入射场垂直和平行于传输线导体的分量,${A_{\rm{T}}}$${A_{\rm{L}}}$计算如下

      $ \left\{ \begin{array}{l} {A_{\rm{T}}} = \left( {{e_x}{x_k} + {e_y}{y_k}} \right) \\ {A_{\rm{L}}} = \left( {\dfrac{{{x_k}}}{{{v_x}}} + \dfrac{{{y_k}}}{{{v_y}}}} \right){e_z} \end{array} \right. $

      (14)

      ${e_x}$${e_y}$${e_z}$为方向矢量,$v_x,{v_y},{v_z}$为波沿坐标轴的传播速度

      $ \left\{ \begin{array}{l} {e_x} = \sin {\theta _{\rm{E}}}\sin {\theta _{\rm{p}}} \\ {e_y} = - \sin {\theta _{\rm{E}}}\cos {\theta _{\rm{p}}}\cos {\phi _{\rm{p}}} - \cos {\theta _{\rm{E}}}\sin {\phi _{\rm{p}}} \\ {e_z} = - \sin {\theta _{\rm{E}}}\cos {\theta _{\rm{p}}}\sin {\phi _{\rm{p}}} + \cos {\theta _{\rm{E}}}\cos {\phi _{\rm{p}}} \end{array} \right. $

      (15)

      $ \left\{ \begin{array}{l} v{}_x = - \dfrac{v}{{\cos {\theta _{\rm{p}}}}} \\ {v_y} = - \dfrac{v}{{\sin {\theta _{\rm{p}}}\cos {\phi _{\rm{p}}}}} \\ {v_z} = - \dfrac{v}{{\sin {\theta _{\rm{p}}}\sin {\phi _{\rm{p}}}}} \end{array} \right. $

      (16)

      其中,${\theta _{\rm{p}}},{\varphi _{\rm{p}}},{\theta _{\rm{E}}}$分别为入射波的入射角、方位角、极化角。

    • 根据文献[7]给出的定理,假设初值问题的单步方法的增量函数$\varphi \left( {x,y;h} \right)$关于xyh是连续的,并对于y满足Lipschitz条件,即$\left| {\varphi \left( {x,y;h} \right) - \varphi \left( {x,z;h} \right)} \right| \leqslant L\left| {y - z} \right|$,那么单步法收敛的充分必要条件是单步法是相容的。为了证明显示RK方法的收敛性,只要证明其相容性以及增量函数满足Lipschitz条件就可以了。

      以经典的2级2阶RK方法(中点公式)为例

      $ \varphi \left( {x,y;h} \right) = f\left( {x + \dfrac{1}{2}h,y + \dfrac{1}{2}hf(x,y)} \right) $

      (17)

      $ \mathop {\lim }\limits_{h \to 0} \varphi \left( {x,y;h} \right) = f\left( {x,y} \right) $

      (18)

      上式表明中点公式满足相容性。

      $ \begin{split} & {\rm{ }}\left| {\varphi \left( {x,{y_1};h} \right) - \varphi \left( {x,{y_2};h} \right)} \right| = \left| {f(x + \dfrac{h}{2},{y_1} + \dfrac{1}{2}hf(x,{y_1})) - f(x + \dfrac{h}{2},{y_2} + \dfrac{1}{2}hf(x,{y_2}))} \right| \leqslant \\ & {L_1}\left| {{y_1} - {y_2} + \dfrac{1}{2}h\left[ {f(x,{y_1}) - f(x,{y_2})} \right]} \right| \leqslant {L_1}\left| {{y_1} - {y_2}} \right| + \dfrac{1}{2}h{L_1} \cdot {L_1}\left| {{y_1} - {y_2}} \right| = {L_1}\left( {1 + \dfrac{1}{2}h{L_1}} \right)\left| {{y_1} - {y_2}} \right| \end{split} $

      (19)

      $h \leqslant {h_0}$,取$L = {L_1}\left( {1 + \frac{1}{2}{h_0}{L_1}} \right)$。此时$\varphi \left( {x,y;h} \right)$y满足Lipschitz条件,从而可知中点公式收敛,其他RK方法皆可根据此方法证明收敛。

    • 为了简单起见,一般把数值方法用于试验方程来进行讨论[8]

      $ y'\left( x \right) = \lambda y\left( x \right) $

      (20)

      其中,${\rm{Re}} \left( \lambda \right) < 0$,若$\lambda $为实数,则$\lambda < 0$

      将经典的4级4阶龙格—库塔方法应用至试验方程得到

      $ {y_{n + 1}} = \left[ {1 + \lambda h + \frac{{{{\left( {\lambda h} \right)}^2}}}{2} + \frac{{{{\left( {\lambda h} \right)}^3}}}{6} + \frac{{{{\left( {\lambda h} \right)}^4}}}{{24}}} \right]{y_n} $

      (21)

      $ E\left( {\lambda h} \right) = {1 + \lambda h + \frac{{{{\left( {\lambda h} \right)}^2}}}{2} + \frac{{{{\left( {\lambda h} \right)}^3}}}{6} + \frac{{{{\left( {\lambda h} \right)}^4}}}{{24}}} $

      (22)

      如果$\left| {E\left( {\lambda h} \right)} \right| < 1$,那么$\lambda h$的区域为经典的4级4阶龙格库塔方法的绝对稳定性区域,见图1。其中横轴为实数轴,竖轴为虚数轴。

      图  1  龙格-库塔方法的稳定性区域

      Figure 1.  Stability region of the Runge-Kutta method

    • RK数值解法需要存储一个系数矩阵F,以及4个中间变量${k_1},{k_2},{k_3},{k_4}$。为了尽量减小计算内存,采用稀疏矩阵的存储方式对F进行存储,每得到一个中间变量就把它放入累加器中,可以得到

      $ {\rm{R}}{{\rm{K}}_4} - {\rm{H}}{{\rm{O}}_4}\text{数值解法所需内存} = \left[ {\left( {2N} \right) \times 12 + 2N + 2N} \right] \times 4 $

      (23)
    • 采用经典算例进行计算[8-10],一根半径${r_{\rm{w}}} = 0.254\;{\rm{ mm}}$、长$L = 1\;{\rm{ m}}$的导线位于无限大接地平面以上,距地高度$h = 2\;{\rm{ cm}}$。终端的阻性负载分别为${R_{\rm{S}}} = 500\;{\rm{ \Omega }}$${R_{\rm{L}}} = 1\;000\;{\rm{ \Omega }}$。入射的均匀平面波从顶部入射,${\theta _{\rm{E}}} = 0^\circ $${\theta _{\rm{p}}} = 0^\circ $${\phi _{\rm{p}}} = 0^\circ $。入射电场$\xi \left( t \right)$的时域波形为具有不同上升时间的梯形周期脉冲序列,重复频率$1\;{\rm{ MHz}}$,幅值$1\;{\rm{V/m}}$,占空比50%。脉冲上升时间设置为$10\;{\rm{ ns}}$,如图2所示。

      图  2  传输线模型

      Figure 2.  Transmission line model

      计算结果如图3所示,可以看出当空间步长为$0.02\;{\rm{m}}$时,RK数值方法的计算结果与FDTD方法的计算结果在此算例上相同,验证了本文提出的RK方法的正确性。当空间步长为$0.02\;{\rm{m}}$时,此时RK数值方法依旧稳定,但是由于步长过大导致了计算结果误差较大。

      图  3  传输线左端电压

      Figure 3.  Voltage at left end of twin conductor transmission line

    • 针对多导体平行均匀传输线,设置以下算例,如图4(a)图4(b)所示。导线长2 m,导线半径为0.1905 mm,绝缘层厚度${r_{\rm{w}}} = 0.254\;{\rm{mm}}$,导线中心之间的距离为1.3 mm。均匀平面波沿着z方向入射,电场的极化方向沿x方向。单位长度的参数矩阵$L = \left[\! {array}{l} 0.748\;5\quad {\rm{0}}{\rm{.240\;8}} \\ {\rm{0}}{\rm{.240\;8}}\quad 0{\rm{.748\;5}} \\ {array} \!\right]{\rm{ \mu H/m}}$$C = \left[\! {array}{l} 24.982\quad {\rm{ - 6}}{\rm{.266}} \\ {\rm{ - 6}}{\rm{.266}}\quad 24{\rm{.982}} \\ {array} \!\right]{\rm{ pF/m}}$。传输线终端接$500\;{\rm{\Omega }}$电阻,${R_{\rm{S}}} = {R_{\rm{L}}} = \left[\! {array}{l} {\rm{500}} \quad 0 \\ {\rm{0 }} \quad 500 \\ {array} \!\right]{\rm{ }}\Omega $,入射电场$\xi \left( t \right)$的时域波形如图4(c)所示。入射波入射角度分别为${\theta _{\rm{E}}} = {90^ \circ },{\theta _{\rm{p}}} = {90^ \circ },$${\phi _{\rm{p}}} = - {90^ \circ }$

      图  4  多导体带状电缆模型

      Figure 4.  Model of multiconductor transmission lines

      计算结果如图5所示,结果验证了改进RK方法针对场激励多导体传输线计算的正确性。

      图  5  传输线左端电压

      Figure 5.  Voltage at left end of multiconductor transmission line

    • 针对集中激励源激励的传输线进行计算,采用图6(a)所示的传输线模型,采用图6(b)所示的电压源,线长$0.8\;{\rm{ m}}$${L_0} = 309\;{\rm{ nH/m}}$${C_0} = 144\;{\rm{ pF/m}}$,${R_{\rm{s}}} = 50\;{\rm{ \Omega }},{R_{\rm{L}}} = 50\;{\rm{ \Omega }}$

      图  6  无损传输线模型

      Figure 6.  Lossless transmission line model

      根据传输线基本理论

      $ {Z_0} = \sqrt {{L_0}/{C_0}} $

      (24)

      $ {u_{{R_{\rm{s}}}}}\left( t \right) = \frac{{{Z_0}}}{{{R_{\rm{s}}} + {Z_0}}}{V_{\rm{s}}}\left( t \right) $

      (25)

      $ \beta = \frac{{{u_{\rm{L}}}}}{{{u_{{R_{\rm{s}}}}}}} = \frac{{2{R_{\rm{L}}}}}{{{R_{\rm{L}}} + {Z_0}}} $

      (26)

      其中,${Z_0}$为传输线的特征阻抗,${u_{{R_{\rm{s}}}}}\left( t \right)$为激励源加载至传输线上的电压,$\;\beta $为透射系数,${u_{\rm{L}}}$为负载电压,对式(18)进行求解可以得到负载电压(计算过程中保留有效小数点后8位)

      $ {u_{\rm{L}}}\left( t \right) = \frac{{{Z_0}}}{{{R_{\rm{s}}} + {Z_0}}}{V_{\rm{s}}}\left( t \right)\frac{{2{R_{\rm{L}}}}}{{{R_{\rm{L}}} + {Z_0}}} $

      (27)

      根据式(20)可以画出负载电压波形图,如图7所示。

      图  7  负载电压波形

      Figure 7.  Waveform of load

      FDTD方法选取$\Delta {\textit{z}} = 0.000\;8\;{\rm{ m}}$$\Delta t = 5 \times {10^{ - 12}}\;{\rm{ s}}$,RK数值方法选取$\Delta {\textit{z}} = 0.005\;{\rm{ m}}$$\Delta t = 1 \times {10^{ - 11}}\;{\rm{ s}}$,计算结果如图8所示。

      图  8  终端电压响应

      Figure 8.  Terminal voltage response

      RK方法的计算时间为1.1 s,FDTD的计算时间为4.8 s。从图8可以看出,RK方法与FDTD方法在波形转换的前沿处都存在振荡现象,其中RK方法的振荡现象更为严重,这是由于RK方法的空间步长比FDTD方法的空间步长更大,且RK方法对于激励源的光滑性有更高的要求。从波形平稳处截取10 ns点的波形与解析解进行对比,RK方法与解析解的误差为0.000 449 79,FDTD方法与解析解的误差为0.016 700 00。

      上述算例的计算结果表明,当激励源的函数足够光滑,在同等的计算时间下,RK方法的计算精度更高。

    • 为了改善传输线方程的求解精度,本文将RK方法引入至外界场激励下传输线方程的求解。提高了空间展开阶数,扩展了RK方法应用范围,并对其收敛特性与绝对稳定性做了分析。理论分析表明,激励源的函数足够光滑时,在相同的计算条件下,此方法的计算精度高于传统FDTD方法,通过数值算例验证,同等精度下,改进RK方法消耗的计算时间少于传统FDTD方法。

参考文献 (10)

目录

    /

    返回文章
    返回