如果假定地表的溫度由下述公式決定：

Ts=T0+beta*h;

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

在地表各處溫度時不相同的。這也導致了地表熱流的不同。

下面的程序用來計算地表地形有起伏，上表面溫度由上述公式決定，下表面溫度為1365°的等溫面時，地表熱流值隨著地形的變化。

使用http://www.cnblogs.com/heaventian/archive/2012/11/23/2784488.html文中的程序，來進行穩態溫度場的計算。

?

首先，代碼coord3.m進行網格劃分，單元，節點的編號與連接。

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

使用下面的程序，可以查看網格剖分情況：

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

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

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

結果如下圖：

?

使用u_d.m程序，將地表溫度設置為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程序確定橫向熱流邊界條件。本文橫向熱流邊界條件設為絕熱邊界條件。

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程序確定泊松方程右側，即溫度載荷。由于不考慮內生熱，無熱流，無對流換熱，因此，右側定義為0.

即ΔT+qv/λ=0;

其中，qv為內生熱，λ為熱導率。

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確定單元剛度矩陣：

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，利用有限元方程，計算溫度場：

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')

計算得到的溫度場如下圖

?

地殼范圍內更細致的溫度場圖如下：

計算地表的熱流值：

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，這個很低。

徑向平均溫度如下圖：

為一直線，這和理論預測一致。

這一溫度和熱流都比真實的地殼的溫度要低，原因是沒有考慮內生熱。

熱流低可能是因為巖石圈200km較厚的原因。因此，也還好。

?

考慮內生熱:

假定體積內生熱qv=質量內生熱*密度=9.6e-10W/kg*3e3kg/m^3=1.1520e-06

地殼的熱導率K=2.5W/K/m.

則qv/K=?4.6080e-07

?

帶入到計算熱載荷項的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) );

此時，地下溫度場如下：

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

看看地殼溫度吧：

這顯然太高了。

再先看下熱流吧：

這個也明顯太高了。

?

?

轉載于:https://www.cnblogs.com/heaventian/archive/2012/11/30/2795599.html

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。