如果假定地表的温度由下述公式决定：
Ts=T0+beta*h;
其中，T0是0海拔处地表的平均温度，h为海拔，beta为空气的温度梯度，一般为-6.5K/km。则
在地表各处温度时不相同的。这也导致了地表热流的不同。
下面的程序用来计算地表地形有起伏，上表面温度由上述公式决定，下表面温度为1365°的等温面时，地表热流值随着地形的变化。
使用/heaventian/archive//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 nextxmin=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)]; endendsave coordinates.dat Coord -ascii% the element indexk=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];endend%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 interpxlabel('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的岩石圈都加了地壳的生热率。
看看地壳温度吧：
这显然太高了。
再先看下热流吧：
这个也明显太高了。
