永冻土层热传导问题数学建模

参考文献

[1] 张辰熙,于明.季节冻土层下人工冻土墙温度场实验研究[J].低温建筑

术.2010(6):144-146.

[2]文斌,吴青柏,蒋观利,张鹏.模拟退火优化算法的冻土热传导参数反分析[J].

岩 土力学.2013(8):15-16.

[3] 王建伟,李书民，李振卿.冻土层等效导热系数的探讨[J].防渗技

术.2001(3):36-38.

[4]李述训,吴通化.冻土温度状况研究方法和应用分析[J].冰川冻土.2004(4)：

142-148.

[5] 候彬彬,卢德唐,陈永辉,孔祥言.非达西渗流模型模拟青藏铁路抛石路基冻

土层温度[J].水工渗流研究与应用进展第五届全国水利工程渗透学术研讨会论文集.2006.

[6] 史册.热传导方程有限差分法的MATLAB实现[J].咸阳师范学院学

报.2009(4):24-26.

[7] 张世民.青藏铁路多年冻土路基热-力稳定性数值仿真分析[J].土木建筑与

环境工程.2012(34):3-7.

[8] 穆彦虎,马巍,牛富俊,李国玉,王大雁,刘永智.青藏铁路多年冻土区普通路

基热状况监测分析[J].冰川冻土.2014(2):10-13.

[9] 田立慧，凌贤长,王力娜,张锋.青藏铁路高温多年冻土区列车行驶路基长期

永久变形数值模拟研究[J].地震工程学报.2014(4):21-25.

13

附录

问题一

外界环境温度随时间变化的函数关系式拟合程序 creatFit.m

function [fitresult, gof] = createFit(x, y) %CREATEFIT(X,Y) % Create a fit. %

% Data for 'temperature' fit: % X Input : x % Y Output: y % Output:

% fitresult : a fit object representing the fit. % gof : structure with goodness-of fit info. %

% 另请参阅 FIT, CFIT, SFIT.

% 由 MATLAB 于 17-Aug-2015 19:27:00 自动生成

%% Fit: 'temperature'.

[xData, yData] = prepareCurveData( x, y );

% Set up fittype and options. ft = fittype( 'fourier2' );

opts = fitoptions( 'Method', 'NonlinearLeastSquares' ); opts.Display = 'Off'; opts.Normalize = 'on';

opts.StartPoint = [0 0 0 0 0 1.9316882339819];

% Fit model to data.

[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.

figure( 'Name', 'temperature' ); h = plot( fitresult, xData, yData );

legend( h, 'y vs. x', 'temperature', 'Location', 'NorthEast' ); % Label axes xlabel( 'x' ); ylabel( 'y' ); grid on

14

nihe.m x=1:24;

y=[12 12 13 13 13 13 13 13 12 13 13 15 17 18 19 21 22 22 21 20 16 13 12 12]; [fitresult, gof] = createFit(x, y)

偏微分方程求解程序

function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典显式格式求解抛物型偏微分方程

%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %

%方程：u_t=C*u_xx 0 <= x <= uX,0<= t <= uT %初值条件：u(x,0)=phi(x)

%边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t) %

%输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……

% x -空间变量 % t -时间变量

%输入参数：uX -空间变量x的取值上限 % uT -时间变量t的取值上限

% phi -初值条件，定义为内联函数 % psi1 -边值条件，定义为内联函数 % psi2 -边值条件，定义为内联函数 % M -沿x轴的等分区间数 % N -沿t轴的等分区间数 % C -系数，默认情况下C=1 %

%应用举例： % uX=1;uT=0.2;

% M=15;N=100;C=1;

% phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');

%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C); %

inline('263+15.11+2.216*cos(1.757*t)+3.687*sin(1.757*t)-1.895*cos(2*1.757*t)+1.728*sin(2*1.757*t)') %设置参数C的默认值 if nargin==7 C=1; end

%计算步长

dx=uX/M;%x的步长

15

dt=uT/N;%t的步长

x=(0:M)*dx; t=(0:N)*dt;

r=C*dt/dx/dx;%步长比 r1=1-2*r;

if r > 0.5

disp('r > 0.5,不稳定') end

%计算初值和边值 U=zeros(M+1,N+1); for i=1:M+1

U(i,1)=phi(x(i)); end

for j=1:N+1

U(1,j)=psi1(t(j)); U(M+1,j)=psi2(t(j)); end

%逐层求解 for j=1:N

for i=2:M

U(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); end end

U=U';

%作出图形 mesh(x,t,U);

title('温度分布图像') ylabel('空间变量 x') xlabel('时间变量 t ') zlabel('温度 U')

return;

16