FEM之一维静态热传导

边界条件：

 |--------------------------------------------------------------|

x=0m，0 摄氏度                                                                   x=100000m，q=K*dT/dx=0.014（表示地底下有稳定热源供给）

参数：

% header 
nodes = 5;  %总共nodes 点，nodes -1 个单元
K  = 2;
Ho = 1e-9;
rho = 2700;
Lambda = 20000;
z0   = 0;
zmax = 100000;
qlower = 0.014;

% ********** Below is only code, all variables defined above  ***************
% calc nodal locations, with more towards top  计算节点的位置，对数分布
% this is a bit arcane, and only necessary for this type of problem
% some sort of mesh generation is typically required in more complex
% problems

Z = zeros(1,nodes);
dels = 1/(nodes-1);     % this is the spacing in log10 space
s = 0:dels:1+dels/2;    % make sure 1 is included
Z = (10.^s);
offset = Z(1)-z0;
Z = Z-offset;

z = Z*zmax/9;

X=z;


U=zeros(nodes,1);

element=[];  %单元与节点关系
for i=1:length(X)-1
    element=[element;i,i+1;];
end


  syms x;
  N=[1-x x];
  dN=diff(N);
  dx=diff(X);
 node_number=length(X);
 element_number=length(element);

%单元矩阵
 Ke=[1 -1;  
     -1 1];

A = zeros(node_number,node_number);
%组装到总体矩阵
 for i=1:element_number
     enodes=element(i,:);   
      Ke=1/dx(i)*[1 -1;
     -1 1];
  A(enodes,enodes)=A(enodes,enodes)+Ke; 
 end

% calc forcing functions
q  = rho* Ho * exp(-z/Lambda);  % heat released at the nodal depths

% calc basic forcing vector 计算右边驱动力向量 ，具体见neil的讲义
% note this is the simplest form of the integration, minor improvement can
% be achieved by improving this integration analytically
 for i=2:element_number
  U(i)=q(i)*(dx(i-1)+dx(i))/(2*K); %积分时面积为dx/2 ,而且从2:element_number的每个节点联系了两个形函数所以要相加(dx(i-1)+dx(i))/2
 end

%上下两个点特殊处理，因为只有一段形函数
U(1)     = q(1)*dx(1)/(2*K);
U(nodes) = q(nodes)*dx(nodes-1)/(2*K);

%BC 边界条件，
A(1,:)=0;
A(1,1)=1;
U(1)  = 0;

U(nodes) = U(nodes) + qlower/K; %qlower/K是右边的边界条件

T=A\U;


figure(1);
subplot(2,2,1);
plot(T,z/1000);
axis ij;

%%%%% from Prof.neil  为了对比用
for n = 2:nodes
    L(n-1) = z(n)-z(n-1);     % coud use the matlab function 'diff' instead
end


% calc basic A matrix
AA = zeros(nodes,nodes);
for n= 2:nodes-1
    AA(n,n-1) = -1/L(n-1);
    AA(n,n+1) = -1/L(n);
    AA(n,n)   = -AA(n,n-1)-AA(n,n+1);
end
AA(1,1) = 1/L(1);               AA(1,2)          = -1/L(1);
AA(nodes,nodes) = 1/L(nodes-1); AA(nodes,nodes-1)= -1/L(nodes-1);


% calc forcing functions
q  = rho* Ho * exp(-z/Lambda);  % heat released at the nodal depths


% calc basic forcing vector
% note this is the simplest form of the integration, minor improvement can
% be achieved by improving this integration analytically
f = zeros(nodes,1);
for n = 2:nodes-1
    f(n) = q(n)*(L(n)+L(n-1))/(2*K);
end
f(1)     = q(1)*L(1)/(2*K);
f(nodes) = q(nodes)*L(nodes-1)/(2*K);


%f = [ q(1)*(l1)/(2*K);  % this shows the next step in improving the
%integration
%    rho* Ho * (exp(-0.5*(z(1)+z(2))/Lambda)*l1+exp(-0.5*(z(2)+z(3))/Lambda)*l2)/(2*K); 
%    rho* Ho * (exp(-0.5*(z(2)+z(3))/Lambda)*l2+exp(-0.5*(z(3)+z(4))/Lambda)*l3)/(2*K); 
%    rho* Ho * (exp(-0.5*(z(3)+z(4))/Lambda)*l3+exp(-0.5*(z(4)+z(5))/Lambda)*l4)/(2*K); 
%    rho* Ho * exp(-0.5*(z(4)+z(5))/Lambda)*(l4)/(2*K)];


% add BCs
% T1 already known, so first line of A matrix is just a 1 and 0s
AA(1,:)=0;
AA(1,1)=1;
f(1)  = 0;
AA


% gradient BC added to f(nodes)
f(nodes) = f(nodes) + qlower/K;
f


T = AA\f
subplot(2,2,2);
plot(T,z/1000,'*k-');
axis ij;