一维PDE传热的问题大家一起学习吧-111.zip

一维PDE传热的问题大家一起学习吧-111.zip 本帖最后由 remie 于 2013-3-11 23:49 编辑     最近忙一个课题，是传热传质的一个pde,现在遇到问题了，写个帖子请教请教大家，顺便把这几天对pde的理解和大家分享分享。     本人上课学pde，也只是学了些皮毛，原本是想自己写个差分方程来的，后来发现强大的matlab 似乎是可以解全部的一维PDE 的问题的，这是因为matlab 的帮助系统里面把 一维的pde问题的表示得很宽泛，如下： 1.png 形如2-2的方程，貌似都能解的。对应于doc pdepe  里面的例子（math大哥提示某个会员看的），本人的c、f、s描述的m文件为：function [c,f,s] = pdex1pde c=1;% x:mm       t:s G=1.46; rhov=0.712; R=8.314; muv=12.22*10^; rhol=954;mul=265*^;A=12400;B=610/999;mmax=100/G-66.7; f=*G*permeav11*rhov*R* 273.15)/muv ...    <font color="red">permeal11*rhol/mul*A*B*-20)))^/-20)))^</font>)*DuDx; s=0;复制代码f就是本人那个非常恶心的函数，对应着方程2-2来看。其中涉及3个m文件，permeav11.m 、permeal11.m和 t10585.m，待会一起给出。 要这3个m文件是表示pde里面的变参数，tmd都是2维的变参数，我本想用拟合的方法把这参数用函数表示的，但是拟合的效果并不好，变参数的图如下： 9 空间完整的温度图像.jpg 6 permeability vapor.jpg 8 渗透性插值.jpg 试了下，f里面是可以用m文件来表达的，所以也就这样用来。m文件的调用里面的输入有一些倍数关系，如t10585，这个t/3600，只是表征物理单位不一样，其他有些地方也是同样的含义。 回来说这些m文件，其实是基于数据的插值，二维函数用的自然是interp2插值了。interp2插值，有正确的输入一般都是有输出的，这个可以参考相关的书籍和matlab的帮助系统。m.文件比较多，待会一起给出。 说到这里该把之前有些东西交代下，我有个startup.m，是在启动matlab时就做一些设置的，内容如下%% Set grids on all axis. set; set; set; %% Set the display format format compact; format short; %% And finally put down anything ye have brought forth clear all; %% project worked on open T105_85_file.m open pdex1.m open pdex1pde.m open pdex1ic.m open pdex1bc.m global darvmm darvxx Dar_v ttt2 XXX2 TTT darlmm darlxx Dar_l load myproject.mat commandwindow复制代码前面几行都是系统设置的，后面才是与这次课题相关的一些内容。load了一个mat文件，也就是我的实验数据，并且申明了全局变量，我昨天就在论坛里面问了个全局变量的问题的，这些开机就自动运行的代码其实主要是为调用变参数的m文件服务的。 mat.的文件如下mat文件好像不是论坛上传的允许文件，所以我就打包了下。。： 新建文件夹.zip 对照doc pdepe里面的描述，解pde的主程序m文件如下：function pdex1 clear; load myproject.mat m = 0; x = linspace; t = linspace; sol = pdepe; % Extract the first solution component as u. u = sol; % A surface plot is often a good way to study a solution. figure;surf hold on; plot3 title xlabel ylabel') 复制代码可以看到里面是有些变量如mmmm,tttt的，这些都存在mat里面，是有startup.m自动调入matlab的workspace的。 边界条件m文件：function [pl,ql,pr,qr] = pdex1bc %pl = ul-0.045 ; %ql = 0; %pr = ur-0.045; %qr = 0; %下面的边界条件对应的x范围是2.78~47.22 pl = ul- 17.39*exp) ; % R^2=0.9916 ql = 0; pr = ur- 0.003091*exp) ; %R^2=0.9959 qr = 0;复制代码带%的内容可以不看，是我给自己坐的标记。 initial value（matlab里面应该叫initial condition吧？？）的m文件，这个比较简单function u0 = pdex1ic u0 = 1.4;复制代码接下来是3个变参数的m文件，语法什么的很接近的。 ------------------------------------------------------------------------function T_fun=t10585 global ttt2 XXX2 TTT T_fun = interp2; 复制代码---------------------------------------------------------------------------function perm_fun=permeav11 global Dar_v darvmm darvxx perm_fun = interp2;复制代码-----------------------------------------------------------------------function perm_fun=permeal11 global Dar_l darlmm darlxx perm_fun = interp2;复制代码------------------------------------------------------------------------------- 要给的程序基本就是这样的了。但是这样是无法解出来的，出错的信息是在调用permeal11.m时有复数的输入。就终止了运算。 我自己是想，在解一些隐式房产时可能出现复数解，然后再调入m文件就出问题了，不知有没有人能分享分享自己的观点。 因为遇到了问题，今天本人已经修改了一天了。把描述微分方程的那个函数，也就是pdex1pde.m里面红色的字去掉是可以解得，解得结果如下图： 10 比较理想的结果1.jpg 10 比较理想的结果2.jpg 这两个图是比较理想的结果，*是实验数据。（当然这只是我的pde的一部分，不完全）。 以上就是我所遇到的问题和一些实用的经验，希望对读者有启发，更希望有人能指出我的不足、错误的地方。 为方便大家学习和测试，我把要用到的所有的m. .mat文件打包，111.zip是可以得到结果的，222.zip就是本人有问题的程序，希望大家多多指点。。 有学习交流的同志可以加本人qq:794145960  注明论坛交流  本人在校大三学生，商业活动、违法活动勿扰。 希望大家多多指点，感激不尽！！ 111.zip 222.zip