clear;
syms xi theta t real;
% % 参数赋值
L = 0.1; % 圆柱壳长度 (m)
h = 5e-4; % 厚度 (m)
R = 0.05; % 半径 (m)
rho = 2600; % 密度 (kg/m³)
zeta = 0.0005; % 阻尼比
mu = 0.3; % 泊松比
E = 7e10; % 弹性模量 (Pa)
Omega_bar = 0; % 无量纲转速
L_r = L / R; % 无量纲长度
H_r = h / R; % 无量纲厚度
f_bar =5.2e-9;
%参数2
R=0.05;
L_r = 5; % 无量纲长度
H_r = 0.002; % 无量纲厚度
% % 参数3（算例2参数 - 来自文献数据）
% L3 = 0.256; % 圆柱壳长度 (m)
% R3 = 0.16; % 平均半径 (m)
% h3 = 0.0025; % 厚度 (m)
% E3 = 110e9; % 弹性模量 (Pa) 110 GPa
% mu3 = 0.31; % 泊松比
% rho3 = 4480; % 密度 (kg/m³)
% Omega_rpm = 20000; % 转速 (r/min)
% Omega3 = Omega_rpm * 2*pi/60; % 转换为 rad/s
%矩阵无量纲化需在能量方程统一除以rho * h * R * L * h/(gamma^2)
% 计算 Q11, Q22, Q12, G
Q11 = E / (1 - mu^2);
Q22 = E / (1 - mu^2);
Q12 = (mu * E) / (1 - mu^2);
G = E / (2 * (1 + mu));
gamma = R*sqrt(rho*(1-mu^2)/E);
% omega_d_bar= gamma*;
% tau= ;
% 高斯基函数参数
Sx = 10; % 高斯中心点数
bls = 2; % 每侧边界扩展层数（仅右侧有意义）
az = 1; % 物理长度
cxz = az / Sx; % 基函数带宽（物理坐标）
acz_end = az + bls * cxz; % 扩展后的区间终点（右侧）
acz_start = 0-bls * cxz; % 起始点固定为0（不允许负轴）
NumX = Sx + 2*bls + 1; % 总点数 = 原始点数 + 右侧扩展层数 +1
qxs = linspace(acz_start, acz_end, NumX); % 物理坐标中心点
FX = exp(-(xi - qxs).^2 / (2 * cxz^2))'; % 无量纲基函数
% overlap = exp(-(qxs(2)-qxs(1))^2/(2*cxz^2)); % 若<0.3则带宽不足
%三角函数
n = 6; % 周向波数
FY = [cos(n*theta),sin(n*theta)]'; % 周向模态函数
NX = length(FX);
NY = length(FY);
% 计算积分矩阵 H_i
[F00x, F01x,F02x,F11x,F22x,F00y,F01y,F02y,F11y,F22y] = Integral(FX, FY);
[H1, H2, H3, H4, H5, H6, H7, H8, H9, H10, H11, H12, H13, H14, H15, H16, H17, H18] = ...
ComputeAllIntegralMatrices(F00x, F01x, F02x, F11x, F22x, F00y, F01y, F02y, F11y, F22y);
% --- 构造无量纲矩阵 ---
% (1) 质量矩阵 M_bar
M_bar = blkdiag(H1, H7, H13);
% (4) 二次刚度矩阵 K_bar
K11 = (1 / L_r^2) * H2 + (Omega_bar^2 + (1 - mu)/2) * H3;
K12 = (1 / L_r) * (mu * H4 + (1 - mu)/2 * H5);
K13 = (mu / L_r) * H6;
K22 = (1 + Omega_bar^2) * H10 + ((1 - mu) / (2 * L_r^2)) * H8;
K23 = H11 + Omega_bar^2 * H11 - Omega_bar^2 * H12;
K33 = H13 + (H_r^2 / (12 * L_r^4)) * H14 + ...
(mu * H_r^2 / (12 * L_r^2)) * H15 + (H_r^2 / 12) * H16 + ...
((1 - mu) * H_r^2 / (6 * L_r^2)) * H17 + Omega_bar^2 * H18;
K_bar = [K11, K12, K13;
K12', K22, K23;
K13', K23', K33];
% 定义人工弹簧的刚度值
Asp = 1e11;
% 定义边界条件类型
PDBC = 3;
if PDBC == 1
% % % % 固支边界条件 (CC)
ku1 = Asp; kv1 = Asp; kw1 = Asp;
pu1 = 0; pv1 = 0; pw1 = Asp;
ku2 = Asp; kv2 = Asp; kw2 = Asp;
pu2 = 0; pv2 = 0; pw2 = Asp;
elseif PDBC == 2
% % % 一端固支一端自由边界条件 (CF)
ku1 = Asp; kv1 = Asp; kw1 = Asp;
pu1 = 0; pv1 = 0; pw1 = Asp;
ku2 = 0; kv2 = 0; kw2 = 0;
pu2 = 0; pv2 = 0; pw2 = 0;
elseif PDBC == 3
% % % 简支边界条件 (SS)
ku1 = 0; kv1 = Asp; kw1 = Asp;
pu1 = 0; pv1 = 0; pw1 = 0;
ku2 = 0; kv2 = Asp; kw2 = Asp;
pu2 = 0; pv2 = 0; pw2 = 0;
end
dFX_dx = diff(FX,xi);
FX1 = subs(FX,xi,0); FX2 = subs(FX,xi,1);
dFX_dx_a1 = subs(dFX_dx,xi,0); dFX_dx_a2 = subs(dFX_dx,xi,1);
F001x = double(conj(FX1)*FX1');
F002x = double(conj(FX2)*FX2');
F00d1x = double(conj(dFX_dx_a1)*dFX_dx_a1');
F00d2x = double(conj(dFX_dx_a2)*dFX_dx_a2');
keduu = ku1*kron(F001x,F00y) + ku2*kron(F002x,F00y);
kedvv = kv1*kron(F001x,F00y) + kv2*kron(F002x,F00y);
kedww = kw1 *kron(F001x,F00y) + kw2*kron(F002x,F00y) + (1/(L_r^2*R^2))*pw1*kron(F00d1x,F00y) + (1/(L_r^2*R^2))*pw2*kron(F00d2x,F00y);
% Ked_dim = -(mu^2 - 1)*blkdiag(keduu,kedvv,kedww);
Ked_dim = ((1-mu^2)/(E*H_r*L_r))*blkdiag(keduu,kedvv,kedww);
% 合并边界条件
K_global = (K_bar + Ked_dim);
% (3) 陀螺矩阵 G_bar
G_bar = 2 * Omega_bar * [zeros(size(H1)), zeros(size(H1)), zeros(size(H1));
zeros(size(H7)), zeros(size(H7)), H9;
zeros(size(H13)), -H9', zeros(size(H13))];
% %二次项刚度矩阵
% [Kx111y000, Kx100y011, Kx001y110, Kx011y100, ...
% Kx000y111, Kx110y001, Kx011y000, Kx110y000, ...
% Kx000y001, Kx000y011, Kx000y110, Kx101y010] = Integral3D(FX, FY);
% %二次刚度矩阵
% Knon2uu = zeros(NX*NY,NX*NY);
% Knon2uv = zeros(NX*NY,NX*NY);
% Knon2uw = (h*R/(2*L^2))* Q11* Kx111y000 + (h/(2*R))* Q12* Kx100y011+ h/R* G* Kx001y110;
% Knon2vv = zeros(NX*NY,NX*NY);
% Knon2vw = h/(2*L)* Q12* Kx011y100 + h*L/(2*R^2)* Q22* Kx000y111 + h/L* G* Kx110y001;
% Knon2ww = (h*R/(2*L^2))* Q11* Kx111y000 + h/(2*L)* Q12* Kx110y001 + h/(2*L)* Q12* Kx011y000...
% + h/L* Q12* Kx110y000 + h/(2*R)* Q12* Kx001y110 + h*L/(2*R^2)* Q22* Kx000y111...
% + h*L/(2*R^2)* Q22* Kx000y011 + h*L/(R^2)* Q22* Kx000y110...
% + h/R* G* Kx100y011 + h/L* G* Kx101y010;
% K_nonlinear2 = [Knon2uu, Knon2uv, Knon2uw;
% Knon2uv', Knon2vv, Knon2vw;
% Knon2uw', Knon2vw', Knon2ww];
% % K_nonlinear2dimless=h*K_nonlinear2/(rho * h * R * L * h/(gamma^2))/h^2;
% K_nonlinear2dimless=h^2*K_nonlinear2/(rho * h * R * L * h/(gamma^2));
% %三次刚度矩阵
% [Kx1111y0000, Kx0000y1111, Kx1100y0011, Kx0011y1100] = Integral4D(FX, FY);
% Knon3uu = zeros(NX*NY,NX*NY);
% Knon3uv = zeros(NX*NY,NX*NY);
% Knon3uw = zeros(NX*NY,NX*NY);
% Knon3vv = zeros(NX*NY,NX*NY);
% Knon3vw = zeros(NX*NY,NX*NY);
% Knon3ww = (h*R/(2*L^3))* Q11* Kx1111y0000 + (h*L/(2*R^3))* Q22* Kx0000y1111 ...
% + (h/(2*R*L)*Q12 + h/(R*L)*G)*(Kx1100y0011 + Kx0011y1100);
% K_nonlinear3 = [Knon3uu, Knon3uv, Knon3uw;
% Knon3uv', Knon3vv, Knon3vw;
% Knon3uw', Knon3vw', Knon3ww];
% K_nonlinear3dimless = h^3 * K_nonlinear3 / (rho * h * R * L * h/(gamma^2));
% %无量纲外激励
% W=kron(FX,FY);
% F=[zeros(size(W)); zeros(size(W)); W*f_bar*cos(omega_d_bar*tau)/(H_r^2*L_r)];
% --- 合并为全局矩阵 ---
N = size(M_bar, 1);
A = [zeros(N), eye(N);
-M_bar\K_global, -M_bar\G_bar]; % 状态空间矩阵
% 计算特征频率
[Te,eig_vals] = eigs(A,1e2,'smallestabs');
frequencies=imag(diag(eig_vals));
% save('tezhengxiangliang.mat',"frequencies",'Te');
NX = length(FX);
NY = length(FY);
DispDofs=3*NX*NY;
ModeTruncation=1;
VsN10=real(Te(1:DispDofs, ModeTruncation));
AlphXL=FX;
B0ysL=FY;
ux0=AlphXL;
vx0=AlphXL;
wx0=AlphXL;
ZerNx=zeros(size(ux0));
A0xuL=[ux0;ZerNx;ZerNx];
A0xvL=[ZerNx;vx0;ZerNx];
A0xwL=[ZerNx;ZerNx;wx0];
% % 分别形成两个方向对于的向量函数（也即变量分离）
J1mx=ones(size(A0xuL));
J1ny=ones(size(B0ysL));
AlphXLextendUx=kron(A0xuL,J1ny); %%% U在x方向的振型函数向量
AlphXLextendVx=kron(A0xvL,J1ny); %%% V在x方向的振型函数向量
AlphXLextendWx=kron(A0xwL,J1ny); %%% W在x方向的振型函数向量
ModeFU=AlphXLextendUx'*VsN10; %%% U在x方向的振型函数
ModeFV=AlphXLextendVx'*VsN10; %%% V在x方向的振型函数
ModeFW=AlphXLextendWx'*VsN10; %%% W在x方向的振型函数
% 创建符号函数
ModeFUi = matlabFunction(ModeFU);
ModeFVi = matlabFunction(ModeFV);
ModeFWi = matlabFunction(ModeFW);
% 先计算函数值，然后归一化
xa=0;xb=1;
xs=linspace(xa,xb,1e4);
Uy0 = ModeFUi(xs);
Vy0 = ModeFVi(xs);
Wy0 = ModeFWi(xs);
% 找到最大值的位置进行归一化
[maxU, idxU] = max(abs(Uy0));
[maxV, idxV] = max(abs(Vy0));
[maxW, idxW] = max(abs(Wy0));
ModeFU = ModeFU / maxU;
ModeFV = ModeFV / maxV;
ModeFW = ModeFW / maxW;
save('ModeFUVW.mat','ModeFU','ModeFV','ModeFW');
ModeFUi = matlabFunction(ModeFU);
ModeFVi = matlabFunction(ModeFV);
ModeFWi = matlabFunction(ModeFW);
xa=0;xb=1;
xs=linspace(xa,xb,1e3);
Uy1 = ModeFUi(xs);
Vy1 = ModeFVi(xs);
Wy1 = ModeFWi(xs);
plot(xs,Uy1,'r-');
hold on;
plot(xs,Vy1,'g-');
hold on;
plot(xs,Wy1,'b-');
hold off;
legend('U','V','W');
% 指定目标文件夹路径
% folder_path = 'H:\旋转圆柱壳\旋转圆柱壳代码'; % Windows示例
% 拼接完整的文件路径
% full_file_path = fullfile(folder_path, 'ModeFUVW.mat');
% 保存文件到指定位置
% save(full_file_path, "ModeFU", 'ModeFV', 'ModeFW');
clear
syms xi theta t real;
n = 6; % 周向波数
FY = [cos(n*theta),sin(n*theta)]'; % 周向模态函数
load('ModeFUVW.mat');
NX=length(ModeFU);NY=length(FY);
% 参数赋值
L = 0.1; % 圆柱壳长度 (m)
h = 5e-4; % 厚度 (m)
R = 0.05; % 半径 (m)
rho = 2600; % 密度 (kg/m³)
zeta = 0.0005; % 阻尼比
mu = 0.3; % 泊松比
E = 7e10; % 弹性模量 (Pa)
Omega_bar = 0; % 无量纲转速
L_r = L / R; % 无量纲长度
H_r = h / R; % 无量纲厚度
%参数2
R=0.05;
L_r = 5; % 无量纲长度
H_r = 0.002; % 无量纲厚度
%矩阵无量纲化需在能量方程统一除以rho * h * R * L * h/(gamma^2)
% 计算 Q11, Q22, Q12, G
Q11 = E / (1 - mu^2);
Q22 = E / (1 - mu^2);
Q12 = (mu * E) / (1 - mu^2);
G = E / (2 * (1 + mu));
gamma = R*sqrt(rho*(1-mu^2)/E);
% 计算积分矩阵 H_i
[Fuu00x,Fuu11x] = Integraluux(ModeFU);
[Fvv00x,Fvv11x] = Integraluux(ModeFV);
[Fww00x,Fww11x,Fww22x,Fww02x] = Integralwwx(ModeFW);
[Fuv01x,Fuv10x] = Integraluvx(ModeFU,ModeFV);
[Fuw10x] = Integraluwx(ModeFU,ModeFW);
[Fvw00x] = Integralvwx(ModeFV,ModeFW);
[F00y,F01y,F02y,F11y,F22y] = Integraly(FY,0,2*pi);
[H1, H2, H3, H4, H5, H6, H7, H8, H9, H10, H11,H12, H13, H14, H15, H16, H17, H18] = ...
ComputeAllIntegralMatricest_2(Fuu00x, Fuu11x, Fuv10x, Fuv01x, Fuw10x, Fvv00x,...
Fvv11x, Fvw00x, Fww00x, Fww22x, Fww02x, Fww11x, F00y, F01y, F02y, F11y, F22y);
% --- 构造无量纲矩阵 ---
% (1) 质量矩阵 M_bar
M_bar = blkdiag(H1, H7, H13)
% (4) 刚度矩阵 K_bar
K11 = (1 / L_r^2) * H2 + (Omega_bar^2 + (1 - mu)/2) * H3;
K12 = (1 / L_r) * (mu * H4 + (1 - mu)/2 * H5);
K13 = (mu / L_r) * H6;
K22 = (1 + Omega_bar^2) * H10 + ((1 - mu) / (2 * L_r^2)) * H8;
K23 = H11 + Omega_bar^2 * H11 - Omega_bar^2 * H12;
K33 = H13 + (H_r^2 / (12 * L_r^4)) * H14 + ...
(mu * H_r^2 / (12 * L_r^2)) * H15 + ...
(H_r^2 / 12) * H16 + ...
((1 - mu) * H_r^2 / (6 * L_r^2)) * H17 + ...
Omega_bar^2 * H18;
K_bar = [K11, K12, K13;
K12', K22, K23;
K13', K23', K33]
% 合并边界条件
K_global = K_bar;
% (3) 陀螺矩阵 G_bar
G_bar = 2 * Omega_bar * [ zeros(size(H1)), zeros(size(H1)), zeros(size(H1));
zeros(size(H7)), zeros(size(H7)), H9;
zeros(size(H13)), -H9', zeros(size(H13))];
%二次项刚度矩阵
[KxU1W1W1yU0W0W0, KxU1W0W0yU0W1W1, KxU0W0W1yU1W1W0] = Integral3Duw(ModeFU, ModeFW, FY, FY)
[KxV0W1W1yV1W0W0, KxV0W0W0yV1W1W1, KxV1W1W0yV0W0W1] = Integral3Dvw(ModeFV, ModeFW, FY, FY)
[KxW1W1U1yW0W0U0, KxW1W1V0yW0W0V1, KxW0W1W1yW0W0W0,...
KxW1W1W0yW0W0W0, KxW0W0U1yW1W1U0, KxW0W0V0yW1W1V1,...
KxW0W0W0yW0W1W1, KxW0W0W0yW1W1W0, KxW1W0U0yW0W1U1,...
KxW1W0V1yW0W1V0] = Integral3Dww(ModeFU, ModeFV, ModeFW, FY, FY, FY)
% (4) 刚度矩阵 K_bar
Knon2uu = zeros(NX*NY,NX*NY);
Knon2uv = zeros(NX*NY,NX*NY);
Knon2uw = (h*R/(2*L^2))* Q11* KxU1W1W1yU0W0W0 + (h/(2*R))* Q12* KxU1W0W0yU0W1W1+ h/R* G* KxU0W0W1yU1W1W0;
Knon2vv = zeros(NX*NY,NX*NY);
Knon2vw = h/(2*L)* Q12* KxV0W1W1yV1W0W0 + h*L/(2*R^2)* Q22* KxV0W0W0yV1W1W1 + h/L* G* KxV1W1W0yV0W0W1;
Knon2ww = (h*R/(2*L^2))* Q11* KxW1W1U1yW0W0U0 + h/(2*L)* Q12* KxW1W1V0yW0W0V1 + h/(2*L)* Q12* KxW0W1W1yW0W0W0...
+ h/L* Q12* KxW1W1W0yW0W0W0 + h/(2*R)* Q12* KxW0W0U1yW1W1U0 + h*L/(2*R^2)* Q22* KxW0W0V0yW1W1V1...
+ h*L/(2*R^2)* Q22* KxW0W0W0yW0W1W1 + h*L/(R^2)* Q22* KxW0W0W0yW1W1W0...
+ h/R* G* KxW1W0U0yW0W1U1 + h/L* G* KxW1W0V1yW0W1V0;
K_nonlinear2 = [Knon2uu, Knon2uv, Knon2uw;
Knon2uv', Knon2vv, Knon2vw;
Knon2uw', Knon2vw', Knon2ww]
K_nonlinear2dimless=h^2*K_nonlinear2/(rho * h * R * L * h/(gamma^2))
% Knon2ww = h/(2*L)* Q12* KxW0W1W1yW0W0W0...
% + h/L* Q12* KxW1W1W0yW0W0W0 ...
% + h*L/(2*R^2)* Q22* KxW0W0W0yW0W1W1 + h*L/(R^2)* Q22* KxW0W0W0yW1W1W0;
K_nonlinear2dimless = vpa(simplify(expand(K_nonlinear2dimless)), 2); % 8位有效数字
disp('小数系数矩阵:');
K_nonlinear2dimless
%三次刚度矩阵
[Kx1111y0000, Kx0000y1111, Kx1100y0011, Kx0011y1100] = Integral4D(ModeFW, FY)
Knon3uu = zeros(NX*NY,NX*NY);
Knon3uv = zeros(NX*NY,NX*NY);
Knon3uw = zeros(NX*NY,NX*NY);
Knon3vv = zeros(NX*NY,NX*NY);
Knon3vw = zeros(NX*NY,NX*NY);
Knon3ww = (h*R/(2*L^3))* Q11* Kx1111y0000 + (h*L/(2*R^3))* Q22* Kx0000y1111 ...
+ (h/(2*R*L)*Q12 + h/(R*L)*G)*(Kx1100y0011 + Kx0011y1100);
K_nonlinear3 = [Knon3uu, Knon3uv, Knon3uw;
Knon3uv', Knon3vv, Knon3vw;
Knon3uw', Knon3vw', Knon3ww];
K_nonlinear3dimless=h^3*K_nonlinear3/(rho * h * R * L * h/(gamma^2))
K_nonlinear3dimless = vpa(simplify(expand(K_nonlinear3dimless)), 2); % 8位有效数字
disp('小数系数矩阵:');
K_nonlinear3dimless
%阻尼矩阵
c_u = 2 * zeta * rho * h * 0.0231; % u方向
c_bar=sqrt((1-mu^2)/(rho*E))*c_u;
C_bar= 1/H_r*c_bar*blkdiag(H1, H7, H13);
% % %无量纲外激励
W=kron(ModeFW,FY);
% 如果ModeFW和FY是函数句柄
% 计算在xi=0.5, theta=0时的值
W_value = subs(W,{xi,theta},[0.5, 0]);
disp('在xi=0.5, theta=0时的外激励向量:');
disp(W_value);
f_bar =5.2e-9;
F=[zeros(size(W)); zeros(size(W)); f_bar*double(W_value)/(H_r^2*L_r)]
%
% --- 合并为全局矩阵 ---
N = size(M_bar, 1);
A = [zeros(N), eye(N);
-M_bar\K_global, -M_bar\G_bar]; % 状态空间矩阵
% 计算特征频率
[Te,eig_vals] = eigs(A,1e2,'smallestabs');
frequencies=imag(diag(eig_vals));
% save('MKGnonlinear.mat',"M_bar",'K_global','G_bar','K_nonlinear2dimless','K_nonlinear3dimless','F');
% 指定目标文件夹路径
folder_path = 'H:\旋转圆柱壳\旋转圆柱壳代码(无量纲)\ThickShellCODE\ThickShellCODE'; % Windows示例
% 拼接完整的文件路径
full_file_path = fullfile(folder_path, 'MKGnonlinear.mat');
% 保存文件到指定位置
save(full_file_path, "M_bar", 'K_global','C_bar', 'G_bar', 'K_nonlinear2dimless', 'K_nonlinear3dimless', 'F');
这是我的全部代码，质量矩阵和刚度矩阵肯定没问题，因为里面和振型回代求解的频率一致，并且与已有文献结果一致
最新发布
扫一扫
5757
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
