《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 470|回复: 5
收起左侧

[数学] 科学吃西瓜(把西瓜放冰箱里,多长时间才能冰透)

  [复制链接]

该用户从未签到

发表于 2017-8-4 10:17:10 | 显示全部楼层 |阅读模式
2017年8月2日,一个叫做“毕导”的网友发了一篇文章“吃了几十个西瓜后,我终于发明了能避开所有瓜籽的科学吃瓜法!”,全文参见如下网址:
https://mp.weixin.qq.com/s/5yS7O1dxq-DjUV_ZFIpnyg
文中提到了一个有趣的问题:把西瓜放冰箱里,多长时间才能冰透?文中用MATLAB对这个问题做出了解答(文中没有附程序源码),其本质是求解热传导方程:
  1. % 求解冰镇西瓜过程中的温度扩散问题(热传导方程)
  2. %       1   ∂T     1   ∂        ∂
  3. %      ---.---- = ---.----( r2.----)
  4. %       a   ∂t     r2  ∂r       ∂r
  5. %
  6. %   初始条件:
  7. %            t = 0, 0 ≤ r ≤ R, T = T0
  8. %   边值条件:
  9. %            t > 0, r = 0, ∂T/∂r = 0
  10. %            t > 0, r = R, h(T - Tinf) = -k.∂T/∂r
复制代码


% 相关参数
h = 5;        % 西瓜与静止冷空气的对流换热系数 [w/(㎡.K)]
k = 0.48;     % 西瓜的导热系数 [w/(m.K)]
p = 918;      % 西瓜的密度 [kg/m3]
Cp = 3990;    % 热容 [J/(kg.K)]
a = k/(p*Cp); % 西瓜的热扩散系数 [㎡/s]
T0 = 32;      % 初始温度 [℃]
Tinf = 6;     % 终极温度 [℃]

模型假设:
  1. %   1. 假设西瓜的外形是一个半径为R的完美的球体
  2. %   2. 假设西瓜没有瓜皮,没有瓜籽,瓜瓤也是完全均匀的,其密度、导热能力、比热容
  3. %      全是均匀的,一切物性参数都不随温度变化而变化
  4. %   3. 冰箱是一个恒温为Tinf ℃的环境
复制代码


本人对这个问题比较感兴趣,就用MATLAB进行了实现。
首先运行watermelon函数,绘制一个西瓜,如下图所示:
西瓜.jpg

然后运行SolvePde_FrozenWatermelon函数,求解热传导偏微分方程,得到温度随时间和距离变化的三维面图如下:
冰镇西瓜示意图.jpg
两横坐标轴分别代表离瓜心的距离以及冷却时间,纵坐标代表温度。从这张图中可以读出冷却过程中任意时刻、瓜内部任意位置的温度。
运行SolvePde_FrozenWatermelon函数同时得到冰镇西瓜内部温度随时间扩散的动画效果图如下:
冰镇西瓜温度扩散动画.gif

watermelon函数源码如下:
  1. function watermelon
  2. % 绘制一个西瓜
  3. %   August 4  2017, editted by ZhongHua-Xie, Tianjin University of Science and Technology.

  4. % ----绘制西瓜面----
  5. % 西瓜相关参数
  6. a = 10;
  7. b = 10;
  8. c = 15;
  9. % 通过球面坐标定义西瓜外表面网格坐标
  10. Pj = linspace(-pi/2,pi/2,60);
  11. Tj = linspace(0,2*pi,60);
  12. [P,T] = meshgrid(Pj,Tj);
  13. X = a*cos(P).*cos(T);
  14. Y = b*cos(P).*sin(T);
  15. Z = c*sin(P);
  16. id = (Z > 0) & (T > 3*pi/2 & T < 2*pi);
  17. Z(id) = NaN;  % 镂空处理
  18. C = sin(T*10 + rand(size(T)));  % 西瓜外表面颜色索引矩阵
  19. % 自定义颜色矩阵
  20. ColorMat = linspace(0,0.2,64)';
  21. ColorMat = [ColorMat,-3.5*ColorMat + 1,ColorMat];
  22. surf(X,Y,Z,C,'FaceColor','interp','EdgeColor','none')  % 绘制西瓜外表面
  23. colormap(ColorMat)  % 设置西瓜外表面颜色
  24. hold on;
  25. view(31,22)
  26. axis equal;
  27. xlabel('X');
  28. ylabel('Y');
  29. zlabel('Z');

  30. % ----绘制瓜秧(三维螺线)----
  31. Z2 = linspace(0,pi,100);
  32. X2 = Z2.*cos(5*Z2)/2;
  33. Y2 = Z2.*sin(5*Z2)/2;
  34. plot3(X2,Y2,Z2 + 15,'Color',[0.2,0.7,0.2],'LineWidth',3);

  35. % ----绘制西瓜外表面镂空部分对应的内切面----
  36. [P,R] = meshgrid(linspace(0,pi/2,30),linspace(1,0,30));
  37. X = a*cos(P).*R*cos(3*pi/2);
  38. Y = b*cos(P).*R*sin(3*pi/2);
  39. Z = c*sin(P).*R;
  40. surf(X,Y,Z,'FaceColor','r','EdgeColor','none','FaceAlpha',1)
  41. rng(3)
  42. id = randi(900,5,1);
  43. plot3(X(id),Y(id),Z(id),'k.','MarkerSize',12)  % 绘制西瓜籽

  44. X = a*cos(P).*R*cos(2*pi);
  45. Y = b*cos(P).*R*sin(2*pi);
  46. Z = c*sin(P).*R;
  47. surf(X,Y,Z,'FaceColor','r','EdgeColor','none','FaceAlpha',1)
  48. plot3(X(id),Y(id),Z(id),'k.','MarkerSize',12) % 绘制西瓜籽

  49. [T,R] = meshgrid(linspace(3*pi/2,2*pi,30),linspace(1,0,30));
  50. X = a*cos(T).*R;
  51. Y = b*sin(T).*R;
  52. Z = 0*R;
  53. surf(X,Y,Z,'FaceColor','r','EdgeColor','none','FaceAlpha',1)
  54. plot3(X(id),Y(id),Z(id),'k.','MarkerSize',12) % 绘制西瓜籽

  55. % ----绘制坐标轴线----
  56. quiver3(0,0,0,1,0,0,a,'Color',[0.5,0.5,0.2],'ShowArrowHead','off')
  57. quiver3(0,0,0,0,-1,0,b,'Color',[0.5,0.5,0.2],'ShowArrowHead','off')
  58. quiver3(0,0,0,0,0,1,c,'Color',[0.5,0.5,0.2],'ShowArrowHead','off')

  59. caxis([-1,1])
  60. % camlight
  61. % lighting phong
复制代码


SolvePde_FrozenWatermelon函数源码如下:
  1. function SolvePde_FrozenWatermelon
  2. % 把西瓜放冰箱里,多长时间才能冰透?
  3. % 求解冰镇西瓜过程中的温度扩散问题(热传导方程)
  4. %       1   ∂T     1   ∂        ∂
  5. %      ---.---- = ---.----( r2.----)
  6. %       a   ∂t     r2  ∂r       ∂r
  7. %
  8. %   初始条件:
  9. %            t = 0, 0 ≤ r ≤ R, T = T0
  10. %   边值条件:
  11. %            t > 0, r = 0, ∂T/∂r = 0
  12. %            t > 0, r = R, h(T - Tinf) = -k.∂T/∂r
  13. %
  14. %   1. 假设西瓜的外形是一个半径为R的完美的球体
  15. %   2. 假设西瓜没有瓜皮,没有瓜籽,瓜瓤也是完全均匀的,其密度、导热能力、比热容
  16. %      全是均匀的,一切物性参数都不随温度变化而变化
  17. %   3. 冰箱是一个恒温为Tinf ℃的环境
  18. %
  19. %   August 4  2017, editted by ZhongHua-Xie, Tianjin University of Science and Technology.

  20. % 相关参数
  21. h = 5;        % 西瓜与静止冷空气的对流换热系数 [w/(㎡.K)]
  22. k = 0.48;     % 西瓜的导热系数 [w/(m.K)]
  23. p = 918;      % 西瓜的密度 [kg/m3]
  24. Cp = 3990;    % 热容 [J/(kg.K)]
  25. a = k/(p*Cp); % 西瓜的热扩散系数 [㎡/s]
  26. T0 = 32;      % 初始温度 [℃]
  27. Tinf = 6;     % 终极温度 [℃]

  28. m = 2;                                     % 方程中的m参数
  29. r = linspace(0,9,50)/100;                  % 定义距离向量r
  30. t = linspace(0,10,50)*3600;                % 定义时间向量t
  31. Td = pdepe(m,@pdefun,@pdeic,@pdebc,r,t);   % 方程求解

  32. % 结果可视化
  33. figure;
  34. surf(r*100,t/3600,Td);                     % 绘制温度T关于距离r和时间t的三维曲面
  35. colorbar;
  36. zlim([min(Td(:))-2,max(Td(:))+2]);
  37. title('T(r,t)');
  38. xlabel('Distance r');
  39. ylabel('Time t')
  40. view(116,33);
  41. text(8,-0.5,34.5,'开始冷藏','FontSize',18,'Rotation',40)
  42. text(0,5,26,'瓜心','FontSize',18)
  43. text(7,2,10.5,'瓜皮','FontSize',18)

  44. % 绘制西瓜等效切面温度扩散动画
  45. num = size(Td,2);
  46. thetai = linspace(0,2*pi,50);
  47. [R,Theta] = meshgrid(r,thetai);
  48. X = 100*R.*cos(Theta);
  49. Y = 100*R.*sin(Theta);
  50. Tdi = ones(size(thetai'))*Td(1,:);
  51. X_bac = linspace(min(X(:))-1,max(X(:))+1,50);
  52. Y_bac = linspace(min(Y(:))-1,max(Y(:))+1,50);
  53. [X_bac,Y_bac] = meshgrid(X_bac,Y_bac);
  54. figure;
  55. % 绘制背景温度面
  56. surf(X_bac,Y_bac,Tinf*ones(size(X_bac)),'FaceColor','interp','EdgeColor','none');
  57. hold on;
  58. % 绘制零时刻西瓜等效温度切面
  59. hs = surf(X,Y,Tdi,'FaceColor','interp','EdgeColor','none');
  60. axis equal;
  61. axis([min(X_bac(:)),max(X_bac(:)),min(Y_bac(:)),max(Y_bac(:))])
  62. xlabel('X');
  63. ylabel('Y');
  64. caxis([Tinf,max(Td(:))]);
  65. colorbar;
  66. view(2);
  67. ht = title(['冷藏',num2str(t(1)/3600,'%2.1f'),'小时后']);

  68. % filename = '冰镇西瓜温度扩散动画.gif';
  69. % f = getframe(gcf);
  70. % IM = frame2im(f);
  71. % [IM,map] = rgb2ind(IM,256);
  72. % imwrite(IM,map,filename,'gif', 'Loopcount',inf,'DelayTime',0.2);

  73. for i = 2:num
  74.     Tdi = ones(size(thetai'))*Td(i,:);
  75.     set(hs,'ZData',Tdi);
  76.     set(ht,'String',['冷藏',num2str(t(i)/3600,'%2.1f'),'小时后']);
  77.     pause(0.2);
  78.    
  79. %     f = getframe(gcf);
  80. %     IM = frame2im(f);
  81. %     [IM,map] = rgb2ind(IM,256);
  82. %     imwrite(IM,map,filename,'gif','WriteMode','append','DelayTime',0.2);   
  83. end

  84.     function [c,f,s] = pdefun(r,t,T,dT)
  85.         % 偏微分方程函数
  86.         c = 1/a;
  87.         f = dT;
  88.         s = 0;
  89.     end

  90.     function u0 = pdeic(r)
  91.         % 初始条件函数
  92.         u0 = T0;
  93.     end

  94.     function [pa,qa,pb,qb] = pdebc(ra,Ta,rb,Tb,t)
  95.         % 边值条件函数
  96.         pa = 0;
  97.         qa = 1;
  98.         pb = h*(Tb - Tinf);
  99.         qb = k;
  100.     end

  101. end
复制代码

吃了几十个西瓜后,我终于发明了能避开所有瓜籽的科学吃瓜法!.pdf

9.45 MB, 下载次数: 29, 下载积分: 贝壳 -1

该用户从未签到

发表于 2017-8-4 18:25:52 | 显示全部楼层
怎么都没有人回复,这个论坛很好啊,顶顶楼主
回复 支持 反对

使用道具 举报

签到天数: 598 天

[LV.9]以坛为家II

发表于 2017-8-19 09:14:12 | 显示全部楼层
厉害了楼主
回复 支持 反对

使用道具 举报

签到天数: 12 天

[LV.3]偶尔看看II

发表于 2017-9-15 22:21:37 | 显示全部楼层
大家都好久没回来发帖子啦。都太忙了。

顶谢老师。
回复 支持 反对

使用道具 举报

签到天数: 1 天

[LV.1]初来乍到

发表于 2017-10-6 15:30:54 | 显示全部楼层
厉害了,围观看看
回复 支持 反对

使用道具 举报

该用户从未签到

发表于 2017-12-6 15:21:36 | 显示全部楼层
楼主厉害,生活常识都有大理论,理论得到验证。厉害啊。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 陕ICP备08102094号  

GMT+8, 2017-12-14 18:02 , Processed in 1.098168 second(s), 29 queries , Gzip On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表