课题：斯特林发动机机械系统动力学仿真

一、初始设计参数与热力学计算

对于我们设计的β型斯特林发动机，提出了如下的设计目标：

物理参数	数值（单位）
输出功率	0.5W

在我们初步设计的斯特林发动机（模型如下图所示）中，相关的尺寸参数如下：

物理参数	数值（单位）
排气活塞行程h1	42mm
做功活塞行程h2	45mm
相位角α	85°
气缸内径r	10mm
排气活塞半径r0	8mm
气缸内气体压强最小值Pmin	101300Pa（与环境大气压一致）

将设计好的模型导入Ansys软件中进行静态热力学的仿真（如下图所示），可以得到气体温度的状态参数如下：

压缩空间气体温度Tc	439K
膨胀空间气体温度Te	611K

根据史密特理论的相关计算公式，可以编写相应程序，由以上参数为基础计算并绘制P-V图以及单次循环所作功，代码如下：

afa=2*85*pi/360;
theta=0:0.01:pi*2;
tc=439;
te=611;
l1=0.021;
h1=l1*2;
l2=0.0225;
h2=l2*2;
r=0.01;
r0=0.008;
vse=r0*r0*pi*h1;
vsc=r*r*pi*h2;
ve=vse.*(1-cos(theta))./2;
vb=(vse+vsc)./2-sqrt((vse.*vse+vsc.*vsc)./4-vse.*vsc.*cos(afa)./2);
vc=vse.*(1+cos(theta))./2+vsc.*(1-cos(theta-afa))./2-vb;
vr=(r*r-r0*r0)*0.02*pi;
tao=tc/te;
k=vsc/te;
xb=vb/vse;
x=vr/vse;
fai=atan(k.*sin(afa)./(1-tao-k.*cos(afa)));
s=tao+4.*tao.*x./(1+tao)+k+1-2.*xb;
b=sqrt(tao.*tao+2.*k.*(tao-1).*cos(afa)+k.*k-2.*tao+1);
deta=b./s;
pmin=101300;
p=pmin.*(1+deta)./(1-deta.*cos(theta-fai));
plot(rad2deg(theta),p);
xlabel('角度θ（°）');
ylabel('压强P（Pa）');
title('α=85°时θ-P图线');
v=ve+vr+vc;
figure;
plot(rad2deg(theta),v);
xlabel('角度θ（°）');
ylabel('体积V（m^3）');
title('α=85°时θ-V图线');
figure;
plot(v,p);
xlabel('体积V（m^3）');
ylabel('压强P（Pa）');
title('α=85°时P-V图线');
w=10000000*pmin.*vse.*pi.*deta.*(1-tao).*sin(fai).*sqrt(1-deta)./((1+sqrt(1-deta.*deta)).*sqrt(1+deta));
disp(w);

通过运行上述代码，可绘制出如下热力学数据图线，并计算出单次循环做功为**0.0523J**。

因此，若要达到设计目标的0.5W功率要求，需要转速达到rmin=0.5*60/0.0523≈573.6rpm。

此外，以上设计参数所得到的P-V图线偏扁圆形，与常见的P-V图线形状有一定差异，这主要是与设定的初始相位角有关，若将相位角改为45°，则可以得出如下P-V图线，并可计算得出此时对应的单次循环做功为0.0608J，较先前有所提高；但在β型斯特林发动机中，相位角是由相关零件的设计直接确定的，故在后面的仿真中仍然保持相位角α=85°的设定。

二、Adams动力学仿真

在不考虑各类摩擦的情况下，对于基本的曲柄连杆传动机构来说，有如下基本公式：

根据此公式可得到如下代码，绘制转矩变化曲线如下图：

ap=pi*l2*l2;
   fp=ap.*(p-pmin);
   tq=(sin(theta)-cos(theta)).*fp.*r;
   t_qm=w/2/pi;
   t_qm_(1,1:629)=t_qm;
   figure;
   plot(rad2deg(theta),tq,rad2deg(theta),t_qm_);
   xlabel('角度θ（°）');
   ylabel('转矩');
   title('转矩变化曲线');

根据如下仿真步骤，将Fusion360建模软件中建立完成的模型导出为STEP格式，进入Adams仿真软件中进行进一步的动力学仿真。

在第一次仿真时，仅仅将原有模型中设计到传动的部分保留，并将简化后的模型导入仿真，主要反映出两大问题，第一时间进行了修改（以上给的参数均为该次修改后确定的）：

（1）传动部分设计失误，主要表现为各个曲柄的转动不同轴而导致角速度不一致，从而无法达到稳定的压缩与膨胀之间的状态转换，即活塞体系无法完成循环；

（2）连杆设计尺寸不合适，导致部件之间出现穿模问题。

通过修改转动轴连接方式，重新捋清传动循环原理并对不合适的尺寸进行修改（将连接排气活塞的连杆长度缩短），最终实现了在给定初始转速下动力学仿真模型的稳定运转（演示效果见PPT）。

在验证Adams模型与连接建立设置的可行性后，为进行下一步仿真，设置了系统单元变量（推力，对应的力矩，飞轮转动角度与角速度），并将前两者作为输入，后两者作为输出，利用Controls插件导出对应文件以供与MATLAB中Simulink的联合仿真。

三、Simulink联合仿真

打开Adams导出的m文件并运行从而加载对应系统变量，在命令行中输入adams_sys指令即可进入Simulink仿真界面。根据计算与仿真需求，建立Simulink原理图如下：

通过上述分析可知，对于整个发动机来说，我们需要根据不同时刻气缸中密封气体的状态来计算其对活塞的推力，并将其转换为转矩施加在飞轮上，从而形成对整个发动机循环工作的驱动。因此，在Adams建立系统单元时，我们选择气体推力及其等效飞轮转矩（无摩擦情况下）作为整个系统的输入变量，而将飞轮转动的实时角度作为输出量，同时设置飞轮转动的角速度即转速作为最终需要监测的关键输出变量。其中，气体对活塞产生的推力本质上与飞轮转动的实时角度息息相关，也正因此需要将该角度作为Adams仿真系统中的输出变量，以便于在Simulink中回传至模型中进行迭代仿真。而实现角度信息到推力乃至转矩数据的转换，就需要在Simulink中插入S-Function进行实时计算处理。在本原理图中，左侧的plant框体即为我们设置的S-Function函数。

事实上，计算推力乃至转矩关键就要计算气缸内密封气体对于活塞产生的压强，也就是气压值P，而这一值在前面的热力学分析中已有提及，并有相应的公式可以根据实时的飞轮转动角度给出对应气压P的计算。于是，通过对于上述热力学斯特林循环代码的简单改写，结合推力与转矩的转换，可以得到plant函数代码如下：

function [sys,x0,str,ts,simStateCompliance] = plant(t,x,u,flag)
switch flag,
case 0,
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
case 1,
sys=mdlDerivatives(t,x,u);
case 2,
sys=mdlUpdate(t,x,u);
case 3,
sys=mdlOutputs(t,x,u);
case 4,
sys=mdlGetTimeOfNextVarHit(t,x,u);
case 9,
sys=mdlTerminate(t,x,u);
otherwise
DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
end
end
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 2;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1; 
sys = simsizes(sizes);
x0 = [];
str = [];
ts = [0 0];
simStateCompliance = 'UnknownSimState';
end

function sys=mdlDerivatives(t,x,u)
sys = [];
end
function sys=mdlUpdate(t,x,u)
sys = [];
end
function sys=mdlOutputs(t,x,u)
theta=u;
afa=2*85*pi/360;
tc=439;
te=611;
h1=0.042;
h2=0.045;
r=0.01;
r0=0.008;
vse=r0*r0*pi*h1;
vsc=r*r*pi*h2;
vb=(vse+vsc)./2-sqrt((vse.*vse+vsc.*vsc)./4-vse.*vsc.*cos(afa)./2);
vr=(r*r-r0*r0)*0.02*pi;
tao=tc/te;
k=vsc/te;
xb=vb/vse;
x1=vr/vse;
fai=atan(k.*sin(afa)./(1-tao-k.*cos(afa)));
s=tao+4.*tao.*x1./(1+tao)+k+1-2.*xb;
b=sqrt(tao.*tao+2.*k.*(tao-1).*cos(afa)+k.*k-2.*tao+1);
deta=b./s;
pmin=101300;
P=pmin.*(1+deta)./(1-deta.*cos(theta-fai));
Ap=pi*r*r;
Fp=Ap.*(P-pmin);
%对曲柄转矩 N.m
if t<0.1
force=0;
torque=-0.4;
else
Tq=(sin(theta)-cos(theta)).*Fp*r;
force=0;
torque=Tq;
end
sys(1)=force;
sys(2)=torque;
end

function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1; % Example, set the next hit to be one second later.
sys = t + sampleTime;
end
function sys=mdlTerminate(t,x,u)
sys = [];
% end mdlTerminate
end

除此之外，由于两个仿真系统之间的运算单位不统一，在将Adams中角度单位改为rad（弧度），长度单位改为m后，仍然需要对相应的角度与转速数据进行单位上的换算处理（如下图所示），从而可以得到数量级正确的数据图线。

经过以上处理，并将adams_sub单元中的仿真模式（改为interactive）与步长（根据情况改为0.0001或0.0005较为合适）进行对应调整后，可以得到仿真结果的数据波形图如下（仿真时间设置为5s)：

**图线结果分析**：

（1）转速开始突增的原因：初始转矩设置为0.4（负方向），相对于后续整体产生的转矩值都较大，因此在开始时转速在迅速增长后又回落；

（2）可以看到转速在短暂突增后逐渐趋于平稳，但稳定在3000rpm左右，这个转速较快，在现实中可能较难实现，会超过实际材料的承受极限；

（3）但转速在稳定后仍有较大波动，初步判定是由于飞轮质量较轻导致其转动产生的震荡较大，初步观察图线可以发现其在平均值附近最大有将近±300rpm的浮动，需要在后续处理中予以优化；

（4）角度的变化较为规律，且在单位换算调整后能够基本稳定在-180°到180°之间（偶尔有小的浮动），说明仿真运行基本正常；

（5）转矩的变换也与前面动力学分析时绘制的转矩图线基本吻合，随角度变化也呈周期性变化，符合理想情况。

四、迭代改进及其仿真

在上述章节中，我们已经完成了从Adams力学建模仿真到接入Simulink建立联合仿真并得到转速变化图线结果的全过程，但仍然存在一些问题，需要进行进一步修改以达到预期效果并指导发动机对应尺寸的设计修改。

1.初次修改——调节设计参数

在之前的初步仿真中，得出的转速会在短暂突增后稳定在3000rpm左右，这在实际中对于我们设计的结构规模来说是不太可能实现的，同时也远远超出了实际材料的结构强度与承受范围，因此为了将稳定时的转速降低，同时进行了如下两处改动以达到效果：

（1）减小活塞与气缸半径

减小活塞与气缸半径实质上是减小了活塞受气体推力作用的横截面积，同时也减小了气缸中容纳的气体体积，但由于密封前仍与外界空气联通，故初始气压（即Pmin）仍与大气压一致，也因此该改动实际上是降低了相同温度状态下气体的压强，从而导致气体对活塞推力的减小，传动到飞轮上的转矩也随之减小，最终导致了转速的降低。

在本次改动中，我们对活塞与气缸半径均进行了减小处理，并在plant函数中对应的结构参数部分做出了对应的调整，调整后的参数如下表所示：

物理参数	数值（单位）
气缸内径r	5mm
排气活塞半径r0	4mm

该设计参数的调节也同步影响到了单次循环做功的数值，通过重新代入修改后的参数进行热力学参数的计算，得到此参数下单次斯特林循环的做功数值为**0.0327J**。

（2）增大飞轮质量（转动惯量）

在之前的仿真中为便于贴合加工处理将所有部件的材质设置为铝（密度较低，质量轻），但是结果显示转速在大致稳定后仍然出现了一定程度的震荡浮动，这是我们想要避免的。因此，一方面为了降低转速的浮动幅度，另一方面也是为了通过增加转动惯量加大负载从而降低转速，我们希望将飞轮的质量增大；同时，为了尽可能减少我们的修正对于整体建模结构的影响，我们尽量不去改动飞轮本身的尺寸，而是在Adams中对于飞轮部件的质量特性做出调整，通过将其材质由铝转换为密度更大的不锈钢从而实现增加飞轮质量的效果。

综合以上两点调整后，再次进入Simulink中运行仿真，同时由于调整结构设计参数降低了气体对活塞做功施加到飞轮上的等效转矩的整体数值，因此也将初始驱动力矩降低为-0.2（尽可能贴近气体膨胀压缩产生等效力矩的浮动范围）以尽可能降低初始力矩对于最终稳定转速的影响，得到仿真结果如下图（演示视频见PPT）：

**图线结果分析**：

（1）转速在仿真开始阶段的突增消失；

（2）转速在短暂的增加后迅速达到平稳状态，且稳定在1450rpm左右，较先前有了明显的降低，同时结合更改参数后热力学计算得到的单次斯特林循环做功数值为0.0327J，可以得出此时斯特林发动机运行的理论功率为0.7906W，满足初始设定的设计要求；

（3）但转速在稳定后的波动明显减小，观察图线可以发现其在平均值附近仅有±50rpm左右的浮动，运行基本处于稳定状态；

（4）角度的变化较为规律，且在单位换算调整后能够基本稳定在-180°到180°之间（偶尔有小的浮动），说明仿真运行基本正常；

（5）转矩的变换也与前面动力学分析时绘制的转矩图线基本吻合，随角度变化也呈周期性变化，符合理想情况。

2.汇报后的更正——活塞推力的给予方式改变以及传动机构间摩擦的添加

在12月5号课程的汇报中，老师对于我们的仿真提出了一些指导意见，其中最主要的就是提到将气体对活塞推力转化为飞轮转矩的合理性，以及随之带来的传动过程中的摩擦问题。同时老师还提到，不要轻易的去改动原先设计的尺寸参数，事实上在第一阶段的改动中仅仅是在热力学计算以及对应推力与转矩计算中改动了对应的参数，而由于修改模型再重新导入Adams设置连接与参数过于麻烦所以并未对应进行调整，所以我们对于气缸与活塞半径的修改进行了回退处理。

事实上，基于斯特林发动机整体的运作机理考虑，首先需要在预热的同时人为拨动飞轮，本质上与给予初始转矩相对应，因此在t<0.1s的启动阶段，仍然保留了初始力矩值torque的设定，而在后续的循环运作阶段传出实时气体对活塞的推力Fp作为force的值，因此改动plant函数对应部分如下：

pmin=101300;
P=pmin.*(1+deta)./(1-deta.*cos(theta-fai));
Ap=pi*r*r;
Fp=Ap.*(P-pmin);
if t<0.1
force=0;
torque=-0.4;
else
force=Fp;
torque=0;
end
sys(1)=force;
sys(2)=torque;
end

同时为监测对应的推力改变情况，修改Simulink系统原理图如下：

经过对于各种参数的多次调整与运行，我们进行了如下的补充仿真：

（1）给定初始力矩对于转速的影响

我们分别设定t<0.1s时给定的初始转矩Torque=-0.5和-0.05，运行仿真并得出转矩图线：

通过观察上述两图并对比可以发现，事实上初始转矩的大小是决定运行初期转速变化情况的主要因素，具体体现为若初始转矩较大则会出现转速突增现象而降低初始转矩后则该现象消失，同时初始转矩较低时稳定后的转速也有对应降低，但降低幅度有限，大致还是稳定在1200rpm左右；同时在取消转换为传动后力矩这一操作后，转速也较之前有了明显的降低，说明各部件之间连接在动力学实际运作过程中的影响不可忽略，主要是由于各传动部件存在质量且相互的配合之间会产生摩擦。

（2）增加飞轮与支架在旋转时产生的摩擦

综合各因素的调整与尝试，发现在各部件的传动摩擦中只有飞轮与支架之间由于飞轮旋转而产生的摩擦对于整体装置运作的影响相对较大，因此加上原本设置的气缸——做功活塞与做功活塞——排气活塞间摩擦（动摩擦系数均设置为0.5），运行Simulink仿真有如下结果：

观察图线可以发现，增加传动摩擦后，转速在短暂增长（由于初始转矩）后迅速衰减，并在0.34s左右停止转动。

为改变这一现象，计划采取润滑的方式降低动摩擦系数，因此在设置时将动摩擦系数降为0.2，重新运行Simulink仿真，转速图线结果如下：

观察图线可知，当摩擦系数降低后，转速衰减速率也明显降低，但仍有明显衰减且在2.5s左右发生反转现象，且转速逐渐增大。

（3）反转带来的反思——改变初始力矩方向（符号）

根据上述图线，我们发现在0.1s初始力矩作用后正向转速达到最大值，此后却一直衰减直至反向，而飞轮反向旋转后转速却一直增大，且在5s仿真末尾阶段有逐渐平稳的趋势。因此，我们考虑将初始力矩的符号改变（即方向更改）观察仿真效果如下（设置torque=0.05)：

可以发现，在初始力矩反向之后，转速的方向也随之改变，且在初始力矩施加的0.1s后迅速趋于平稳，稳定后转速平均值在**1100rpm左右，仍然大于热力学计算中理论所需的最低转速573.6rpm，满足设计要求**。

该情况下对应的推力force与转角theta图线如下图所示：

五、感悟与总结

通过完成本次动力学仿真，对于斯特林发动机的热力循环过程以及机械转动部分的运作过程有了更加细致和深刻的理解，同时也通过函数的编写与变量的选取对于定量计算中的细节有了更深刻的把握。特别是Adams仿真过程中在建立变量（系统单元）时对于角度和角速度测量的坐标系选取，非常考验对于体系运动过程中细节的把握，究竟是哪个点绕着哪个点转才是正确的角速度，只有正确选取才可能得到正确的结论；单位换算也是仿真过程中的一大难点，不仅要在Adams中对于默认的单位进行修改，还要在MATLAB中编写对应的函数完成转换（Simulink图中的MATLAB Function图标中的函数内容），才可以得到正确的数值结果。除此之外，整个仿真过程中由于涉及到两个软件的联动，对于电脑的性能与配置也有不小的要求（具体体现为把MATLAB和Adams都重装了至少三遍最后还是被迫在别人的电脑上完成，说来也神奇，同样的模型，同样的代码，传到别人电脑上就能跑，我的就不行）。

在完成了整套仿真的流程并得到了正确的数值与图线结果后，还需要针对我们所需要的性能，结合实际情况的可行性，对于对应参数进行修改。经过不断的尝试和摸索，目前有如下的初步结论：

（1）调大飞轮的质量（改变材质，增大材料密度而不改变原本飞轮的尺寸设计），可以使得转速在稳定后的周期性波动幅度明显减小；

（2）为贴近真实情况，应该直接采用气体膨胀压缩对活塞的推力作为输入变量进入Adams仿真系统中进行仿真，转速会比转换为等效转矩作用于飞轮上的更低（因为传动部件本身具有的质量不可忽略，且各部件间存在摩擦）；

（3）传动部分的摩擦损耗主要体现在飞轮在旋转过程中与支架的摩擦力；

（4）如果代入铝材间原有的摩擦系数，发动机无法正常工作，需要通过加润滑油等方式来实现摩擦系数的降低；

（5）初始力矩对应的是发动机在预热阶段用手拨动飞轮产生的驱动力，该转矩会直接影响飞轮在发动机工作初期的转速，并对后续的稳定转速产生一定的影响；因此给定的初始力矩要与后续气体推力对应的等效转矩浮动范围尽可能贴近，以免出现开始时转速突增的现象；

（6）改变气缸与活塞的设计尺寸也能显著改变稳定后的转速，但并不建议这么做，因为需要每次重新导入对应修改后尺寸参数的新模型进行Adams中连接的创立，同时对于加工方面频繁的尺寸修改也不是长久之计。

总而言之，装置的转速与气缸与活塞的尺寸、飞轮质量（转动惯量）、摩擦系数等多方面因素有关。这次的仿真也指导了我们在后续的设计与加工过程中进行进一步的改进，并为发动机整体功率是否满足设计要求从计算上给出了合理的科学依据。