地表地形对地下温度及地表热流的影响
如果假定地表的温度由下述公式决定：
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:Ny for i2=1:Nx k=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-1 for i2=1:Nx-1 k=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 ); return end
使用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; end end % % 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的岩石圈都加了地壳的生热率。
看看地壳温度吧：
这显然太高了。
再先看下热流吧：
这个也明显太高了。
