matlalb与python编程进行动力总成悬置模态计算对比——困惑待解

lijil168

已于 2023-05-07 18:02:58 修改

阅读量519

点赞数 1

分类专栏：数学文章标签： python numpy

于 2022-08-11 12:47:36 首次发布

本文链接：https://blog.csdn.net/lijil168/article/details/126281868

版权

数学专栏收录该内容

120 篇文章 20 订阅

订阅专栏

TOC

1、困惑

以前用python编程算动力总成悬置系统模态时，有时有少量误差，困惑的很，为了搞清原因，将整个理论建模及算法重新研究了一遍，并用matlab编程进行计算对比，奇怪的是，python程序开机重启后计算结果又是完全正确的了。

？关注变量的引用

2 参考模态计算与仿真计算对比

3 python计算结果及代码

在这里插入图片描述

# _*_ coding:UTF-8 _*_
import numpy as np
from numpy import cos
def M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz):
    return np.array([[m,0,0,0,0,0],\
                  [0,m,0,0,0,0],\
                  [0,0,m,0,0,0],\
                  [0,0,0,Ixx,-Ixy,-Ixz],\
                  [0,0,0,-Ixy,Iyy,-Iyz],\
                  [0,0,0,-Ixz,-Iyz,Izz]])
def K_sys_byMatrix(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci):
    #参数均采用数组格式
    n=np.size(Ai)
    ki=np.zeros((n,3,3))
    Ti=np.zeros_like(ki)
    BBi=np.zeros((n,3,6))
    KKi=np.zeros((n,6,6))
    K=np.zeros((6,6))
    for i in range(n):
        ki[i]=np.array([[k_pi[i],0,0],\
                        [0,k_qi[i],0],\
                        [0,0,k_ri[i]]])
        Ti[i]=np.array([[cos(theta_pi[i]),cos(phi_pi[i]),cos(psi_pi[i])],\
                        [cos(theta_qi[i]),cos(phi_qi[i]),cos(psi_qi[i])],\
                        [cos(theta_ri[i]),cos(phi_ri[i]),cos(psi_ri[i])]])
        BBi[i]=np.array([[1,0,0,0,Ci[i],-Bi[i]],\
                        [0,1,0,-Ci[i],0,Ai[i]],\
                        [0,0,1,Bi[i],-Ai[i],0]])
        KKi[i]=BBi[i].T@Ti[i].T@ki[i]@Ti[i]@BBi[i]
    K=np.sum(KKi,axis=0)
    return K
def Energy_DistributionJ(M,Value,Vector_Matrix):
    #3维能量分布：（阶次、Dof、Dof)
    n=M.shape[0]
    KE_klj=np.zeros((n,n,n))    
    for j in range(n):
        Value_j=Value[j] #0维数组：特征值 频率的平方
        Vector_r=Vector_Matrix[:,j] #一维数组：行向量
        Vector_c=Vector_r.reshape(n,1) #改成二维数组：列向量
        KE_klj[j]=0.5*Value_j*M*Vector_c*Vector_r
    return KE_klj
def Energy_percent(KE_klj): 
    n=KE_klj.shape[0]
    #KE_klj 3维能量分布：（阶次、Dof、Dof=页、行、列）
    #将每行Dof的能量合并缩维，得到（各Dof总能量占比,阶）：
    Energy_perDof=np.sum(KE_klj,2)  #按行合并，第3维压缩掉，成2维数组：矩阵（阶次、各Dof能量）
    Energy_AllDof_r=np.sum(Energy_perDof,1) #按行合并，第2维压缩掉，成一维数组：行向量[1阶总能量、2阶总能量、...]
    Energy_AllDof_c=Energy_AllDof_r.reshape(n,1) #改成二维数组：列向量[[1阶总能量],[2阶总能量]、...]
    Energy_percent=100.0*Energy_perDof/Energy_AllDof_c #2维数组：矩阵（阶次、各Dof能量占比）
    Energy_percent=Energy_percent.T ########二维数组转置，得到矩阵（各Dof总能量占比,阶）
    return Energy_percent
###########input###########################################################################
[m,Ixx,Iyy,Izz,Ixy,Ixz,Iyz]=[1224.4,79.385,523.73,491.332,-5.866,-76.403,-4.698]
PG=[0.266249,-0.011964,0.186296]
ABC=[[-0.525,-0.413,0.110],[-0.525,0.413,0.110],[0.675,-0.380,0.151],[0.675,0.380,0.151]]
ABC=np.array(ABC)-np.array(PG)
Ai,Bi,Ci=ABC[:,0],ABC[:,1],ABC[:,2]
theta_pi=[0,0,0,0] #3
theta_qi=[90,90,90,90]
theta_ri=[90,90,90,90] #93
phi_pi=[90,90,90,90]
phi_qi=[0,0,0,0]
phi_ri=[90,90,90,90]
psi_pi=[90,90,90,90] #87
psi_qi=[90,90,90,90]
psi_ri=[0,0,0,0] #3
[theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri]=\
    np.array([theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri])*np.pi/180
k_pi=[1386000,1386000,2745000,2745000]
k_qi=[354000,354000,367500,367500]
k_ri=[1068000,1068000,2025000,2025000]
###########end input###########################################################################
M_sys=M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz)
K_sys=K_sys_byMatrix(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci)
Value, Vector_Matrix = np.linalg.eig(np.linalg.inv(M_sys)@K_sys)
#np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
np.set_printoptions(precision=3, suppress=True)
print(M_sys)
print(K_sys)
print(np.sqrt(Value)/2/np.pi)
print(Vector_Matrix)
#print(Vector_Matrix.T@M_sys@Vector_Matrix)
#print(Vector_Matrix.T@K_sys@Vector_Matrix)
KE_klj=Energy_DistributionJ(M_sys,Value,Vector_Matrix)
print(Energy_percent(KE_klj))

4 matlab计算结果及代码

注意对多维数组，用python时我习惯用页、行、列表示3个维度，用matlab时我习惯用行、列、页表示3个维度。
在这里插入图片描述

%M_sys%
function M=M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz)
M=[m,0,0,0,0,0;
   0,m,0,0,0,0;
   0,0,m,0,0,0;
   0,0,0,Ixx,-Ixy,-Ixz;
   0,0,0,-Ixy,Iyy,-Iyz;
   0,0,0,-Ixz,-Iyz,Izz];
end
%%%-------------------

%K_sys%
function K=K_sys(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci)
    %参数均采用数组格式
    n=length(Ai);
    ki=zeros(3,3,n);
    Ti=zeros(3,3,n);
    BBi=zeros(3,6,n);
    KKi=zeros(6,6,n);
    %K=zeros(6,6);
    for i=1:n
        ki(:,:,i)=[k_pi(i),0,0;
                   0,k_qi(i),0;
                   0,0,k_ri(i)];
        Ti(:,:,i)=[cos(theta_pi(i)),cos(phi_pi(i)),cos(psi_pi(i));
                   cos(theta_qi(i)),cos(phi_qi(i)),cos(psi_qi(i));
                   cos(theta_ri(i)),cos(phi_ri(i)),cos(psi_ri(i))];
        BBi(:,:,i)=[1,0,0,0,Ci(i),-Bi(i);
                    0,1,0,-Ci(i),0,Ai(i);
                    0,0,1,Bi(i),-Ai(i),0];
        KKi(:,:,i)=BBi(:,:,i)'*Ti(:,:,i)'*ki(:,:,i)*Ti(:,:,i)*BBi(:,:,i);
    end
    K=sum(KKi,3);
end
%%%-------------------
function KE_klj=Energy_DistributionJ(M,Value,Vector_Matrix)
    %3维能量分布：（Dof、Dof、阶次)
    n=size(M,1);
    KE_klj=zeros(n,n,n);
    for j=1:n
        Value_j=Value(j); %0维数组：特征值 频率的平方
        Vector_c=Vector_Matrix(:,j);%二维数组：列向量
        Vector_r=Vector_c'; %改成一维数组：行向量
        Vector_ci=repmat(Vector_c,1,n); %列向量展成n行
        Vector_ri=repmat(Vector_r,n,1); %行向量展成n列
        KE_klj(:,:,j)=0.5*Value_j*M.*Vector_ci.*Vector_ri;
    end
end
%%%-------------------
function Energy_perc=Energy_percent(KE_klj)
    n=size(KE_klj,3);
    %KE_klj 3维能量分布：（Dof、Dof、阶次=行、列、页）
    Energy_perDof=sum(KE_klj,2);  %按行合并，第2维为空但是不压缩掉，成3维数组：（各Dof能量和、、阶次）
    Energy_AllDof=sum(Energy_perDof,1); %按列合并，第1、2维为空但是不压缩掉，成3维数组：(各阶能量和、、阶次）
    Energy_perc=zeros(n);
    for j=1:n
        Energy_perc(:,j)=Energy_perDof(:,:,j)./Energy_AllDof(:,:,j); %2维数组：矩阵（各Dof能量占比,阶次）
    end
%%%--------------------------------------------------
%%%--------------------------------------------------
%Mount_egg%
%%%%input:
[m,Ixx,Iyy,Izz,Ixy,Ixz,Iyz]=deal(1224.4,79.385,523.73,491.332,-5.866,-76.403,-4.698);
PG=[0.266249,-0.011964,0.186296];
ABC=[-0.525,-0.413,0.110;-0.525,0.413,0.110;0.675,-0.380,0.151;0.675,0.380,0.151];
ABC=ABC-[PG;PG;PG;PG];
[Ai,Bi,Ci]=deal(ABC(:,1),ABC(:,2),ABC(:,3));
theta_pi=[0,0,0,0]*pi/180; %#3
theta_qi=[90,90,90,90]*pi/180;
theta_ri=[90,90,90,90]*pi/180; %93
phi_pi=[90,90,90,90]*pi/180;
phi_qi=[0,0,0,0]*pi/180;
phi_ri=[90,90,90,90]*pi/180;
psi_pi=[90,90,90,90]*pi/180; %87
psi_qi=[90,90,90,90]*pi/180;
psi_ri=[0,0,0,0]*pi/180; %3
k_pi=[1386000,1386000,2745000,2745000];
k_qi=[354000,354000,367500,367500];
k_ri=[1068000,1068000,2025000,2025000];
%%%-------------------
%%egg:
M=M_sys(m,Ixx,Ixy,Ixz,Iyy,Iyz,Izz);
K=K_sys(k_pi,k_qi,k_ri,theta_pi,phi_pi,psi_pi,theta_qi,phi_qi,psi_qi,theta_ri,phi_ri,psi_ri,Ai,Bi,Ci);
%[eig_vec,eig_val]=eig(inv(M)*K);
[eig_vec,eig_val]=eig(M\K);%比inv计算更快更精确?
wn=sqrt(diag(eig_val));
fn=wn/2/pi;
KE_klj=Energy_DistributionJ(M,diag(eig_val),eig_vec);
Energy_perc=100*Energy_percent(KE_klj);