電阻率法數值模擬實驗和激發極化法數字模擬實驗報告
求matlab代碼
資料越簡單易懂越好
uj5u.com熱心網友回復:
哈哈你也是新學的吧,這里是我的代碼,探討一下。dy=1
ny=101
ymin=-50
dx=1
nx=101
xmin=-50
x=xmin:dx:(xmin+(nx-1)*dx)
y=ymin:dy:(ymin+(ny-1)*dy)
[X,Y]=meshgrid(x,y)
p1=10
u=0.02
r=5
h=9
%中間梯度裝置
ps=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2)./(h^2+x.^2).^5/2]
plot(x,ps),title('中間梯度主刨面例外')
grid on;
hold on
j=0
figure
for y=-5:1:5
j=j+1
ps=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2+y^2)./(h^2+x.^2+y^2).^5/2]
q(j)=plot(x,ps)
hold on
end
grid on;
title('中間梯度隨刨面位置變化')
figure
ps=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*X.^2+Y.^2)./(h^2+X.^2+Y.^2).^5/2]
pcolor(X,Y,ps);shading interp
title('中間梯度等值線圖')
% 聯合刨面
MN=4
AM=10
AN=14
BN=10
BM=14
dx=1
nx=101
x=xmin:dx:(xmin+(nx-1)*dx)
KK=AM*AN/MN
da=((x-2-AM).^2+h^2).^(1/2)
db=((x+2+BN).^2+h^2).^(1/2)
rm=((x-2).^2+h^2).^(1/2)
rn=((x+2).^2+h^2).^(1/2)
cosa=(da.^2+rm.^2-AM^2)./(da.*rm.*2)
cosb=(da.^2+rn.^2-AN^2)./(da.*rn.*2)
psA=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(cosa./(da.^2+rm.^2)-cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN^2)./(db.*rn.*2)
psB=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(-cosc./(db.^2+rm.^2)+cosd./(db.^2+rn.^2))]
figure
xx=plot(x,psA)
hold on
yy=plot(x,psB)
legend("psA","psB")
title("聯合刨面低阻球體上方視電阻率區線")
%對稱四極
figure
for x=0:1:100
AB=20:10:1000
MN=4
AM=(AB-MN)/2
BN=(AB-MN)/2
AN=AM+MN
BM=BN+MN
K=AM.*AN/MN
da=((x-2-AM).^2+h^2).^(1/2)
db=((x+2+BN).^2+h^2).^(1/2)
rm=((x-2).^2+h^2)^(1/2)
rn=((x+2).^2+h^2)^(1/2)
cosa=(da.^2+rm^2-AM.^2)./(2*da.*rm)
cosb=(da.^2+rn^2-AN.^2)./(2*da.*rn)
psA1=p1*[1+2*(u-1)./(2*u+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM.^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN.^2)./(db.*rn.*2)
psB1=p1*[1+2*(u-1)/(2*u+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
psAB=(psA1+psB1)./2
plot(log10(AB/2),psAB)
hold on
end
title("視電阻率測深曲線")
figure
x=-50:1:50
ab=10:4:100;
[AB,X]=meshgrid(ab,x)
MN=4
AM=(AB-MN)/2
BN=(AB-MN)/2
AN=AM+MN
BM=BN+MN
K=AM.*AN/MN
da=((X-2-AM).^2+h^2).^(1/2)
db=((X+2+BN).^2+h^2).^(1/2)
rm=((X-2).^2+h^2).^(1/2)
rn=((X+2).^2+h^2).^(1/2)
cosa=(da.^2+rm.^2-AM.^2)./(2*da.*rm)
cosb=(da.^2+rn.^2-AN.^2)./(2*da.*rn)
psA1=p1*[1+2*(u-1)./(2*u+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM.^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN.^2)./(db.*rn.*2)
psB1=p1*[1+2*(u-1)/(2*u+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
psAB=(psA1+psB1)./2
pcolor(X,log10(AB./2),psAB),shading interp
set(gca,'YDir','reverse');
colorbar
xlabel('X(m)'),ylabel('AB/2(m)');
title("擬斷面圖")
%這是數值模擬實驗
% 聯合刨面
r=5
h=10
dy=1
ny=101
ymin=-50
dx=1
nx=101
xmin=-50
p1=10
p2=8
MN=1
AM=10
AN=14
BN=10
BM=14
dx=1
nx=101
n=0.3
u=p2/p1
y=ymin:dy:(ymin+(ny-1)*dy)
x=xmin:dx:(xmin+(nx-1)*dx)
[X,Y]=meshgrid(x,y)
KK=AM*AN/MN
da=((x-0.5-AM).^2+h^2).^(1/2)
db=((x+0.5+BN).^2+h^2).^(1/2)
rm=((x-0.5).^2+h^2).^(1/2)
rn=((x+0.5).^2+h^2).^(1/2)
cosa=(da.^2+rm.^2-AM^2)./(da.*rm.*2)
cosb=(da.^2+rn.^2-AN^2)./(da.*rn.*2)
psA=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(cosa./(da.^2+rm.^2)-cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN^2)./(db.*rn.*2)
psB=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(-cosc./(db.^2+rm.^2)+cosd./(db.^2+rn.^2))]
p3=p2/(1-n)
u=p3/p1
psA2=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(cosa./(da.^2+rm.^2)-cosb./(da.^2+rn.^2))]
psB2=p1*[1+2*KK*(u-1)/(2*u+1)*r^3*(-cosc./(db.^2+rm.^2)+cosd./(db.^2+rn.^2))]
nsA=1-psA./psA2
nsB=1-psB./psB2
figure
plot(x,nsA)
hold on
plot(x,nsB)
title("低阻球體上方連剖視極化率曲線")
legend("nsA","psB")
xlabel("X"),ylabel("ns")
%中間梯度裝置
p1=10
p2=8
u=p2/p1
ps1=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2)./(h^2+x.^2).^5/2]
p3=p2/(1-n)
u=p3/p1
ps2=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2)./(h^2+x.^2).^5/2]
ns=1-ps1./ps2
figure
plot(x,ns)
title("中間梯度視極化率主刨面圖")
xlabel("X"),ylabel("ns")
hold on
figure
for y=-5:1:5
u=p2/p1
ps1=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2+y^2)./(h^2+x.^2+y^2).^5/2]
p3=p2/(1-n)
u=p3/p1
ps2=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*x.^2+y^2)./(h^2+x.^2+y^2).^5/2]
ns=1-ps1./ps2
plot(x,ns)
hold on
end
grid on;
title('隨刨面位置變化')
xlabel("X"),ylabel("ns")
figure
ps=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*X.^2+Y.^2)./(h^2+X.^2+Y.^2).^5/2]
p3=p2/(1-n)
u=p3/p1
ps2=p1.*[1+2*(u-1)./(2*u+1).*r^3*(h^2-2*X.^2+Y.^2)./(h^2+Y.^2+X.^2).^5/2]
ns=1-ps1./ps2
pcolor(X,Y,ns);shading interp
colorbar
title('中間梯度等值線圖')
%對稱四極
figure
for x=0:1:50
u=p2/p1
AB=2:4:100
MN=1
AM=(AB-MN)/2
BN=(AB-MN)/2
AN=AM+MN
BM=BN+MN
K=AM.*AN/MN
da=((x-0.5-AM).^2+h^2).^(1/2)
db=((x+0.5+BN).^2+h^2).^(1/2)
rm=((x-0.5).^2+h^2)^(1/2)
rn=((x+0.5).^2+h^2)^(1/2)
cosa=(da.^2+rm^2-AM.^2)./(2*da.*rm)
cosb=(da.^2+rn^2-AN.^2)./(2*da.*rn)
psA1=p1*[1+2*(u-1)./(2*u+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM.^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN.^2)./(db.*rn.*2)
psB1=p1*[1+2*(u-1)/(2*u+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
p3=p2/(1-n)
u2=p3/p1
psA2=p1*[1+2*(u2-1)./(2*u2+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
psB2=p1*[1+2*(u2-1)/(2*u2+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
nsA=1-psA1./psA2
nsB=1-psB1./psB2
ns=(nsA+nsB)/2
plot(log10(AB/2),ns)
hold on
end
grid on
title("對稱四級視極化率")
xlabel("AB/2"),ylabel("ns")
figure
x=-50:1:50
ab=10:4:100;
[AB,X]=meshgrid(ab,x)
MN=1
AM=(AB-MN)/2
BN=(AB-MN)/2
AN=AM+MN
BM=BN+MN
K=AM.*AN/MN
da=((X-0.5-AM).^2+h^2).^(1/2)
db=((X+0.5+BN).^2+h^2).^(1/2)
rm=((X-0.5).^2+h^2).^(1/2)
rn=((X+0.5).^2+h^2).^(1/2)
cosa=(da.^2+rm.^2-AM.^2)./(2*da.*rm)
cosb=(da.^2+rn.^2-AN.^2)./(2*da.*rn)
psA1=p1*[1+2*(u-1)./(2*u+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
cosc=(db.^2+rm.^2-BM.^2)./(db.*rm.*2)
cosd=(db.^2+rn.^2-BN.^2)./(db.*rn.*2)
psB1=p1*[1+2*(u-1)/(2*u+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
p3=p2/(1-n)
u=p3/p1
psA2=p1*[1+2*(u-1)./(2*u+1)*r.^3*(K.*cosa./(da.^2+rm.^2)-K.*cosb./(da.^2+rn.^2))]
psB2=p1*[1+2*(u-1)/(2*u+1)*r^3*(-K.*cosc./(db.^2+rm.^2)+K.*cosd./(db.^2+rn.^2))]
nsA=1-psA1./psA2
nsB=1-psB1./psB2
ns2=(nsA+nsB)/2
pcolor(X,AB,ns2); shading interp
set(gca,'YDir','reverse');
xlabel('X(m)'),ylabel('AB/2(m)');
colorbar
title("擬斷面圖")
%這是激發極化法
用不同的資料運行一下。
uj5u.com熱心網友回復:
謝謝兄弟的代碼,我回去好好研究一下
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/52604.html
標籤:其他
