% smt2djxfwd01.m
% 二維大地電磁有限元數(shù)值模擬 矩形+線性插值
clear
close all
clc
format short g
%
% 讀頻率數(shù)據(jù)
%%?

[nfreq,freq]=sinfreq('freq.dat');
%
% 讀模型數(shù)據(jù)
%
[idx,ney,nez,dy,dz,res,resbottom]=sinmt2dmodjx1('mod2dfem1')

t1=cputime;
for ifreq=1:nfreq
? ? ifreq;
? ? frek=freq(ifreq);
? ? [gk,s]=sglobalksrectlinear(frek,idx,ney,nez,dy,dz,res,resbottom);
? ? u=gk\s;
? ? %
? ? % 計(jì)算視電阻率
? ? %
? ? [es,hs,ksair] = sfldsurf(idx,frek,u,ney,nez,dy,dz,res);
? ? [z,r,p]=sImpsurf(idx,frek,es,hs,ney);
? ? zimp(ifreq,1:ney+1) = z(1:ney+1,1);
? ? rhoa(ifreq,1:ney+1) = r(1:ney+1,1);
? ? phaz(ifreq,1:ney+1) = p(1:ney+1,1);
end
t2=cputime-t1
stemp = sprintf('花費(fèi)時(shí)間: %g s',t2);
disp(stemp);

clear x
x(1)=0;
for iy=2:ney+1
? ? x(iy) = x(iy-1)+dy(iy-1);
end
for iy=1:ney+1
? ? x(iy) = x(iy) - x(ney+1)/2;
end

clear x1
for i=6:16
? ? x1(i-6+1) = x(i)/1000;
end
clear x
x = x1;
nx = length(x);

clear y
for i=1:nfreq
? ? y(i) = log10(freq(i));
end

clear x1 y1 z1
for i=1:nx
? ? for j=1:nfreq
? ? ? ? x1(i,j) = x(i);
? ? ? ? y1(i,j) = y(j);
? ? ? ? z1(i,j) = log10(rhoa(j,i-1+6));
? ? end
end

figure(1);
colormap(jet);
contourf(x1,y1,z1);
title('TM模式的視電阻率擬斷面圖');
xlabel( ' x /km');
ylabel(' log10(f) /Hz');
colorbar('vertical');

clear y
y(1)=0;
for i=2:nez+1
? ? y(i) = y(i-1) + dz(i-1);
end
for i=1:nez+1
? ? y(i) = y(i)/1000;
end

clear x1 y1 z1;
for i=1:nx
? ? for j=1:nez-6
? ? ? ? x1(i,j) = x(i);
? ? ? ? y1(i,j) = -y(j);
? ? ? ? z1(i,j) = res(j,i-1+6);
? ? end
end

figure(2);
colormap(jet);
pcolor(x1,y1,z1);
title('TM模式網(wǎng)格剖分示意圖');
xlabel(' x/km');
ylabel(' z/km');
axis equal;
axis tight;
% colorbar('vertical');

D140

?

總結(jié)

以上是生活随笔為你收集整理的二维大地电磁有限元数值模拟矩形+线性插值的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。