键合图-多体及多物理系统建模利器

1、 键合图的广义状态变量

键合图是多体及多物理系统建模的利器，是数学抽象与形象表示的完美接合。

1.1势是流的原因

1.2流是势的原因

1.3动量是势对时间的累积效应

1.4位移是流对时间的累积效应

1.5功率是势和流的相干效应

1.6动能（磁能）是流对冲量（磁通量）的累积效应

1.7势能（电能）是势对位移（电荷）的累积效应

2、键合图到方框图的转换

2.1键合图元转方框图元

2.2 RC电路键合图转方框图

3、键合图列写系统状态方程

建立简单的二自由度振动机械系统的键合图，并转成方框图，列出系统状态方程，然后对三种建模方法进行仿真计算。

4、AMESIM键合图建模

AMEsim的键合图模型非常好用，封装的十分漂亮。m1和m2的位移、m1的速度见下图：

5、Matlab Simulink或AMESIM方框图建模

AMEsim的方框图模型也很好用，根据键合图转换成方框图，m1的速度如下，和键合图的结果一致。

6、Adams微分方程建模

Adams不仅可以进行3d动力学仿真，也可以进行一维动力学仿真，可以建立传递函数模型、状态空间方程模型和通用的微分方程模型，可以建立显式微分方程，也可以建立隐式微分方程，十分方便。这样，对于一些不太复杂的一维多物理场和3d动力学的耦合模型，我们可以直接在Adams中完成，不必进行与Simulink或者AMEsim的联合仿真。

仿真的结果和前面用键合图和方框图完全一样，m1和m2的质量块位移如下：

7、键合图场

建立结型结构，其它键合图元仅与结型结构相连，储能场按因果关系分成独立储能场（积分）和非独立储能场（微分）。

8、非线性系统键合图及状态方程的推导

8.1、常见非线性系统键合图及状态方程的推导

按前一节键合图场的理论建立系统键合图，除耗能场为线性系统外，其他系统采用非线性系统，并推导出系统状态方程，推导的主要难点是复合向量函数的求导，特别雅克比函数矩阵的计算：1、对向量函数可直接求其雅克比函数矩阵；2、对函数矩阵*向量，需要考虑先乘得到向量函数，再求向量函数的雅克比矩阵，比如[J(x)i1 J(x)i2 …J(x)ip ] [U1 U2 … Up]‘,其雅克比矩阵为
[[dJ(x)i1/dx1 dJ(x)i2/dx1 …dJ(x)ip/dx1 ] [U1 U2 … Up]’ … [dJ(x)i1/dxn dJ(x)i2/dxn …dJ(x)ip/dxn ] [U1 U2 … Up]']

8.2、常见非线性系统键合图及状态方程的推导示例

示例来自《键合图理论及其在系统动力学中的应用》，一个凸轮传动机构，3个质量块，本来是3*2个自由度，但是由于凸轮的约束，一储能场为非独立储能场，所以系统有5个独立自由度，键合图的设计比书上稍简单（书上对第二个弹簧加了一个1-节点，实际可以不用加，只需要将13键的方向换一下即可），书上对Jii(Xi)的推导有点错误，也许是印刷错误，我买的是复印版。不同的场函数矩阵有不同的特点，可以检验推导是否有误，这一点很不错。
推导过程虽然繁琐，但是整个过程逻辑性很强，并不困难，经过推导，对系统的认识有较大的提高。

8.3、AMEsim键合图建模验证

AMEsim可能没有质量块微分子模型，用方框模型代替，计算结果和书上略有差异，不清楚是模型差异还是计算的问题，非线性问题不好算准。
su
凸轮位移：

从动件位移：

凸轮速度：

8.4、Adams建模验证

Adams计算的结果也略有差异，X轴与书上设置的相反，注意设置好初始状态的各零件的位置和弹簧预紧力、初始速度，两外，转盘用等效平动质量代替，保证连接弹簧平动伸缩。

凸轮位移：

从动件位移：

8.5、系统状态方程建模验证

matlab编程建模验证，用matlab的符号运算得到系统状态方程，解的结果和AMEsim一样，证明书上不是印刷错误，是写错了。

图上的绿色线是k2弹簧的压缩变形，即q2,其初值应设为-20,使其和q1一样有初始载荷，而x=0是他们的平衡位置！！

 %-- 2024/1/2 21:55 --%
syms k1 k2 k3 m1 J2 m3 RT r q1 q2 q3 p6 p8 dFt Ft
Jii=[0 0 0 1 0;
0 0 0 1 r;
0 0 0 -cos(0.05*q1) 0;
-1 -1 cos(q1/20) 0 0;
0 -r 0 0 0];
Jid=[0;0;0;cos(q1/20);0];
JiL=zeros(5,1);
JiL(4)=-1;
JiU=zeros(5,1);
JiU(4)=1;
Jdi=[0 0 0 -cos(q1/20) 0];
JdU=[0];
JLi=[0 0 0 1 0];
JLd=0;
JLU=0;
Ad=[p6/20*sin(q1/20)/m1 0 0 0 0];
AU=zeros(1,5);
J_PHIsi=diag([k1 k2 k3 1/m1 1/J2]);
J_PHIsd=m3;
R=RT;
JLL=0;
Md=(Jid+JiL*inv((eye(length(R))-R*JLL))*R*JLd)*J_PHIsd;
Mdf=Md*(Jdi*J_PHIsi+Ad+AU);
Xi=[q1;q2;q3;p6;p8];
PHIsi=J_PHIsi*Xi;
dU=dFt;
U=Ft;
dXi=inv(eye(length(Mdf))-Mdf)*(Jii*PHIsi+Md*JdU*dU+JiU*U+JiL*inv((eye(length(R))-R*JLL))*R*(JLi*PHIsi+JLU*U));
df_sym=subs(dXi,{'q1','q2','q3','p6','p8'},{'x(1)','x(2)','x(3)','x(4)','x(5)'});
df_sym=subs(df_sym,{'m1','J2','m3','k1','k2','k3','r','Ft','RT'},{0.02,50,0.01,15,2,2,10,18,0.15});
n=length(df_sym);
i=1;
ss='[';   %先定义好很重要，否则再循环体中定义时，每一循环ss不累加。
while i<n
ss=strcat(ss,char(df_sym(i)),';');
i=i+1;
end
ss=strcat(ss,char(df_sym(i)));
ss=strcat(ss,']');
f=inline(ss,'t','x');
[t,x]=ode45(f,[0 1],[-20,-20,20*(1+sin(1)),0,0]);
plot(t,[x(:,1),x(:,2),x(:,3)]);

用pretty()写系统状态方程如下：

9、键合图解动力学逆问题

9.1、混合因果关系线性系统源输入计算的统一方程

键合图解动力学逆问题，由系统的响应计算系统的激励，对解决系统控制问题很有用，当要求系统按设定的目标状态进行动作时，需要知道力该如何控制。当系统的载荷不好测量时，也可以根据系统的运动状态对载荷进行估算。

9.2、示例

9.3、Matlab符号计算推导

 %-- 2024/1/9 18:36 --%
syms k1 k2 k3 m1 m2 m3 RT p7 p14 q2 q8 q11 FM vt dp7 dp14 dq2 dq8 dq11 dFM dvt m
Xi=[p7;p14;q2;q8;q11];
dXi=[dp7;dp14;dq2;dq8;dq11];
Fi=diag([1/m1,1/m2,k3,k1,k2]);
Fd=m3;
R=RT;
U1=0;
dU1=0;
U2=[FM;vt];
dU2=[dFM;dvt];

Jii=[0 0 -m -1 1;
0 0 0 0 -1;
m 0 0 0 0;
1 0 0 0 0;
-1 1 0 0 0];
Jid=[-m;0;0;0;0];
JiL=[-1;-1;0;0;0];
JiU1=zeros(5,1);
JiU2=[0 0;
    1 0;
    0 1;
    0 0;
    0 0];
Jdi=[m 0 0 0 0];
JdU1=0;
JdU2=[0 0];
JLi=[-1 1 0 0 0];
JLd=0;
JLL=0;
JLU1=0;
JLU2=[0 0];

IL=eye(length(R));
In=eye(length(Fi));
A1=Jid*Fd*JdU2+JiL*R*inv(IL-JLL*R)*JLd*Fd*JdU2;
A2=Jid*Fd*JdU1+JiL*R*inv(IL-JLL*R)*JLd*Fd*JdU1;
A3=JiL*R*inv(IL-JLL*R)*JLU2+JiU2;
A4=JiL*R*inv(IL-JLL*R)*JLU1+JiU1;
A5=Jii*Fi+JiL*R*inv(IL-JLL*R)*JLi*Fi;
A6=Jid*Fd*Jdi*Fi+JiL*R*inv(IL-JLL*R)*JLd*Fd*Jdi*Fi-In;

df_left_sym=A1*dU2;
df_right_sym=A3*U2+A2*dU1+A4*U1+A5*Xi+A6*dXi;

df_left_val=subs(df_left_sym,{'m1','m2','m3','k1','k2','k3','m','RT'},{0.02,50,0.01,15,2,2,2,0.15});
df_right_val=subs(df_right_sym,{'m1','m2','m3','k1','k2','k3','m','RT'},{0.02,50,0.01,15,2,2,2,0.15});

其方程为：

9.4、AMEsim键合图建模验证

先由18N的力计算系统状态，再根据系统的状态计算力，计算的力=-18N，是因为力的方向定义的不同。