动力学第三章
空气动力学第三章
（3.13）
γ /( γ −1)
（3.14）
⎡ ⎤ ⎥ γ +1 c ⎢ = ⎢ ⎥ c* ⎢ 2(1 + γ − 1 M 2 ) ⎥ ⎣ ⎦ 2
1/2
（3.15）
考虑能量方程：
V = 2c p (T0 − T ) = 2γ R (T0 − T ) γ −1
& m G * = ( ) max A
R (1 + γ − 1 M 2 )(γ +1)/[2(γ −1)] 2 & p γ 2 (γ +1)/(γ −1) m = *= 0 ( ) A T0 R r + 1
γ
M
A G 1 2 γ − 1 2 ( γ +1)/[ 2( γ −1)] M )] = = [( )(1 + * A G M γ +1 2
γ − 1 *2 M γ +1
马赫数和临界马赫数的关系曲线如图3.6所示：
当M<1时,M*<1; 当M=1时,M*=1; 当M>1时,M*>1; 当M趋近无穷时;
M* = r +1 r −1
• 3.4 由马赫数表示的质量流流率
& m G = = ρV A
ρ = p / RT
c = γ RT
V γ G = p( ) c RT
V2 = M2 γ RT
T0 γ −1 2 = 1+ M T 2
（3.4 ）
cp =
γR γ −1
公式（3.4）实用于绝热流动和等熵流动。
对于完全气体的等熵流动，其压力和密度与温度的关系 为： p0 T0 γ /(γ −1) ρ0 T0 1/(γ −1) =( ) =( ) T ρ p T 将上述公式与（3.4）结合起来，可以得到压力和密度由 马赫数来表示的关系式如下：
化学动力学
第三章化学动力学(总15页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--第三章 化学动力学3-1．在1 100 K 时，3NH (g)在金属钨丝上发生分解。
实验测定，在不同的3NH (g)的初始压力0p 下所对应的半衰期12t ，获得下列数据 0/Pa p ×104 ×104 ×104 12/min t 试用计算的方法，计算该反应的级数和速率系数。
解： 根据实验数据，反应物3NH (g)的初始压力不断下降，相应的半衰期也不断下降，说明半衰期与反应物的起始浓度（或压力）成正比，这是零级反应的特征，所以基本可以确定是零级反应。
用半衰期法来求反应的级数，根据半衰期法的计算公式12121,121,2n t a t a -⎛⎫= ⎪⎝⎭即 ()12,112,221ln /1ln(/)t t n a a =+把实验数据分别代入，计算得()()12,112,2440,20,1ln /ln 7.6/3.7110ln(/)ln(1.710/3.510)t t n p p --=+=+≈⨯⨯ 同理，用后面两个实验数据计算，得 ()ln 3.7/1.710ln(0.75/1.7)n =+≈所以，该反应为零级反应。
利用零级反应的积分式，计算速率系数。
正规的计算方法应该是分别用3组实验数据，计算得3个速率系数，然后取平均值。
这里只列出用第一组实验数据计算的结果，即010022p at k k == 4310012 3.510Pa 2.310 Pa min 227.6 minp k t -⨯===⨯⋅⨯3-2．某人工放射性元素，能放出α粒子，其半衰期为15 min 。
若该试样有80%被分解，计算所需的时间解：放射性元素的蜕变，符合一级反应的特征。
对于一级反应，已知半衰期的数值，就能得到速率系数的值，因为一级反应的半衰期是与反应物浓度无关的常数。
然后再根据一级反应的定积分式，计算分解80%所需的时间。
第三章 化学反应动力学的计算
(3.15) 式中函数，是变元的函数。若函数族在区间上是确定且可微的，当时， 满足关系式：
则称是微分方程组(3.15)的解。 在真实的化学反应体系中，总能满足上述要求，因此一定存在数值 解，具体的解是用计算机寻找满足初始条件的数值解。 给定的初值 是已知常数。 为了书写方便，一阶微分方程组(3.15)式使用微向量表示，即 初值。 现以 （3.16） 为例说明微分方程组的R-K算法。初始条件： R-K公式为：
9.93522×10-3 1.40291×10-5 5.07574×10-5
9.87084×10-3 1.46652×10-5 1.14494×10-4
9.80689×10-5 1.46078×10-5
0.1×10-7 0.333×10-6 0.356×10-6
0 0.47×10-7 0.44×10-7
开式子： Eular法只取了前二项而忽略了高次项，所以产生了误差。
3． Runge-Kutta方法 .1 常微分方程的Runge-Kutta方法 Runge-Kutta方法是建立在泰勒公式基础上的一种方法。通常采用 的是四
阶R－K公式，即考虑了泰勒公式中四次项，而Eular公式只取了一次 项。故R-K公式比Eular公式有了很大的改进。R-K方法在求解范围大、 精度要求主的情况下是一种比较好的方法，并且计算工作量不算太大， 所以在化学化工中应用颇多。 四阶的R-K公式为：
有已斜率的小线段，即可得方向场的略图(图3-2)
从方向场各点的略图可以推出微分方程的原函数图形。例如画出微
分方程的方向场略图，其解为：在平面上（除原点外）的若干个点，画
第三章自由基聚合(第7周)动力学
R p = kp
☺
fkd kt
1/2
[ I ]1/2[M]
解释：体系粘度随转化率提高，链段重排受到阻碍，
活性末端甚至可能被包埋，双基终止困难，终止速率常 数 kt 显著降低；但此时，体系粘度还不足以严重妨碍单 体扩散，增长速率常数 kp 变化不大,活性链寿命延长， 因此自动加速显著，分子量也同时迅速增加。 转化率继续升高后，粘度大到妨碍单体活动的程 度，增长反应也受扩散控制，此时 kp 开始变小,当综 1/2 合值 kp kt 减小时, 聚合速率开始逐渐降低,直到不 能再继续聚合为止。
18
ln2 0.693 kd = = = 4.375×10-6 s -1 t1 / 2 44 ×3600
fk d 1 2 Rp = k p ( ) c(I)1 2 c(M) kt
1×4.37×10- 1 2 -3 1/ 2 -6 Rp = 145×( ) × (4.0 × 10 ) × 0 . 2 = 4 . 58 × 10 (mol/L• s) 7.0×107
体系体积随聚合反应进行而收缩。实验证明，当一定量单体聚合 时，体系体积收缩与单体转化率成正比，所以测定不同聚合时间 体积，可计算聚合速率。
5
3.7 聚合速率
转化率C(%)与聚合时体积收缩率△V/V0成线性关系：
V 1 X C% = V0 * K
式中，△V为体积收缩值（即聚合物体积与单体体积差）； V0为原始体积值； K为体积变化率。
16
聚合速率对单体浓度呈一级反应是单体自由基形成速率
很快、对引发速率无甚影响的结果。如果初级自由基与
单体的引发反应较慢，可与引发剂分解速率相比拟。
Ri =2fkd[ I ] [M]
Rp = kp fkd kt
第三章_动力学方程的三种基本形式
为计算虚功，可将系统上的力集中到某几个刚体上，如集中到 为计算虚功，可将系统上的力集中到某几个刚体上，如集中到O1O3曲柄上 。
应用力学研究所 李永强 第7页
§3.1 虚功形式的动力学方程－动力学普遍方程
集中后曲柄上的力为：常力偶矩 对它的摩擦力矩为M 集中后曲柄上的力为：常力偶矩M2，轮O1、O2、O3对它的摩擦力矩为 1、M2、M3
2
12
12
δϕ1
M g2
13
1
1
2
g1
x g2
x g3
3
M2' 右转 轴承O3： ϕ 3 ≡ 0 轴承
& 左转， ϕ1 左转，故M13'右转 右转
M 13 = M 13′ = M 1
& & ϕ3r = ϕ31 左转
（3）加惯性力 ） && 左转， && && && 左转， && ϕ1 = ϕ 左转， ϕ 2 = 2ϕ1 左转， ϕ 3 = 0 曲柄O 简化中心为O 曲柄 1O3：简化中心为 1
r r r Fi + N i + Fgi = 0
r r r r & ∆Pi = ( Fi + N i + Fgi ) ⋅ ∆ri = 0
& 在此瞬时，相应的位形上给第 个质点虚速度 r 在此瞬时，相应的位形上给第i个质点虚速度 ∆ri ，第i个质点的虚功率 个质点的虚功率
对于系统可得： 对于系统可得：
& && & ∆PC = ( −mg sin θ − RgC ) ∆vC + ( -M gC ) ∆ϕC = - ( mgr sin θ + 1.5mr 2ϕ ) ∆ϕ
第三章 化学动力学基础
1 d pB υ ν B dt
（4）单位为mol· L-1· s-1
7
影响化学反应速率的内因是：反应物的本性。 影响化学反应速率的外因是：反应的各种条件， 主要包括反应物浓度、温度和催化剂
8
3.2浓度对反应速率的影响——速率方程式
3.2.1 化学反应速率方程 对于一般反应：aA+bByY+zZ 反应速率与反应物浓度间的定量关系为： =kc (A) c (B) 称为化学反应的速率定律或反应的速率方程式
1 k υ α β c Ac B
单位由反应级数而定: 零级反应k的单位为mol· L -1 · s-1, 一级反应k的单位为s-1， 二级反应k的单位为mol -1· L ·s-1 。 它是表征化学反应速率相对大小的物理量。大小 11 与浓度无关，但与温度有关。
2.反应级数 、——c (A)、 c(B)的指数，称为反应级数。 一般有， ≠a、 ≠b 如果=1，表示该反应对A物质为一级反应。 =2 表示该反应对B物质是二级反应。 +——反应总级数。 ★注意： 反应级数由实验确定。可以是零、正整数、分 数和负数。
35
活化络合物：指运动着的两种(或多种)反应物 分子逐渐接近并落入对方的影响范围之内而 形成的处于反应物与产物之间的一种结合状 态。例如下列反应中 NO+ O3─→O2+NO2 NO+ O3─→[O—NOO—O] ─→ O2+NO2 反应物转化为生成物的过程中，分子构型发生 连续变化，生成了活化络合物，它所处的状态 称为过渡状态。
15
3.2.3 确定化学反应速率方程的方法 ------初始速率法
其基本要点为： ①将反应物按不同的组成配制一系列混合物 ②先只改变一种反应物A的浓度，保持其它反 应物浓度不变 ③反应在某一温度下开始进行，记录一定时间间隔内 A的浓度变化，作出图cA-t，确定t=0时的瞬时速率。 若能得到至少两个不同cA条件下瞬时速率，就可确 定A的反应级数。 ④同样的方法，确定其它反应物的反应级数。 这种由反应物初始浓度的变化确定反应速率和 速率方程式的方法，称为初始速率法。
自由度机械系统动力学
1. 解析法
d
t t0 Je 0 Me()
（3.4.6）
若
Me()ab
则
再求出其 反函数
t
t0
Je b
ln ab ab0
f (t)
（3.4.7）
若
d
tt0Je 0abc2
演讲完毕，感谢观 看
（3.4.8）
一、等效力和等效力矩 二、等效质量和等效转动惯量
等效力学模型
等效原则： 等效构件具有的动能＝各构件动能之和
M e
n j 1
m
j
vSj v
2
J
j
j
v
2
J e
n j 1
m
j
vSj
2
J
j
j
2
（3.3.3）
等效质量和等效转动惯量与传动比有关， 而与机械驱动构件的真实速度无关
2W()
Je()
（3.4.3）
若
是以表达式
给出，且为可积函数时，
（3.4.3）可得到解析解。
但是
常常是以线
图或表格形式给出，则只
能用数值积分法来求解。
常用的数值积分法有梯形
法和辛普生法。
运动方程式的求解方法
一、等效力矩是位置的函数时运动方程的求解
二、等效转动惯量是常数、等效力矩是角速度的函数时运动方程
单自由度机械系统可以采用等效力学模型来进行研究，即系统的动力学问题转化为一个等效构件的动力学问题来研究，可以 使问题得到简化。
当取作定轴转动的构件作为等效构件时，作用于系统上 的全部外力折算到该构件上得到等效力矩，系统的全部 质量和转动惯量折算到该构件上得到等效转动惯量。
当取作直线运动的构件作为等效构件时，作用于系统上 的全部外力折算到该构件上得到等效力，系统的全部质 量和转动惯量折算到该构件上得到等效质量。
化学反应动力学-第三章-基元反应动力学
论相比较, Ea却是和温度有关的,即 Ea = f (T)
一般来说,当Ec较大且在低温下反应时,Ea=Ec; 当Ec较小且在高温下反应时,只有采用 ln( k 2 ) 1 1 作图才能得较为理想的直线。
上一内容 下一内容
T
T
返回
五、概率因子 碰撞理论假定分子为刚性硬球,主要考虑了硬
假设A、B双分子的反应是：
A + B A B 产物
只有处于活化状态的（A……B）才能进一步反
应，利用Maxwell的速率分布定律、玻尔兹曼
分布律及统计力学，可得出反应的活化能：
E N0 a
上一内容 下一内容
返回
四、碰撞理论的反应速率公式与讨论 依据分子碰撞理论的两个基本假定，反应速率 公式有二：
以上这些问题在分子碰撞理论中将予解释。
上一内容
下一内容
返回
二、分子碰撞理论和碰撞频率 （一）分子碰撞模型 ⑴分子碰撞的弹性刚球模型： 假定分子是刚性的实心球体,分子占有一定体 积,不考虑分子作用力,分子不能压缩。刚球为 光滑表面,碰撞无摩擦阻力,碰撞时切面方向和 对相对速度不产生任何影响,分子的碰撞是弹 性碰撞。
e
E0 / RT
用此式可求出具有能量ε（≥ 0 ）的分子分数。式 0 E0 / N(N0为 中kB为玻尔兹曼常数，显然 0 Avogadro常数)。这表明反应速率常数与能量大 0 于 （即 E0 / N 0 ）的分子数成正比。 可见，只有具有能量大于 0 的反应物分子才能进 行反应。可见，温度的影响表现在活化能因子 - E0 / RT 上，即活化分子的百分数上，由于T是在 e 指数项上，故其影响显著。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第2章function VTB2(m,c,k,x0,v0,tf,w,f0)%单自由度系统得谐迫振动clcwn=sqrt(k/m);z=c/2/m/wn;lan=w/wnwd=wn*sqrt(1-z^2);A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);t=0:tf/1000:tf;phi1=atan2(x0*wd,v0+z*wn*x0);phi2=atan2(2*z*lan,1-lan^2);B=wn^2*f0/k/sqrt((wn^2-w^2)^2+(2*z*wn*w)^2);x=A*exp(-z*wn*t)、*sin(sqrt(1-z^2)*wn*t+phi1)+B*sin(w*t-phi2); plot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间得关系')function VTB1(m,c,k,x0,v0,tf)%VTB1用来计算单自由度有阻尼自由振动系统得响应%VTB1绘出单自由度有阻尼自由振动系统得响应图%m为质量;c为阻尼;k为刚度;x0为初始位移;v0为初始速度;tf为仿真时间%VTB1(zeta,w,x0,v0,tf)绘出单自由度有阻尼自由振动系统得响应图%程序中z为阻尼系数ξ;wn为固有频率ωn;A为振动幅度;phi为初相位θclcwn=sqrt(k/m);z=c/2/m/wn;wd=wn*sqrt(1-z^2);fprintf('固有频率为%、3g、rad/s、\n',wd);fprintf('阻尼系数为%、3g、\n',z);fprintf('有阻尼得固有频率为%、3g、rad/s、\n',wd);t=0:tf/1000:tf;if z<1A=sqrt(((v0+z*wn*x0)^2+(x0*wd)^2)/wd^2);phi=atan2(x0*wd,v0+z*wn*x0)x=A*exp(-z*wn*t)、*sin(wd*t+phi);fprintf('A=%、3g\n',A);elseif z==1a1=x0;a2=v0+wn*x0;fprintf('a1=%、3g\n',a1);fprintf('a2=%、3g\n',a2);x=(a1+a2*t)、*exp(-wn*t);elsea1=(-v0+(-z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);a2=(v0+(z+sqrt(z^2-1))*wn*x0)/2/wn/sqrt(z^2-1);fprintf('a1=%、3g\n',a1);fprintf('a2=%、3g\n',a2);x=exp(-z*wn*t)、*(a1*exp(-wn*sqrt(z^2-1)*t)+a2*exp(wn*sqrt(z^2-1)*t)); endplot(t,x),gridxlabel('时间(s)')ylabel('位移')title('位移与时间得关系')function jzdd%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-1','wt'); %建立主振型文件fid2=fopen('B-1','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=2;M(2,2)=1、5;M(3,3)=1;%输入刚度矩阵K(1,1)=5;K(1,2)=-2;K(2,1)=-2;K(2,2)=3;K(2,3)=-1;K(3,2)=-1;K(3,3)=1%计算特征值与特征向量D=inv(K)*M; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(3); %B(3)为B中得最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(3);A=B/B(3);endf=sqrt(pp)/2/pi %固有频率单位转换为Hzfprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-1','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-1','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)end%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=1;J2=1;J3=2;k2=1100000;k3=1200000;k4=100000;fid=fopen('chuandi','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R;shita3R=shita2R+1/k3*M2R;M3R=shita2R*(-WN^2*J3)+(1+(-WN^2*J3)/k3)*M2R; shita4R=shita3R+1/k4*M3R;if abs(shita4R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;shita4R];%搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita4R);endfid=fopen('chuandi','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')function cdjz2%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=0、5;J2=1;k2=100000;k3=100000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R;if abs(shita3R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第四个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')function zuoye8%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-2','wt'); %建立主振型文件fid2=fopen('B-2','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=1;M(2,2)=1;M(3,3)=1;%输入刚度矩阵K(1,1)=2;K(1,2)=-1;K(2,1)=- 1;K(2,2)=2;K(2,3)=- 1;K(3,2)=- 1;K(3,3)= 1 %计算特征值与特征向量D=inv(K)*M; %原始动力矩阵A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(3); %B(3)为B中得最后一个元素A=B/B(3);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(3);A=B/B(3);endf=sqrt(pp)fprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-2','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-2','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)endfunction zuoye9%矩阵迭代法求系统得三阶固有频率与主阵型clear allclose allfid1=fopen('A-3','wt'); %建立主振型文件fid2=fopen('B-3','wt'); %建立固有频率文件%输入质量矩阵M(1,1)=4;M(2,2)=2;M(3,3)=1;%输入刚度矩阵K(1,1)=4;K(1,2)=-1;K(2,1)=- 1;K(2,2)=2;K(2,3)=- 1;K(3,2)=- 1;K(3,3)= 1 %计算特征值与特征向量D=inv(K)*M; %原始动力矩阵U=inv(K)A=ones(3,1); %初始振型for i=1:3 %计算三阶固有频率与主振型pp0=0;iB=D*A;pp=1、0/B(1); %B(1)为B中得第一个元素A=B/B(1);while abs((pp-pp0)/pp)>1、e-6pp0=pp;B=D*A;pp=1、0/B(1);A=B/B(1);endf=sqrt(pp)fprintf(fid1,'%20、5f',A); %输入主振型数据fprintf(fid2,'%20、5f',f); %输入固有频率数据D=D-A*A'*M/(A'*M*A*pp);endfid1=fopen('A-3','rt'); %打开主振型文件A=fscanf(fid1,'%f',[3,3]) %主振型写成矩阵fid2=fopen('B-3','rt'); %打开固有频率文件f=fscanf(fid2,'%f',[3,1]) %固有频率写成矩阵t=1:3;h1=figure('numbertitle','off','name','0','pos',[50 200 420 420]);bar(t,f(t,1)),xlabel('频率阶级次'),ylabel('Hz'),title('固有频率'),hold on,gridh1=figure('numbertitle','off','name','1','pos',[50 200 420 420]);bar(t,A(t,1)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('1阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','2','pos',[50 200 420 420]);bar(t,A(t,2)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('2阶主振型'),hold on,gridpause(0、1)h1=figure('numbertitle','off','name','3','pos',[50 200 420 420]);bar(t,A(t,3)),xlabel('自由度(质量块)'),ylabel('振型向量'),title('3阶主振型'),hold on,gridpause(0、1)endfunction cdjz2%chuandijuzhen、m; %传递矩阵方法求固有频率clear all,clear closeJ1=0、5;J2=1;k2=10000;k3=10000;fid=fopen('chuandi3','wt'); %建立(打开)速度文件M1L=0;for WN=0:0、01:2000shita1R=1;M1R=-WN^2*J1;shita2R=shita1R+1/k2*M1R;M2R=shita1R*(-WN^2*J2)+(1+(-WN^2*J2)/k2)*M1R; shita3R=shita2R+1/k3*M2R;if abs(shita3R)<0、005WN %搜索到得固有频率(rad/s),并显示shita=[shita1R;shita2R;shita3R;] %搜索到振型,并显示bar(shita),xlabel('对应得质量块'),ylabel('振型向量')pause(1、0)endfprintf(fid,'%30、15f',shita3R);endfid=fopen('chuandi3','rt');x=fscanf(fid,'%f',[1,200001]);t=1:200001;plot(0、01*t,x);grid,xlabel('频率rad/s'),ylabel('第三个质量块得转角(rad/s)'),title('用传递矩阵法求固有频率')第3章function vtb3(m,c,k,x0,v0,tf,w,f0,delt)%用欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件disp、datfor t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度xd=v0+xdd*delt; %计算速度x=x0+xd*delt; %计算位移xfprintf(fid1,'%10、4f',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间得关系')function vtb4(m,c,k,x0,v0,tf,w,f0,delt)%用改进得欧拉法计算单自由度系统谐迫振动响应wn=sqrt(k/m); %计算固有频率fid1=fopen('disp','wt'); %建立一个位移文件disp、datfor t=0:delt:tf; %delt为时间步长xdd=(f0*sin(w*t)-k*x0-c*v0)/m; %计算加速度x3d=(f0*w*cos(w*t)-k*v0-c*xdd)/m;xd=v0+xdd*delt+x3d*delt^2/2; %计算速度x=x0+xd*delt+xdd*delt^2/2; %计算位移xfprintf(fid1,'%10、4f',x0); %向文件中写数据x0=x;v0=xd;tendfid2=fopen('disp','rt'); %打开disp、dat文件n=tf/delt; %disp、dat文件中位移得个数x=fscanf(fid2,'%f',[1,n]); %将disp、dat文件中文艺写成矩阵t=1:n;plot(t,x),gridxlabel('时间(s)')ylabel('位移(s)')title('位移与时间得关系')function vtb5(tf,delt) %用线性加速度法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp5','wt'); %建立一个位移文件dip5、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度md=inv(m+delt/2*c+1/6*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsex=md*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*x dd0)+delt^2/6*f);%计算位移xdd=6/delt^2*(x-(x0+delt*v0+delt^2/3*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度xdd0=xdd;v0=xd;x0=x;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp5','rt'); %打开disp5、dat文件n=tf/delt; %disp5、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp5、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系') figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') function vtb6(tf,delt) %用线纽马克-β法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp6','wt'); %建立一个位移文件dip6、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度bita=1/6;md=inv(m+delt/2*c+bita*delt^2*k);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsexdd=md*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-bita)*delt^2*xdd0)); %计算加速度xd=v0+delt/2*(xdd0+xdd);%计算速度x=x0+delt*v0+delt^2/2*xdd0+bita*delt^3*(xdd-xdd0)/delt; %计算位移v0=xd;x0=x;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp6','rt'); %打开disp6、dat文件n=tf/delt; %disp6、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp6、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系') figure('numbertitle','off','name','自由度-2得位移','pos',[350 160 400 420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系')420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系') function vtb7(tf,delt) %用威尔逊θ法计算三自由度系统谐迫振动响应,tf为仿真时间,delt为仿真时间步长deltclose all;clcfid1=fopen('disp7','wt'); %建立一个位移文件dip6、datm=2*[1 0 0;0 1 0;0 0 1]; %质量矩阵c=1、5*[2 -1 0;-1 2 -1;0 -1 2]; %阻尼矩阵k=50*[2 -1 0;-1 2 -1;0 -1 2]; %刚度矩阵x0=[1 1 1]'; %初始位移v0=[1 1 1]'; %初始速度theta=1、4;md=inv(k+3*c/theta/delt+6/(theta*delt^2)*m);m1=inv(m);[E,F]=eig(m1*k);diag(sqrt(F)); %计算固有频率(rad/s)for t=0:delt:tf; %delt为时间步长f=[2、00*sin(3、754*t) -2、00*cos(2、2*t) 1、00*sin(2、8*t)]';if t==0; xdd0=m1*(f-k*x0-c*v0); %计算初始加速度elsextheta=md*(m*(2*xdd0+6/theta/delt*v0+6/(theta*delt)^2*x0)+c*(theta*delt /2*xdd0+2*v0+3/theta/delt*x0)+f); %计算(t+θdelt)时刻得速度xddtheta=6/(theta*delt)^2*(xtheta-x0)-6/theta/delt*v0-2*xdd0; %计算(t+θdelt)时刻得加速度xdd=(1-1/theta)*xdd0+1/theta*xddtheta; %计算(t+delt)时刻得加速度xd=v0+delt/2*(xdd0+xdd); %计算(t+delt)速度x=x0+delt*v0+delt^2*(2*xdd0+xdd)/6; %计算(t+delt)位移v0=xd;x0=x;xdd0=xdd;fprintf(fid1,'%10、4f',x0);%向文件中写数据t %显示计算时间步长endendfid2=fopen('disp7','rt'); %打开disp6、dat文件n=tf/delt; %disp7、dat文件中位移得个数x=fscanf(fid2,'%f',[3,n]); %将disp7、dat文件中得位移写成矩阵t=1:n;figure('numbertitle','off','name','自由度-1得位移','pos',[450 180 400 420]);plot(t,x(1,t)),grid,xlabel('时间*0、1秒'),title('自由度-1得位移与时间得关系')420]);plot(t,x(2,t)),grid,xlabel('时间*0、1秒'),title('自由度-2得位移与时间得关系') figure('numbertitle','off','name','自由度-3得位移','pos',[250 140 400 420]);plot(t,x(3,t)),grid,xlabel('时间*0、1秒'),title('自由度-3得位移与时间得关系')。