1、 南京航空航天大学南京航空航天大学 何小祥何小祥 有限元发展历史有限元发展历史有限元的优点有限元的优点二维节点有限元（有限元处理过程）二维节点有限元（有限元处理过程）二维高阶节点有限元二维高阶节点有限元三维节点有限元方法三维节点有限元方法三维高阶节点有限元方法三维高阶节点有限元方法节点有限元的缺点节点有限元的缺点矢量有限元方法矢量有限元方法二维、三维一阶矢量基函数二维、三维一阶矢量基函数高阶矢量基函数高阶矢量基函数有限元中的现代技术有限元中的现代技术时域有限元方法时域有限元方法有限元方程组求解有限元方程组求解-CG1、有限元发展历史、有限元发展历史2、有限元的优点、有限元的优点 非结构网格-可
2、以模拟任意形状的边界 泛函的思想-利于复杂介质特性模拟有限元处理过程有限元处理过程边值问题等价泛函目标建模，网格剖分参数提取矩阵填充局部矩阵计算FEM方程组求解后处理技术差值基函数只需差值基函数只需要在单元内满足要在单元内满足精度要求精度要求 U代表位函数或者某个场分量12一般计算区域比较复杂，积分不容易实现，采用整体基函数很难满足齐次边界条件，采用了广义变分原理3、二维节点有限元（有限元处理过程）、二维节点有限元（有限元处理过程）kuiujuxyo3Patran Ansys所需参数所需参数-坐标，单元编号，节点或棱边编号坐标，单元编号，节点或棱边编号4222211()MMetcttetuuF
3、 uFk udxdyxy 将波导横截面用M个三角形面元离散，则各面元中的u 用线性插值函数表示为iiietuNu31积分区域变为三角形单元内思想、思想、理解理解能量的概念5()/2iiiiNab xc y3N1N2N物理含义123()/2iiiiNab xc y参数提取优点优点222211()MMetcttetuuF uFk udxdyxy 利用横等式可以获得局部矩阵的解析解2102()eeeeTee TeFEAEkEBE1 reeTiiVNNAdxxeeeeeTrVBNNd V6单元一、节点2.4.1方程组矩阵填充方程组矩阵填充7（1）（2）（3）222211()MMetcttetuuF u
4、Fk udxdyxy 单元二、节点2.4.5单元三、节点3.5.2单元四、节点4.5.6单元节点与相邻单元节点必须吻合，过程中使用到各个单元中节点局部编号与整体编号的关系，列向量等可以如上进行填充最后形成的系数矩阵是一个对称的稀疏矩阵最后形成的系数矩阵是一个对称的稀疏矩阵 UBkUAc2可导出线性代数方程组由此方程解出本征值和对应的本征矢。980,1,2,pFpNu222211()MMetcttetuuF uFk udxdyxy 4、二维高阶节点有限元二维高阶节点有限元1(,)(,)Mijmsijmsnuu 当插值函数为高次多项式时，用自然坐标表示的单元(,)()()()tuvijmtiujv
5、mPPP ()kP z11()()KkwNzwP zw1K 10K 一阶基函数二阶基函数三阶基函数2222()()()kpsuuJ uku dxdyxy222211()()()ttllkpeettppuudxdyku dxdyuxyu21 1221 122.(.)pppnnkppppnnK uK uK ukB uB uB u1()()()2ttppesspsietjmjmKdxdy1()()()2ttppesspsjetmimiKdxdy1()()()2ttppesspsmetijijKdxdy()()()tttteeeepspsiipsjjpsmmKKctgKctgKctgttepspseB
6、dxdy 经过化简其中可以事先制定成表格，程序中调用10H20H30H40H2kk p31.87 10375.1 10398.9 103118.0 10kpkp波型 (m)0.0460.0230.0200.01829解析解 (m)0.0460.0230.0200.01834误差0.02%0.35%0.08%0.26%波导尺寸为10mm*23mm。分格数为4*3，应用2阶FEM分析得到特征值如下，三角单元总数为24，场节点总数为63 谐振腔尺寸为12mm*28mm时，应用FEM进行分析，划分网格数为5*4，三角单元总数为40，场节点总数为99。10H20H30H40H2kkp32.06*1038
7、3.2*103106.9*103130.3*10kpkp波型 (m)0.0560.0280.0240.02204解析解 (m)0.0560.0280.0240.02207误 差0.11%0.15%0.25%0.15%01234567805101520253035f(GHz)1stTDFEM 2ndTDFEM Analytical 1stFDFEM 2ndFDFEMTE10 TE20 TE01 TE11 TE30 TE40 TE12 TE22特征值特征值0123456789-25-20-15-10-505101520TE10 TE20 TE01 TE11 TE30 TE40 TE12 TE22e
8、rror%1stTDFEM 2ndTDFEM 1stFDFEM 2ndFDFEM计算误差计算误差5、三维节点有限元方法、三维节点有限元方法6、三维高阶节点有限元方法、三维高阶节点有限元方法 略7、节点有限元的缺点、节点有限元的缺点标量有限元法的缺标量有限元法的缺点：点：计算结果中会出计算结果中会出现不易辨别的非现不易辨别的非物理解或伪解；物理解或伪解；在介质分界面或在介质分界面或导体表面不易强导体表面不易强加边界条件；加边界条件；由于场存在奇异由于场存在奇异性，处理导体及性，处理导体及介质边缘和棱边介质边缘和棱边有困难。有困难。9 二维矢量有限元方法二维矢量有限元方法 用棱边元矢量基函数来表示
9、矢量场有以下特点：（）它能自然满足电场或磁场在介质分界面上的切向连续条件；（）能很容易表示金属表面的切向场条件，在介质、金属边缘或棱角处，也能很好地表示场分布；（）自然满足零散度条件，消除了非物理解；（）由于边棱与边棱之间的耦合弱于节点与节点之间的耦合，所得到的总体矩阵具有较少的非零元，减少了计算量。10 三维矢量有限元方法三维矢量有限元方法其中N为矢量基函数。将(2)式代入(1)式中得到：)(20121eeTeeeTeMeEBEkEAEF 其中：dVNNATeeV1eer dVNNBTeeVeree 应用里兹方法，即对每一个未知棱边场求偏导数，并令其等于零，就可得到本征值方程组 AE=20k
10、 BE 在腔体壁上强加边界条件，解此本征值方程组，即可求出本征值以及对应的本征矢。表 2 对半径为 0.5cm，长为 0.5cm 的金属柱谐振腔计算得到的八个最低非零本征值(0k,1cm),网格剖分为2035 模式 解析解 计算值 误差 010TM 4.810 4.8018 0.170 111TE 7.283 7.3335 7.3335 0.693 0.693 110TM 7.650 7.6280 7.6280 0.288 0.288 011TM 7.840 8.0306 2.431 211TE 8.658 8.9665 8.9665 3.563 3.563 表 5.为半径为 0.5cm，长
11、0.5cm 的金属圆柱腔内底部中心加载一内半径为 0.2cm，外半径为 0.4cm，长 0.16667cm 的介质环，计算出的五个最低本振值(0k,1cm),(16r)模式 波数 1 3.066 2 3.365 3 3.365 4 3.727 5 3.372 0.167 0.4 0.8 表6给出了半径为0.5cm长0.6cm的金属圆柱谐振腔内一端加载内半径为 0.2cm，外半径为 0.4cm，长 0.2cm 的介质环，另一端接一细长金属竿(长 0.2cm，半径 0.1cm),16r 计算出的五个最低本征值(0k,1cm),0.2 0.4 0.2 0.2 0.2 0.1 模式 波数 1 3.43
12、35 2 3.4335 3 3.6534 4 3.6534 5 3.7092 结 论 虽然在用矢量有限元分析问题时不会出现伪解，但是观察到零本征值的存在，而零本征值不对应与物理模式，故应去除。不过,零本征值很容易识别，且可事先预测，所以不成问题。通过以上的数据结果对比与分析可以发现，应用矢量有限元对谐振腔进行数值模拟是非常方便和准确的，这对谐振腔的设计，应用提供了又一非常好的指导途径。矢量有限元性态差矢量有限元性态差11、高阶矢量基函数、高阶矢量基函数略12、有限元中的现代技术、有限元中的现代技术 PML Method DDM Method FEM-BI FEM-PO，IPO etc.13、时
13、域有限元方法、时域有限元方法(一)研究背景时域数值方法时域有限差分法(FDTD)传输线矩阵法(TLM)时域积分方程法(TDIE)时域有限元法(TDFEM)多分辨率时域技术(MRTD)其他时域伪谱方法(PSTD)TDFEM 研究现状J.M.Jin洪 伟金亚秋南京理工大学电子科技大学南京航空航天大学(二)理论分析 对于一般的时变场，在电流源激励下，Maxwell两个旋度方程为：Mt BEt DHJ 由上述两式可以推出关于H的矢量波动方程：22tHHJ 该矢量波动方程的等价泛函为：22()()()()FdtdS HHHHHH JHHn 使用非结构网格离散计算区域，使用矢量基函数Ni()对各个单元内的
14、磁场H进行插值：1,NjjjtutHN 在完纯导体表面满足第二类边界条件：()0 nH 应用里兹方法进行处理，最终可以得到常微分方程：220duTS ufd t,ijijTN N1,ijijSNN 其中，1,if NJ 采用恒稳的Newmark方法离散时间步，得到如下差分方程：212221()()(2(12)()()nnnTStuftTStuTStu 其中，t为时间步长，为Newmark离散参数，un+1为n+1时间步的待求向量，un，un-1是n,n-1时间步的已求向量(三三)数据分析数据分析22200()2()/)exp()/)ttttt J其中中心时间t0=25.9ns，脉冲宽度=4ns
15、 观察点在 0.1250.1125xyxyo1m(40)1m(40)020406080100120-4x107-3x107-2x107-1x10701x1072x1073x1074x107Ey(V/m)t(ns)空 波 导空 波 导图图1 空谐振器中的磁场响应曲线空谐振器中的磁场响应曲线020406080100120-6.0 x107-4.0 x107-2.0 x1070.02.0 x1074.0 x1076.0 x107Hy(V/m)t(ns)图图2 部分介质加载谐振器中磁场响应曲部分介质加载谐振器中磁场响应曲线线 xyo0.55m0.45mr=9.00.00.10.201x1062x106
16、3x1064x1065x1066x106Frequency(GHz)|Hy|(V/m)图图3 部分介质加载谐振器中响应频谱曲线部分介质加载谐振器中响应频谱曲线 0.14038GHz0.20142GHz020406080100120-8.0 x108-6.0 x108-4.0 x108-2.0 x1080.02.0 x1084.0 x1086.0 x1088.0 x108Hy(V/m)t(ns)图图4 全部介质加载全部介质加载(r=2.25)时磁场响应曲线时磁场响应曲线0.00.10.20.30.40.06.0 x1071.2x1081.8x108Frequency(GHz)|Hy|(V/m)图
17、图5 全部介质加载时磁场响应频谱曲线全部介质加载时磁场响应频谱曲线 0.14038GHz0.22583GHz0.28076GHz0.31128GHz0.36011GHz模式精确值(GHz)仿真值(GHz)误差(%)TM110.141420.140380.735TM210.223610.225830.993TM220.282840.280760.735TM130.316230.311280.125TM230.360560.360110.1250102030405060708090 100 1101E-81E-71E-61E-51E-41E-30.010.11误差误差迭代迭代图图6 CG收敛曲线对
18、比收敛曲线对比14 FEM求解器求解器-CG等同于最小余量法了一般频域FEM的系数矩阵不具备自伴性位了求解（11.53）泛函得到如下线性方程组但是一、基函数p未知，二、我们为了求原来的方程组，必须求解新的方程组，无意义初始剩初始剩余向量余向量向量系向量系数数构造新构造新解解构造新的构造新的剩余向量剩余向量构造新的构造新的基函数向基函数向量量Comments on CG 迭代方法-理论上N步收敛 单调收敛 与矩阵性态有关-条件数-谱半径-最大特征值与最小特征值比 预处理技术-SSOR-ICCG 改变内积的定义CG和BICG方法的求解速度与矩阵的性态即条件数由很大关系。通过预处理技术可以有效提高矩阵的性态，从而提高计算速度
