如果假定地表的溫度由下述公式?jīng)Q定：

Ts=T0+beta*h;

其中，T0是0海拔處地表的平均溫度，h為海拔，beta為空氣的溫度梯度，一般為-6.5K/km。則

在地表各處溫度時(shí)不相同的。這也導(dǎo)致了地表熱流的不同。

下面的程序用來計(jì)算地表地形有起伏，上表面溫度由上述公式?jīng)Q定，下表面溫度為1365°的等溫面時(shí)，地表熱流值隨著地形的變化。

使用http://www.cnblogs.com/heaventian/archive/2012/11/23/2784488.html文中的程序，來進(jìn)行穩(wěn)態(tài)溫度場(chǎng)的計(jì)算。

?

首先，代碼coord3.m進(jìn)行網(wǎng)格劃分，單元，節(jié)點(diǎn)的編號(hào)與連接。

View Code % the coordinate index is from 1~Nx(from left to right) for the bottom line % and then from Nx+1~2*Nx (from left to right) for the next line above % and then next xmin=0;xmax=400e3; Nx=201;Ny=201; x0=linspace(xmin,xmax,Nx); ymax=200e3; ymin=2e3*sin(x0/xmax*2*pi);k=0; for i1=1:Nyfor i2=1:Nxk=k+1;Coord(k,:)=[x0(i2),ymin(i2)+(ymax-ymin(i2))*(i1-1)/(Ny-1)]; end end save coordinates.dat Coord -ascii% the element index k=0; elements3=zeros((Nx-1)*(Ny-1)*2,3); for i1=1:Ny-1for i2=1:Nx-1k=k+1;ijm1=i2+(i1-1)*Nx;ijm2=i2+1+(i1-1)*Nx;ijm3=i2+1+i1*Nx;elements3(k,:)=[ijm1,ijm2,ijm3];k=k+1;ijm1=i2+1+i1*Nx;ijm2=i2+i1*Nx;ijm3=i2+(i1-1)*Nx;elements3(k,:)=[ijm1,ijm2,ijm3];end end %elements3=delaunay(Coord(:,1),Coord(:,2)); save elements3.dat elements3 -ascii% The direchlet boundary condition (index of the two end nodes for each boundary line) boundary=zeros(2*(Nx-1),2); temp1=1:Nx-1; temp2=2:Nx; temp3=1:Nx-1; boundary(temp3',:)=[temp1',temp2']; temp1=Ny*Nx:-1:Ny*Nx-Nx+2; temp2=temp1-1; temp3=temp3+Nx-1; boundary(temp3',:)=[temp1',temp2']; save dirichlet.dat boundary -ascii% The Neuemman boundary condition (index of the two end nodes for each boundary line) boundary=zeros(2*(Ny-1),2); temp1=Nx:Nx:Nx*(Ny-1); temp2=temp1+Nx; temp3=1:Nx-1; boundary(temp3',:)=[temp1',temp2']; temp1=1:Nx:Nx*(Ny-2)+1; temp2=temp1+Nx; temp3=temp3+Nx-1; boundary(temp3',:)=[temp1',temp2']; save neumann.dat boundary -ascii

使用下面的程序，可以查看網(wǎng)格剖分情況：

triplot(elements3,Coord(:,1)*1e-3,Coord(:,2)*1e-3)

為了更清楚觀看，采用橫向，縱向均為21個(gè)網(wǎng)格，并且地表地形起伏達(dá)到20km模型（實(shí)際中，由于要減小誤差，網(wǎng)格劃分為201*201，這樣縱向1km一個(gè)網(wǎng)格）。

xmin=0;xmax=400e3;
Nx=31;Ny=31;
x0=linspace(xmin,xmax,Nx);
ymax=200e3;
ymin=20e3*sin(x0/xmax*2*pi);

結(jié)果如下圖：

?

使用u_d.m程序，將地表溫度設(shè)置為6.5*y=-6.5*h

View Code function value = u_d ( u)%*****************************************************************************80%%% U_D evaluates the Dirichlet boundary conditions.%%% The user must supply the appropriate routine for a given problem%%% Parameters:%% Input, real U(N,M), contains the M-dimensional coordinates of N points.%% Output, VALUE(N), contains the value of the Dirichlet boundary% condition at each point.%value = zeros ( size ( u, 1 ), 1 );value(end-length(value)/2+1:end)=1350; temp1=length(value)/2-1;value(1:length(value)/2)=6.5e-3*2e3*sin((0:temp1)/temp1*2*pi);

?

使用g.m程序確定橫向熱流邊界條件。本文橫向熱流邊界條件設(shè)為絕熱邊界條件。

View Code function value = g ( u )%*****************************************************************************80%%% G evaluates the outward normal values assigned at Neumann boundary conditions.%%% This routine must be changed by the user to reflect a particular problem.%% Parameters:%% Input, real U(N,M), contains the M-dimensional coordinates of N points.%% Output, VALUE(N), contains the value of outward normal at each point% where a Neumann boundary condition is applied.%value = zeros ( size ( u, 1 ), 1 );returnend

?

使用f.m程序確定泊松方程右側(cè)，即溫度載荷。由于不考慮內(nèi)生熱，無熱流，無對(duì)流換熱，因此，右側(cè)定義為0.

即ΔT+qv/λ=0;

其中，qv為內(nèi)生熱，λ為熱導(dǎo)率。

View Code function value = f ( u )%*****************************************************************************80%%% F evaluates the right hand side of Laplace's equation. NOtice that F is qv/K instead of qv.%%% This routine must be changed by the user to reflect a particular problem.%%% Parameters:%% Input, real U(N,M), contains the M-dimensional coordinates of N points.%% Output, VALUE(N), contains the value of the right hand side of Laplace's% equation at each of the points.%n = size ( u, 1 );value(1:n) =0;% ones(n,1);%2.0 * pi * pi * sin ( pi * u(1:n,1) ) .* sin ( pi * u(1:n,2) );

下面的stima3.m確定單元?jiǎng)偠染仃?#xff1a;

View Code function M = stima3 ( vertices )%*****************************************************************************80 % %% STIMA3 determines the local stiffness matrix for a triangular element. % % Discussion: % % Although this routine is intended for 2D usage, the same formulas % work for tetrahedral elements in 3D. The spatial dimension intended % is determined implicitly, from the spatial dimension of the vertices. % % Parameters: % % Input, real VERTICES(1:(D+1),1:D), contains the D-dimensional % coordinates of the vertices. % % Output, real M(1:(D+1),1:(D+1)), the local stiffness matrix % for the element. %d = size ( vertices, 2 );D_eta = [ ones(1,d+1); vertices' ] \ [ zeros(1,d); eye(d) ]; M = det ( [ ones(1,d+1); vertices' ] ) * D_eta * D_eta' / prod ( 1:d );

下面為主程序Heat_conduction_steady.m，利用有限元方程，計(jì)算溫度場(chǎng)：

View Code %*****************************************************************************80 % %% Aapplies the finite element method to steady heat conduction equation, % that is Poission's equation % % % The user supplies datafiles that specify the geometry of the region % and its arrangement into triangular or quadrilateral elements, and % the location and type of the boundary conditions, which can be any % mixture of Neumann and Dirichlet. % % The unknown state variable U(x,y) is assumed to satisfy % Poisson's equation (steady heat conduction equation): % -Uxx(x,y) - Uyy(x,y) = F(x,y) in Omega % with Dirichlet boundary conditions % U(x,y) = U_D(x,y) on Gamma_D % and Neumann boundary conditions on the outward normal derivative: % Un(x,y) = G(x,y) on Gamma_N % If Gamma designates the boundary of the region Omega, % then we presume that % Gamma = Gamma_D + Gamma_N % F(x,y) is qv/K,here qv means volumetric heat generation rate % but the user is free to determine which boundary conditions to % apply. % % The code uses piecewise linear basis functions for triangular elements, % and piecewise isoparametric bilinear basis functions for quadrilateral % elements. % % %clear% % Read the nodal coordinate data file. %load coordinates.dat; % % Read the triangular element data file. %load elements3.dat; % % Read the quadrilateral element data file. % % load elements4.dat; % % Read the Neumann boundary condition data file. % I THINK the purpose of the EVAL command is to create an empty NEUMANN array % if no Neumann file is found. %eval ( 'load neumann.dat;', 'neumann=[];' ); % % Read the Dirichlet boundary condition data file. %load dirichlet.dat;A = sparse ( size(coordinates,1), size(coordinates,1) );b = sparse ( size(coordinates,1), 1 ); % % Assembly. %%{for j = 1 : size(elements3,1)A(elements3(j,:),elements3(j,:)) = A(elements3(j,:),elements3(j,:)) ...+ stima3(coordinates(elements3(j,:),:));end %%} %{for j = 1 : size(elements4,1)A(elements4(j,:),elements4(j,:)) = A(elements4(j,:),elements4(j,:)) ...+ stima4(coordinates(elements4(j,:),:));end %} % Volume Forces. % % from the center of each element to Nodes % Notice that the result of f here means qv/K instead of qv %%{for j = 1 : size(elements3,1)b(elements3(j,:)) = b(elements3(j,:)) ...+ det( [1,1,1; coordinates(elements3(j,:),:)'] ) * ...f(sum(coordinates(elements3(j,:),:))/3)/6;end %%} %{for j = 1 : size(elements4,1)b(elements4(j,:)) = b(elements4(j,:)) ...+ det([1,1,1; coordinates(elements4(j,1:3),:)'] ) * ...f(sum(coordinates(elements4(j,:),:))/4)/4;end %} % Neumann conditions. %if ( ~isempty(neumann) )for j = 1 : size(neumann,1)b(neumann(j,:)) = b(neumann(j,:)) + ...norm(coordinates(neumann(j,1),:) - coordinates(neumann(j,2),:)) * ...g(sum(coordinates(neumann(j,:),:))/2)/2;endend % % Determine which nodes are associated with Dirichlet conditions. % Assign the corresponding entries of U, and adjust the right hand side. %u = sparse ( size(coordinates,1), 1 );BoundNodes = unique ( dirichlet );u(BoundNodes) = u_d ( coordinates(BoundNodes,:));b = b - A * u; % % Compute the solution by solving A * U = B for the remaining unknown values of U. %FreeNodes = setdiff ( 1:size(coordinates,1), BoundNodes );u(FreeNodes) = A(FreeNodes,FreeNodes) \ b(FreeNodes); % % Graphic representation. % % show ( elements3, elements4, coordinates, full ( u ) );figure %%{trisurf ( elements3, coordinates(:,1)*1e-3, coordinates(:,2)*1e-3, full ( u )); %%} %{ trisurf ( elements4, coordinates(:,1), coordinates(:,2), u') %} shading interp xlabel('x')

計(jì)算得到的溫度場(chǎng)如下圖

?

地殼范圍內(nèi)更細(xì)致的溫度場(chǎng)圖如下：

計(jì)算地表的熱流值：

Nx=201,
Ny=201;
Nxy=Nx*Ny;
x0=linspace(0,400e3,201);
temp1=(1:Nx)';
temp2=temp1+Nx;
hs=2.5*(u(temp1)-u(temp2))./(coordinates(temp1,2)-coordinates(temp2,2));
figure,plot(x0*1e-3,hs*1e3)

平均只有16.8mw，這個(gè)很低。

徑向平均溫度如下圖：

為一直線，這和理論預(yù)測(cè)一致。

這一溫度和熱流都比真實(shí)的地殼的溫度要低，原因是沒有考慮內(nèi)生熱。

熱流低可能是因?yàn)閹r石圈200km較厚的原因。因此，也還好。

?

考慮內(nèi)生熱:

假定體積內(nèi)生熱qv=質(zhì)量內(nèi)生熱*密度=9.6e-10W/kg*3e3kg/m^3=1.1520e-06

地殼的熱導(dǎo)率K=2.5W/K/m.

則qv/K=?4.6080e-07

?

帶入到計(jì)算熱載荷項(xiàng)的f.m，將其修改如下：

View Code function value = f ( u )%*****************************************************************************80%%% F evaluates the right hand side of Laplace's equation. NOtice that F is qv/K instead of qv.%%% This routine must be changed by the user to reflect a particular problem.%%% Parameters:%% Input, real U(N,M), contains the M-dimensional coordinates of N points.%% Output, VALUE(N), contains the value of the right hand side of Laplace's% equation at each of the points.%n = size ( u, 1 );value(1:n) =4.6080e-7;% ones(n,1);%2.0 * pi * pi * sin ( pi * u(1:n,1) ) .* sin ( pi * u(1:n,2) );

此時(shí)，地下溫度場(chǎng)如下：

只所以，這樣子，原因在于生熱率是地殼的生熱率，而這里卻是將整個(gè)200km的巖石圈都加了地殼的生熱率。

看看地殼溫度吧：

這顯然太高了。

再先看下熱流吧：

這個(gè)也明顯太高了。

?

?

轉(zhuǎn)載于:https://www.cnblogs.com/heaventian/archive/2012/11/30/2795599.html

總結(jié)

以上是生活随笔為你收集整理的地表地形对地下温度及地表热流的影响的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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