clear all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% 大地电磁二维限元正演(TE模式)
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u=4*pi*10^(-7); % 空气中的磁导率
%Z0=[10000,5000,2000,1000,500]; % 空气层
Z0=[10000,5000,3000,2000,1000,500,200,100];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 数据读取 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=load('frequency'); %%
node=load('nodemodel'); %%
survey=load('surveymodel'); sp=survey; [Np,cc]=size(sp); %%
body=load('bodymodle'); sc=struct2cell(body); p=cell2mat(sc); %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 网格参数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=node(2:node(1,1)+1,1)'; %%
Z=node(2:node(1,2)+1,2)'; Z=[Z0,Z]; NK=length(Z0); %%
[cc,Nx]=size(X); [cc,Nz]=size(Z); [Nf,cc]=size(f); %%
NP=(Nx+1)*(Nz+1); %%
Kz1=sparse(NP,NP); Kz2=sparse(NP,NP); Kz3=sparse(NP,NP); %%
[cb,cc]=size(body); [cs,cc]=size(survey); %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 求刚度系数矩阵 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w0=2*pi*1;
pp=100; % 半空间(背景电阻率值)
tic
t=1/(i*w0*u);
p=[10^4*ones(NK,Nx);p];
for h=1:Nx
for j=1:Nz
k=1/p(j,h);
a=(h-1)*(Nz+1)+j; b=(h-1)*(Nz+1)+j+1;
c=(h-1)*(Nz+1)+j+(Nz+1)+1; e=c-1;
K1=Ke1(X(h),Z(j),t); K2=Ke2(X(h),Z(j),k);
Kz1=Kz1+KK1(a,b,c,e,K1,NP);
Kz2=Kz2+KK1(a,b,c,e,K2,NP);
if(j==Nz)
kk=sqrt((-i*w0*u)/pp);
K3=Ke3(Z(j),kk,t);
Kz3=Kz3+KK1(a,b,c,e,K3,NP);
end
end
end
toc
Ue=zeros(NP,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
for r=1:Nf
B=zeros(1,NP);
w=2*pi*f(r);
K=(1/f(r))*Kz1-Kz2+sqrt(f(r))*(1/f(r))*Kz3;
for x=1:(Nx+1)
y=(x-1)*(Nz+1)+1; K(y,y)=K(y,y)*10^10; B(1,y)=K(y,y)*10^10;
end
[L1,U1]=luinc(K,1);
Ue=bicgstab(K,B',1e-15,800,L1,U1,Ue); % 不完全LU分解,稳定双共轭梯度法求解线性方程组
L=Z(1+NK)+Z(2+NK)+Z(3+NK);
for m=1:Np
n=(sp(m)-1)*(Nz+1)+1+NK;
de(m)=(1/(2*L))*((-11)*Ue(n,1)+18*Ue(n+1,1)+(-9)*Ue(n+2,1)+2*Ue(n+3,1));
pe(r,m)=abs(-i*w*u*(Ue(n,1)/de(m))^2);
qe(r,m)=90+atan(imag(Ue(n,1)*i*w*u/de(m))/real(Ue(n,1)*i*w*u/de(m)))*(360/(2*pi));
end
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 保存数据 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
result=[Nf*Np,Nf,Np];
for j=1:Np
result=[result;f,pe(:,j),qe(:,j)];
end
D1=fopen('resultTE.dat','wt');
fprintf(D1,'%.6d %.2d %.2d\n',result');
fclose(D1); % 结果
% 保存结果中第一行数据:第一个为数据个数,第二个位频点数,第三个为测点数
% 从第二行开始,第一列为频率,第二列为电阻率,第三列为相位
没有合适的资源?快使用搜索试试~ 我知道了~
资源推荐
资源详情
资源评论
收起资源包目录
MT2D正演.zip (36个子文件)
MT2D正演
bodymodle.mat 191B
frequency 170B
plotMT2Dmodel.m 2KB
nodemodel 619B
Ke2.m 106B
MT2DMODEL.m 1KB
KK1.m 484B
resultTE.dat 934B
MT2DMODEL.asv 1KB
TEmodel.m 3KB
resultTH.dat 606B
Ke3.m 103B
TMmodel.m 3KB
MT2Dmesh.m 2KB
2Dmodel.fig 38KB
surveymodel 4B
Ke1.m 261B
untitled.fig 14KB
MT2D
bodymodle.mat 191B
frequency 170B
plotMT2Dmodel.m 2KB
nodemodel 619B
Ke2.m 106B
MT2DMODEL.m 1KB
KK1.m 484B
resultTE.dat 934B
MT2DMODEL.asv 1KB
TEmodel.m 3KB
resultTH.dat 606B
Ke3.m 103B
TMmodel.m 3KB
MT2Dmesh.m 2KB
2Dmodel.fig 38KB
surveymodel 4B
Ke1.m 254B
untitled.fig 14KB
共 36 条
- 1
资源评论
- xu123xu12232021-07-07就是注释不全,对于新学者来说不是太清楚
- 安好。。。2023-02-26#内容缺失
Monkey_shuo
- 粉丝: 12
- 资源: 1
上传资源 快速赚钱
- 我的内容管理 展开
- 我的资源 快来上传第一个资源
- 我的收益 登录查看自己的收益
- 我的积分 登录查看自己的积分
- 我的C币 登录后查看C币余额
- 我的收藏
- 我的下载
- 下载帮助
最新资源
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
安全验证
文档复制为VIP权益,开通VIP直接复制
信息提交成功