《matlab计算方法和技巧小球沉降的规律》由会员分享，可在线阅读，更多相关《matlab计算方法和技巧小球沉降的规律（5页珍藏版）》请在金锄头文库上搜索。
1、小球在水中沉降的规律问题小球的质量为m，受到水的浮力为B，受水的粘性力为R = -Kv，式中K是与水的粘性和小球的半径有关的常量，计算小球在水中竖直沉降随时间变化的速度和深度。数学模型如B6.1图所示，小球受到重力mg，方向向下；浮力B，方向向上；粘性力Kv，方向向上。取向下为正方向，根据牛顿第二定律可列小球的运动方程mg - B - Kv = ma。 (6.1.1)xOmgBRB6.1图显然重力要大于浮力mg B。当速率为零时可得最大加速度aM = g - B/m。当加速度为零时可得极限速率vT = (mg - B)/K。由于a = dv/dt，可得微分方程， (6.1.2)分离变量得，积分
2、，可得，化得，速度为。 (6.1.3)当t时vvT，这就是极限速度。当t = m/K时，小球速率为v = vT(1 - e-1) = 0.6321vT。m/K称为特征时间，用表示；0.6321vT是特征速度。在t 时，小球将以极限速度vT匀速下降。在求小球沉降的深度时，则利用关系v = dx/dt，可得微分方程，积分得，即 。 (6.1.4)当时间t足够长时，小球下沉的深度随时间直线增加x = vT(t - )。 (6.1.5)这是一条直线，也是深度曲线的渐近线。算法方法一：用解析式。令t* = t/，t*就是无量纲的时间，也称为约化时间；令v* = v/vT，v* 是无量纲的速度，称为约化速
3、度。约化速度为v* = 1 exp(-t*)。 (6.1.3*)极限速度与特征时间的乘积是特征长度vT，取无量纲深度或约化深度x* = x/vT，则得x* = t* + exp(-t*) - 1。 (6.1.4*)牛顿运动方程是关于时间的二阶微分方程，当微分方程有解析解时，可深入讨论物体的运动规律，也很容易画出运动曲线。程序zqy6_1ode.m部分如下。%小球沉降的速度和深度(用解析公式)clear %清除变量tm=5; %最大时间dt=0.1; %时间间隔t=0:dt:tm; %时间向量v=1-exp(-t); %求速度x=t+exp(-t)-1; %深度figure %创建图形窗口sub
4、plot(2,1,1) %选子图plot(t,v) %画速度曲线grid on %加网格fs=12; %字体大小title(小球沉降的速度曲线,FontSize,fs)%标题xlabel(时间itt/tau,FontSize,fs) %横坐标ylabel(速度itv/vrm_T,FontSize,fs)%纵坐标subplot(2,1,2) %选子图plot(t,x) %画深度曲线grid on %加网格title(小球沉降的深度曲线,FontSize,fs)%标题xlabel(时间itt/tau,FontSize,fs) %横坐标ylabel(深度itx/vrm_Tittau,FontSize
5、,fs)%纵坐标图示(1)如P6.1之上图所示，速度随时间减速增加，趋于极限速度。(2)如P6.1之下图所示，深度开始时随时间加速增加，最后随时间直线增加。方法二：用一个一阶微分方程的数值解。直接用约化物理量，由(6.1.2)式得，此式可化为，其中， = m/K。利用约化时间t* = t/和约化速度v* = v/vT，可得。 (6.1.2*)这是一个关于约化速度的一阶微分方程，左边是导数，右边是导数的表达式。微分方程的初始条件为：当t* = 0时v* = 0。微分方程的数值解就是速度与时间的关系，深度为 P6.1图。通过无量纲速度对无量纲时间的数值积分可求深度与时间的关系，深度的单位是vT。利
6、用常微分方程数值解指令ode23求速度，利用累积梯形指令cumtrapz求深度。程序zqy6_1ode.m部分如下。%小球沉降的速度和深度(求一阶微分方程的数值解)f=inline(1-v,t,v); %一阶微分方程的内线函数v0=0; %初始条件t,v=ode23(f,t,v0); %求一阶微分方程的数值解subplot(2,1,1) %选子图hold on %保持图像plot(t,v,.r) %画速度曲线x=cumtrapz(v)*dt; %用梯形法求深度subplot(2,1,2) %选子图hold on %保持图像plot(t,x,.r) %画深度曲线说明在求一阶常微分方程的数值解时，
7、由于函数特别简单，用内线函数定义就行了。在定义内线函数时，需要指出自变量和函数的符号。常微分方程数值解的常用指令是ode23，就是2、3阶龙格-库塔解法，虽然精度不是很高，但是求解这类问题的精度已经足够了。ode23指令的第一个输入参数是函数名，第二个输入参数是时间向量，第三个输入参数是初始条件。输出的是时间和速度列向量，这里，输出时间列向量与输入时间行向量的数值相同。cumtrapz指令对速度累积求和之后，需要乘以时间间隔才是深度。方法三：用两个一阶微分方程的数值解。根据速度定义v = dx/dt可得，利用约化深度x* = x/vT，方程可化为。设x(1) = x*，x(2) = dx*/d
8、t*，由上式可得，由(6.1.2)式可得。将一个二阶微分方程，化成两个一阶微分方程，即可求数值解。程序zqy6_1ode.m部分如下。%小球沉降的速度和深度(求二阶微分方程的数值解)x0=0,0; %初始条件t,X=ode45(zqy6_1fun,t,x0); %求二阶微分方程的数值解x=X(:,1); %取深度v=X(:,2); %取速度subplot(2,1,1) %选子图plot(t,v,ko) %画速度曲线subplot(2,1,2) %选子图plot(t,x,ko) %画深度曲线程序在执行时将调用函数zqy6_1fun.m。%小球的约化速度和加速度函数(1)function f=fu
9、n(t,x) %函数过程 f=x(2); %约化速度 1-x(2); %约化加速度说明(1)求二阶常微分方程的数值解时，需要将速度和加速度的表达式用文件表示。第一个输入参数是时间向量，第二个输入参数是初始条件。两个初始条件中，第一个是初始深度，第二个是初始速度。fun是函数名，函数名通常与文件名相同，但也可以不同。将主程序的文本名加上fun就形成函数文件名，这是为了便于查找。函数有两行(写在一行也可以)，第一行是速度的表达式，第二行是加速度的表达式，都是根据公式形成的。(2)求常微分方程数值解的最常用的指令是ode45，就是4、5阶龙格-库塔解法，其精度比ode23指令要高。ode45指令输入
10、参数的排列与ode23指令的相同，时间向量也可用时间范围代替。第一个输出参数是时间列向量，第二个输出参数是深度和速度矩阵，矩阵的第一列表示深度，第二列表示速度。方法四：用微分方程的符号解。由(6.1.2*)式可二阶微分方程。当t = 0时，x = 0，dx/dt = 0，即x* = 0，dx*/dt* = 0，根据初始条件可求微分方程的符号解。求微分方程符号解的指令是dsolve，第一个参数是二阶常微分方程字符串，第二个和第三个参数是初始条件。深度的符号解求出来之后，用符号差分指令(与数值差指令在形式上相同)diff求速度。在求深度和速度时的数值时，只要用subs指令将时间符号用时间向量替换就
11、行了。程序zqy6_1ode.m的计算部分如下。%小球沉降的速度和深度(用微分方程的符号解)sx=dsolve(D2x+Dx-1,x(0)=0,Dx(0)=0);%微分方程的符号解sv=diff(sx); %求速度的符号解x=subs(sx,t,t); %深度v=subs(sv,t,t); %速度subplot(2,1,1) %选子图plot(t,v,mh) %画速度曲线subplot(2,1,2) %选子图plot(t,x,mh) %画深度曲线说明求微分方程的符号解时，微分方程必须有解才能求得，这种方法可验证微分方程是否有解。符号解比数值解运行要慢，因此求微分方程的数值解的方法用得更多。湖南大学物电院 周群益5
