内容发布更新时间 : 2024/6/13 16:17:31星期一 下面是文章的全部内容请认真阅读。

清华大学数学实验报告

实验4：常微分方程数值解

化学工程系 分2 安振华 2012011837

【实验目的】

1. 练习数值微分的计算；

2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法； 3. 通过实例学习用微分方程模型解决简化的实际问题；

4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式及稳定性等概念。

【实验内容】 【实验四：习题5】

射性废物的处理：有一段时间，美国原子能委员会（现为核管理委员会）处理浓缩放射性废物时，把它们装入密封性很好的圆桶中然后扔到水深300ft的大海中。这种做法是否会造成放射性污染，自然引起生物学家及社会各界的关注。原子能委员会一再保证，圆桶非常坚固，绝不会破漏，这种做法是绝对安全的。然而一些工程师们却对此表示怀疑，认为圆桶在和海底相撞时有可能发生破裂。于是双方展开了一场笔墨官司。

究竟谁的意见正确呢？原子能委员会使用的是55gal的圆桶，装满放射性废物时的圆桶重量为527.436 lbf,在海水中受到的浮力为470.327 lbf。此外，下沉时圆桶还要受到海水的阻力，阻力与下沉

清华大学数学实验报告

速度成正比，工程师们做了大量实验，测得其比例系数为0.08lbf s/ft。同时，大量破坏性试验发现当圆桶速度超过40ft/s时，就会因与海底冲撞而发生破裂。

（1） 建立解决上述问题的微分方程模型；

（2）用数值和解析两种方法求解微分方程，并回答谁赢了这场官司。【分析】

为了便于分析，首先，假设圆桶下沉过程无障碍且做直线运动，圆桶质量为m，装满放射性废物时的圆桶重量为M,下沉速度为v，极限速度为u，加速度为a, 重力加速度为g，海水深度为H，圆桶下沉时距海面距离为h.，比例系数为k。

受力分析：圆桶下落过程受重力G、浮力F、阻力f三个力的作用，受力情况如图所示

F f

G

圆桶下落过程分析：由于圆桶刚开始时所受重力和浮力均不变，而所受阻力与速度成正比关系，速度增加时，阻力也随之增加，因此圆桶在此过程做加速度减小的加速直线运动。当阻力f增加到G-F-f=0，即加速度a=0时，阻力和速度同时达到最大，此后速度不变，即圆桶做匀速直线运动。

碰撞过程分析：通过以上分析我们可以猜想圆桶与海底碰撞时的速度

清华大学数学实验报告

有两种可能。第一，圆桶在与海底发生碰撞前就已达到最大速度，也就是说，此时圆桶是以最大速度与海底发生碰撞的；第二，圆桶有可能未达到最大速度就已海底发生碰撞。 【解答】

由已知条件换算得： M=527.436 lbf= 239.241kg， m=470.327lbf=213.337kg， g=9.800m/s^2，

k=0.08 lbf.s/ft=0.119 kg.s/m， u=40ft/s=12.192m/s， H=300ft=91.440m

由题意知：G?Mg,F?mg,f?kgv,v?dh/dt,a?dv/dt 且通过受力分析有：G?F?f?Ma

dv(M?m?kv)g综合上式可得：????(1)dtM

v?dh??????(2)dt

通过以上方程，我们可以计算出圆桶到达海底时的速度v，然后与极限速度u作比较，当v>u时，圆桶与海底发生碰撞后破裂，当v

构建函数：

function dv=fun1(t,v) %构造“时间-速度”函数 M=239.241;m=213.337;g=9.8;k=0.119; %根据题意选取模型参数 dv=(M-m-k*v)*g/M; %表示微分方程（1） 清华大学数学实验报告

MATLAB程序: clf ; ts=linspace(0,1500,15000); %考察圆桶下落25分钟内速度变化情况 v0=0; %假设初速度为0 [t,v]=ode45(@fun1,ts,v0); %解微分方程（1） [t,v]; %输出结果 plot(t,v,'b'); %作出“时间-速度”图象 grid on; %加网格线 xlabel('t-时间'); %横坐标表示时间 ylabel('v-速度'); %纵坐标表示速度 title('圆桶下落速度随时间变化情况'); 运行结果：

圆桶下落速度随时间变化情况250200150v-速度1005000500t-时间10001500

结果分析：

图像显示结果与前面分析大致相同，即圆桶在此过程做加速度减小的加速运动，直到当a=0时，v达到最大。

利用微分方程（1）有：a?dv(M?m?kv)g??0dtM

即 M?m?kv?0

清华大学数学实验报告

此时 vmax?217.681ms??12.192ms?u

这就是前面对圆桶与海底碰撞过程分析猜想的第一种可能情况，即圆桶有可能以最大速度与海底发生碰撞，此时圆桶的速度远远大于极限速度，圆桶与海底发生碰撞后会发生破裂。但由于该计算过程并未考虑海底深度，也就是说圆桶有可能并未达到最大速度就与海底发生碰撞，那么在这种情况下圆通是否能够达到极限速度还无法判断，因此还必须分析速度与圆桶下沉距离的关系，求出圆桶下落到海底时的速度，然后与极限速度作比较，进而得出结论。

MATLAB程序:

clf; ts=linspace(0,300,3000); %考察圆桶下落5分钟内速度变化情况 v0=0; %假设初速度为0 [t,v]=ode45(@fun1,ts,v0); %解微分方程（1） [t,v]; %输出结果 h=cumtrapz(t,v); %利用梯形积分公式求圆桶距海面距离 H=91.44;u=12.192; %根据题意选取模型参数 plot(h,v,'r'); %作出“速度-距离”图象 hold on; plot([0,100],[12.192,12.192],'-g'); %设定水平线--表示极限速度 plot([91.44,91.44],[0,15],'-b'); %设定铅垂线--表示海底 grid on; axis([0,100,0,15]); %设定图像范围 text(91.44,6,'海底'); %标注“海底”的位置 text(35,12.192,'极限速度'); %标注“极限速度”的位置 xlabel('h-圆桶到海面的距离'); ylabel('v-速度'); title('速度随距离变化情况');

运行结果：