MoM矩量法(二):阻抗矩阵的填充

MoM矩量法(二):阻抗矩阵的填充

阻抗矩阵的填充是矩量法的核心,也是矩量法中最难的一部分。本文采用的是基于RWG基函数的矩量法,公式的推导不再赘述,直接贴上公式讲代码思路。
这里采用的公式来自于聂在平主编的《目标与环境 电磁散射特性建模——理论、方法与实现(应用篇)》,所有数据均采用国际标准单位。

完整代码已贴在文章最后。文章和代码思路借鉴了https://blog.csdn.net/u014411646/article/details/99689460,谢谢原作者!

由于阻抗矩阵是在两个三角单元上积分而来,当两个三角单元相同时,就会出现积分的奇异性问题,所以阻抗矩阵的计算包括奇异积分和非奇异积分两部分。

一、非奇异积分

在这里插入图片描述

上面是阻抗矩阵的填充公式,我们以 m m m n n n代表两条边,字母右上角的正负号表示的是左右两个三角单元。例如, Z m n + − Z_{mn}^{+-} Zmn+就表示对 m m m边的左边三角单元和 n n n边的右边三角单元进行运算。

上式中各符号的含义:
  1. k k k代表自由空间波数,值为 2 π / λ 2\pi/\lambda 2π/λ,单位为 r a d / m rad/m rad/m Z Z Z代表自由空间波阻抗,值为 376.73031 376.73031 376.73031欧姆;
  2. l m l_m lm l n l_n ln代表两条边的长度;
  3. w k 、 w j w_k、w_j wkwj x k 、 y k 、 z k 、 x j 、 y j 、 z j x_k、y_k、z_k、x_j、y_j、z_j xkykzkxjyjzj分别为采样点加权因子和三角形数值积分采样点参数值,具体取值在后面给出;
  4. p m i p_{mi} pmi p n i p_{ni} pni代表三角单元的三个顶点。
    下表为三角单元数值积分表:
    在这里插入图片描述
核心代码:
for (int k = 0;k < 4;k++) {
		rowvec rmk(3);
		rmk(0) = ptm(0, 0) * xx[k] + ptm(1, 0) * yy[k] + ptm(2, 0) * zz[k];
		rmk(1) = ptm(0, 1) * xx[k] + ptm(1, 1) * yy[k] + ptm(2, 1) * zz[k];
		rmk(2) = ptm(0, 2) * xx[k] + ptm(1, 2) * yy[k] + ptm(2, 2) * zz[k];
		rowvec P = rmk - ptm.row(0);
		for (int j = 0;j < 4;j++) {
			rowvec rnj(3);
			rnj(0) = ptn(0, 0) * xx[j] + ptn(1, 0) * yy[j] + ptn(2, 0) * zz[j];
			rnj(1) = ptn(0, 1) * xx[j] + ptn(1, 1) * yy[j] + ptn(2, 1) * zz[j];
			rnj(2) = ptn(0, 2) * xx[j] + ptn(1, 2) * yy[j] + ptn(2, 2) * zz[j];
			rowvec Q = rnj - ptn.row(0);
			double R = getDist(rmk, rnj);
			Z += ww[k] * ww[j] * (getDot(P, Q) - 4. / (k0 * k0)) * exp(-J * k0 * R) / R;
		}
	}

代码与公式的联系比较直观。 r m k rmk rmk r n j rnj rnj为3列的行矩阵,储存的便是上面公式中的 r m k r_{mk} rmk r n j r_{nj} rnj;而 P P P Q Q Q分别为 ( r m k − p m 1 ) (r_{mk}-p_{m1}) (rmkpm1) ( r n j − p n 1 ) (r_{nj} - p_{n1}) (rnjpn1) R R R即为上式中的 ∣ r m k − r n j ∣ |r_{mk}-r_{nj}| rmkrnj

二、奇异积分

Z m n + − Z_{mn}^{+-} Zmn+的公式就可以看到,它的计算其实就是每次从四个三角单元里面拎出来两个进行运算的过程,而当被拎出来的两个三角单元为同一个三角形时,就会出现奇异问题,这时就要进行相应的奇异积分。
在这里插入图片描述
上式即为奇异值的计算公式。下面来说明一下各符号的含义:

  1. l m 、 l n l_m、l_n lmln w k w_k wk等符号与上一节的一致;
  2. R R R也为上一节中的 ∣ r m k − r n j ∣ |r_{mk}-r_{nj}| rmkrnj
  3. A n A_n An n n n边所对应的三角单元面积, r n c r_{nc} rnc n n n边所对应的三角单元中心;
  4. P 、 R P、R PR l l l的公式见下面代码。
核心代码:
for (int k = 0;k < 7;k++) {
		rowvec rmk(3);
		rmk(0) = ptm(0, 0) * xxx[k] + ptm(1, 0) * yyy[k] + ptm(2, 0) * zzz[k];
		rmk(1) = ptm(0, 1) * xxx[k] + ptm(1, 1) * yyy[k] + ptm(2, 1) * zzz[k];
		rmk(2) = ptm(0, 2) * xxx[k] + ptm(1, 2) * yyy[k] + ptm(2, 2) * zzz[k];
		rowvec P = rmk - ptm.row(0);
		Complex temp1, temp2;
		temp1 = 0;
		for (int j = 0;j < 7;j++) {
			rowvec rnj(3);
			rnj(0) = ptn(0, 0) * xxx[j] + ptn(1, 0) * yyy[j] + ptn(2, 0) * zzz[j];
			rnj(1) = ptn(0, 1) * xxx[j] + ptn(1, 1) * yyy[j] + ptn(2, 1) * zzz[j];
			rnj(2) = ptn(0, 2) * xxx[j] + ptn(1, 2) * yyy[j] + ptn(2, 2) * zzz[j];
			rowvec Q = rnj - ptn.row(0);
			double R = getDist(rmk, rnj);
			temp1 += www[j] * (getDot(P, Q) - 4. / (k0 * k0)) *
				(-k0 * k0 * R / 2. + J * k0 / 6. * (k0 * k0 * R * R - 6));
		}
		temp2 = 0;
		for (int j = 0;j < 3;j++) {
			rowvec lj(3);
			rowvec nj0(3);
			lj = getUnit(ptn.row((j + 1) % 3) - ptn.row(j));
			nj0 = ptn.row(j) + lj * getDot(rmk - ptn.row(j), lj);
			double lmj = getDot(ptn.row(j) - rmk, lj);
			double lpj = getDot(ptn.row((j + 1) % 3) - rmk, lj);
			double rmj = getDist(rmk, ptn.row(j));
			double rpj = getDist(rmk, ptn.row((j + 1) % 3));
			double pj0 = getDist(nj0, rmk);
			if (pj0 > 1.e-15) {
				temp2 += (getDot(P, ptnc - ptn.row(0)) - 4. / (k0 * k0)) * pj0 * log((rpj + lpj) / (rmj + lmj));
			}
		}
		Z += www[k] * (temp1 + temp2 / (2. * TriArea(n)));
	}

与上一节相似,代码与公式的联系比较紧密,所以这里就不再赘述了。

下面贴上所有代码:
/***********************************************************/
/*            该函数对阻抗矩阵进行填充                      */
/*                                                         */
/***********************************************************/
/*                                                         */
/*             Author:Tianshaowenwen                      */
/*                                                         */
/***********************************************************/
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rwg2.h"
#include "armadillo"


using namespace arma;
using namespace std;


extern Complex J;
extern mat Edge_list;//[端点1, 端点2, 对点1, 对点2, 左三角单元, 右三角单元]
extern mat TriCenter;
extern vec EdgeLength;
extern vec TriArea;
extern int TriTotal, EdgeTotal;
extern mat pdata, tdata;
extern double k0, Z0, PI;
extern double xx[4];
extern double yy[4];
extern double zz[4];
extern double ww[4];
extern double xxx[7];
extern double yyy[7];
extern double zzz[7];
extern double www[7];

cx_mat Zmn;



void rwg2_main() {
	rwg1_main();
	complex<double> Zll, Zlr, Zrl, Zrr;
	int ml1, ml2, ml3, mr1, mr2, mr3;
	int nl1, nl2, nl3, nr1, nr2, nr3;
	double lm, ln;
	cx_mat zmn(EdgeTotal, EdgeTotal);
	cout << "\n" << "填充阻抗矩阵..." << endl;
	clock_t start = clock();
	for (int i = 0;i < EdgeTotal;i++) {
		int ml = Edge_list(i, 4);
		int mr = Edge_list(i, 5);
		ml1 = Edge_list(i, 2);
		mr1 = Edge_list(i, 3);
		ml2 = mr3 = Edge_list(i, 0);
		ml3 = mr2 = Edge_list(i, 1);
		lm = EdgeLength(i);
		for (int j = 0;j <= i;j++) {
			int nl = Edge_list(j, 4);
			int nr = Edge_list(j, 5);
			nl1 = Edge_list(j, 2);
			nr1 = Edge_list(j, 3);
			nl2 = nr3 = Edge_list(j, 0);
			nl3 = nr2 = Edge_list(j, 1);
			ln = EdgeLength(j);
			if (ml == nl) {
				Zll = getZmnSingular(ml1, ml2, ml3, nl1, nl2, nl3, nl);
			}
			else {
				Zll = getZmnGeneral(ml1, ml2, ml3, nl1, nl2, nl3);
			}
			if (ml == nr) {
				Zlr = getZmnSingular(ml1, ml2, ml3, nr1, nr2, nr3, nr);
			}
			else {
				Zlr = getZmnGeneral(ml1, ml2, ml3, nr1, nr2, nr3);
			}
			if (mr == nl) {
				Zrl = getZmnSingular(mr1, mr2, mr3, nl1, nl2, nl3, nl);
			}
			else {
				Zrl = getZmnGeneral(mr1, mr2, mr3, nl1, nl2, nl3);
			}
			if (mr == nr) {
				Zrr = getZmnSingular(mr1, mr2, mr3, nr1, nr2, nr3, nr);
			}
			else {
				Zrr = getZmnGeneral(mr1, mr2, mr3, nr1, nr2, nr3);
			}
			zmn(i, j) = J * k0 * Z0 * lm * ln * (Zll + Zrr - Zlr - Zrl) / (4 * PI);
			if (i != j) {
				zmn(j, i) = zmn(i, j);
			}
		}
	}
	clock_t end = clock();
	cout << "\n" << "阻抗矩阵填充时间: " << (end - start) / 1000.0 << "s" << endl;
	zmn = luDcmp(zmn);
	/*zmn.save("Zmn.txt", raw_ascii);*/
	Zmn = zmn;
}

/***********************************************************/
/*            Function:Complex getZmnSingular()           */
/***********************************************************/
/* 1.计算三角单元奇异积分                                   */
/***********************************************************/
Complex getZmnSingular(int m1, int m2, int m3, int n1, int n2, int n3, int n) {
	Complex Z = 0;
	mat ptm(3, 3), ptn(3, 3);
	ptm.row(0) = pdata.row(m1);
	ptm.row(1) = pdata.row(m2);
	ptm.row(2) = pdata.row(m3);
	ptn.row(0) = pdata.row(n1);
	ptn.row(1) = pdata.row(n2);
	ptn.row(2) = pdata.row(n3);
	rowvec ptnc = TriCenter.row(n);
	for (int k = 0;k < 7;k++) {
		rowvec rmk(3);
		rmk(0) = ptm(0, 0) * xxx[k] + ptm(1, 0) * yyy[k] + ptm(2, 0) * zzz[k];
		rmk(1) = ptm(0, 1) * xxx[k] + ptm(1, 1) * yyy[k] + ptm(2, 1) * zzz[k];
		rmk(2) = ptm(0, 2) * xxx[k] + ptm(1, 2) * yyy[k] + ptm(2, 2) * zzz[k];
		rowvec P = rmk - ptm.row(0);
		Complex temp1, temp2;
		temp1 = 0;
		for (int j = 0;j < 7;j++) {
			rowvec rnj(3);
			rnj(0) = ptn(0, 0) * xxx[j] + ptn(1, 0) * yyy[j] + ptn(2, 0) * zzz[j];
			rnj(1) = ptn(0, 1) * xxx[j] + ptn(1, 1) * yyy[j] + ptn(2, 1) * zzz[j];
			rnj(2) = ptn(0, 2) * xxx[j] + ptn(1, 2) * yyy[j] + ptn(2, 2) * zzz[j];
			rowvec Q = rnj - ptn.row(0);
			double R = getDist(rmk, rnj);
			temp1 += www[j] * (getDot(P, Q) - 4. / (k0 * k0)) *
				(-k0 * k0 * R / 2. + J * k0 / 6. * (k0 * k0 * R * R - 6));
		}
		temp2 = 0;
		for (int j = 0;j < 3;j++) {
			rowvec lj(3);
			rowvec nj0(3);
			lj = getUnit(ptn.row((j + 1) % 3) - ptn.row(j));
			nj0 = ptn.row(j) + lj * getDot(rmk - ptn.row(j), lj);
			double lmj = getDot(ptn.row(j) - rmk, lj);
			double lpj = getDot(ptn.row((j + 1) % 3) - rmk, lj);
			double rmj = getDist(rmk, ptn.row(j));
			double rpj = getDist(rmk, ptn.row((j + 1) % 3));
			double pj0 = getDist(nj0, rmk);
			if (pj0 > 1.e-15) {
				temp2 += (getDot(P, ptnc - ptn.row(0)) - 4. / (k0 * k0)) * pj0 * log((rpj + lpj) / (rmj + lmj));
			}
		}
		Z += www[k] * (temp1 + temp2 / (2. * TriArea(n)));
	}
	return Z;
}

/***********************************************************/
/*            Function:Complex getZmnGeneral()            */
/***********************************************************/
/* 1.计算三角单元非奇异积分                                 */
/***********************************************************/
Complex getZmnGeneral(int m1, int m2, int m3, int n1, int n2, int n3) {
	Complex Z = 0;
	mat ptm(3, 3), ptn(3, 3);
	ptm.row(0) = pdata.row(m1);
	ptm.row(1) = pdata.row(m2);
	ptm.row(2) = pdata.row(m3);
	ptn.row(0) = pdata.row(n1);
	ptn.row(1) = pdata.row(n2);
	ptn.row(2) = pdata.row(n3);
	for (int k = 0;k < 4;k++) {
		rowvec rmk(3);
		rmk(0) = ptm(0, 0) * xx[k] + ptm(1, 0) * yy[k] + ptm(2, 0) * zz[k];
		rmk(1) = ptm(0, 1) * xx[k] + ptm(1, 1) * yy[k] + ptm(2, 1) * zz[k];
		rmk(2) = ptm(0, 2) * xx[k] + ptm(1, 2) * yy[k] + ptm(2, 2) * zz[k];
		rowvec P = rmk - ptm.row(0);
		for (int j = 0;j < 4;j++) {
			rowvec rnj(3);
			rnj(0) = ptn(0, 0) * xx[j] + ptn(1, 0) * yy[j] + ptn(2, 0) * zz[j];
			rnj(1) = ptn(0, 1) * xx[j] + ptn(1, 1) * yy[j] + ptn(2, 1) * zz[j];
			rnj(2) = ptn(0, 2) * xx[j] + ptn(1, 2) * yy[j] + ptn(2, 2) * zz[j];
			rowvec Q = rnj - ptn.row(0);
			double R = getDist(rmk, rnj);
			Z += ww[k] * ww[j] * (getDot(P, Q) - 4. / (k0 * k0)) * exp(-J * k0 * R) / R;
		}
	}
	return Z;
}
/***********************************************************/
/*            Function:cx_mat luDcmp(cx_mat zmn)          */
/***********************************************************/
/* 1.对复矩阵进行LU分解                                     */
/***********************************************************/
cx_mat luDcmp(cx_mat zmn) {
	int i, j, k;
	Complex dum, sum, temp;
	for (j = 0;j < EdgeTotal;j++) {
		for (i = 0;i < j;i++) {
			sum = zmn(i, j);
			for (k = 0;k < i;k++) {
				sum -= zmn(i, k) * zmn(k, j);
			}
			zmn(i, j) = sum;
		}
		for (i = j;i < EdgeTotal;i++) {
			sum = zmn(i, j);
			for (k = 0;k < j;k++) {
				sum -= zmn(i, k) * zmn(k, j);
			}
			zmn(i, j) = sum;
		}
		if ((zmn(j, j).real() < 1.e-20) && (zmn(j, j).real() > -1.e-20) &&
			(zmn(j, j).imag() < 1.e-20) && (zmn(j, j).imag() > -1.e-20)) {
			zmn(j, j) = 1.e-20;
		}
		if (j != EdgeTotal) {
			dum = 1.0 / zmn(j, j);
			for (i = j + 1;i < EdgeTotal;i++) {
				zmn(i, j) *= dum;
			}
		}
	}
	return zmn;
}

代码中还包含一部分是对复矩阵的LU分解,这是为了后面解矩阵方程时提高效率。如果你是用已经写好的库函数来解矩阵方程的话,这一部分就可以忽略掉。

  • 7
    点赞
  • 21
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值