辛方法计算轮轨力

clear; clc;
%% 列车参数
V0=[100250300380]; V1=V0*1000/3600;
Mc=31.6e3;Jc=1.55e6;
Mt=3200;Jt=1752;
Mw=2000;
Kz1=1200e3;Kz2=115e3;
Cz1=25e3;Cz2=12e3;
l1=1.25;l2=9;
Rw=0.43;

Mv=diag（[McJcMtJtMtJtMwMwMwMw]）;

Ccc=diag（[2*Cz22*Cz2*l2^2]）; Cww=diag（[Cz1Cz1Cz1Cz1]）;
Ctt=diag（[Cz2+2*Cz12*Cz1*l1^2Cz2+2*Cz12*Cz1*l1^2]）;
Cct=[-Cz20-Cz20;Cz2*l20-Cz2*l20];
Ctw=[-Cz1-Cz100;Cz1*l1-Cz1*l100;00-Cz1-Cz1;00Cz1*l1-Cz1*l1];
OO=zeros（24）;
Cv=[CccCctOO;Cct‘CttCtw;OO‘Ctw‘Cww];

Kcc=diag（[2*Kz22*Kz2*l2^2]）; Kww=diag（[Kz1Kz1Kz1Kz1]）;
Ktt=diag（[Kz2+2*Kz12*Kz1*l1^2Kz2+2*Kz12*Kz1*l1^2]）;
Kct=[-Kz20-Kz20;Kz2*l20-Kz2*l20];
Ktw=[-Kz1-Kz100;Kz1*l1-Kz1*l100;00-Kz1-Kz1;00Kz1*l1-Kz1*l1];
Kv=[KccKctOO;Kct‘KttKtw;OO‘Ktw‘Kww];

%赫兹