% meshfree_vibration_analysis.m
% Free vibration analysis of laminated closed conical/cylindrical shells and annular plates
% using TRPIM meshfree strong-form method (FSDT) in MATLAB
% pinv
clear; clc; close all;
%% 1. Parameters
% Geometry
R1 = 1.0; % base radius
L = 4.0; % axial length
h = 0.1; % thickness
phiSpan = 3/2*pi; % circumferential span (full shell)
% Material (orthotropic laminate)
E1 = 150e9; E2 = 10e9; mu12 = 0.25; rho = 1500;
G12 = 5e9; G13 = 5e9; G23 = 6e9; ks = 5/6;
layup = [0,90,0]; % fiber angles [deg]
Nk = numel(layup);
% Discretization
nx = 16; ntheta = 16; % nodes in x and theta directions
dist_func = @(a,b,n) 0.5*( (b-a)*(1 - cos(pi*(0:n-1)/(n-1))) ) + a; % 'gauss_lobatto'
% dist_func = @(a,b,n) (0:n-1).^2 / (n-1)^2 * (b - a) + a; % 平方渐变示例
X=linspace(0,L,nx);
T=linspace(0,phiSpan,ntheta);
% T=linspace(0,phiSpan,ntheta+1); T(end)=[];
[xg, tg] = meshgrid(X,T);
yg = R1 * cos(tg); zg = R1 * sin(tg);
nnodes = [xg(:), yg(:), zg(:), tg(:)];
nodes = [xg(:), tg(:)]; N = size(nodes,1);
% TRIPM settings
rbf_type = 'MQ'; % 'MQ','EXP','TPS'
dc = 2*(L/(nx-1)+phiSpan/(ntheta-1))/2; % characteristic length
q = 1.03; % MQ exponent\meta = 2; % TPS exponent
mp=0; mq=0;
m = mp*mq; % Tchebychev polynomial total order
bc_type='C';
% Support domain: K-nearest neighbors
ns = 25; % support size
idxNN = knnsearch(nnodes(:,1:3),nnodes(:,1:3),'K',ns);
% in=nnodes(idxNN(1,:),:);
% f1=draw_shell(in);
f=draw_shell(nnodes);
hold on
surf(xg,yg,zg, 'FaceAlpha',0.3, 'EdgeColor','none');
% 标记目标点（红色）
target_idx = 1;
scatter3(nnodes(target_idx,1), nnodes(target_idx,2), nnodes(target_idx,3), ...
100, [0 0 0], 'filled', 'MarkerFaceAlpha',0.8)
% 绘制欧氏距离邻居（蓝色线段）
for i = 1:ns
scatter3(nnodes(idxNN(target_idx,i),1), nnodes(idxNN(target_idx,i),2), nnodes(idxNN(target_idx,i),3), ...
100, [1 0.3 0.3], 'filled', 'MarkerFaceAlpha',0.8)
plot3([nnodes(target_idx,1) nnodes(idxNN(target_idx,i),1)], ...
[nnodes(target_idx,2) nnodes(idxNN(target_idx,i),2)], ...
[nnodes(target_idx,3) nnodes(idxNN(target_idx,i),3)], ...
'b-', 'LineWidth',0.8, 'Color',[1 0.6 0])
end
%% 2. Precompute basis matrices (G matrix inverse) for each node
Ginv = cell(N,1);
supp = cell(N,1);
for i = 1:N
Xi = nodes(i,1); Ti = nodes(i,2);
nbr = idxNN(i,:); supp{i} = nbr;
% build R0 and Tm at support nodes
Xs = nodes(nbr,1); Ts = nodes(nbr,2);
% compute pairwise radial basis
R0 = zeros(ns);
for a=1:ns
for b=1:ns
% Tab(a,b)=Ts(a)-Ts(b);
% Tan1=wrapToPi(Tab);
% rab = sqrt((Xs(a)-Xs(b))^2 - R1^2*(2*cos(Ts(a) - Ts(b))-2));
rab = sqrt((Xs(a)-Xs(b))^2 + (R1*wrapToPi(Ts(a)-Ts(b)))^2);
R0(a,b) = radialRB(rab, rbf_type, dc, q);
end
end
% Tchebychev basis at support nodes
xp = 2*(Xs - 0)/L - 1; tp = 2*(Ts - 0)/phiSpan - 1;
Tm = tchebyBasis(xp,tp,mp,mq);
% assemble G
G = [R0, Tm; Tm', zeros(size(Tm,2))];
% fprintf('cond(G) = %g\n', cond(G));
% %给 G 矩阵加正则化,加一点微小阻尼
% epsReg = 1e-8 * max(abs(diag(G)));
% G = G + eye(size(G)) * epsReg;
Ginv{i} = inv(G);
% [U,S,V] = svd(G);
% s = diag(S);
% tol = max(size(G)) * eps(max(s));
% s_inv = s;
% s_inv(s > tol) = 1./s(s > tol);
% Ginv{i} = V * diag(s_inv) * U';
end
%% 3. Collocation: assemble global K and M
ndof = 5; totalDOF = ndof*N;
K = sparse(totalDOF,totalDOF);
M = sparse(totalDOF,totalDOF);
% material stiffness matrices A,B,D and transverse shear
[A,B,D_mat,As,I] = laminateABD(layup,E1,E2,mu12,G12,G13,G23,h,rho,ks);
Dstruct = [A,B,zeros(3,2);B,D_mat,zeros(3,2);zeros(2,6),As];
Phi=cell(N,1);
for i = 1:N
% evaluation at node i
Xi = nodes(i,1); Ti = nodes(i,2);
nbr = supp{i}; nn = numel(nbr);
% compute shape function and derivatives at eval point
Phi{i} = computeShape(nodes(nbr,:), [Xi,Ti], Ginv{i}, rbf_type, dc, q, mp,mq, L, phiSpan, R1);
% local stiffness and mass kernels
[Ki_cell, Mi_cell] = localKM(Phi{i}, A, B, D_mat, As, R1, I(1), I(2), I(3));
% assemble
rows = (i-1)*ndof + (1:ndof);
for k = 1:nn
col = (nbr(k)-1)*ndof + (1:ndof);
K(rows,col) = K(rows,col) + Ki_cell(:,:,k);
M(rows,col) = M(rows,col) + Mi_cell(:,:,k);
end
end
figure;
subplot(1,2,1);
spy(K);
subplot(1,2,2);
spy(M); %可视化稀疏结构，排查是否有整行全为0
%% 4. Apply boundary springs (clamped, simply supported, free as in Table 1)
[K,M]=springBC(A,B,D_mat,As,R1,K, M, nodes, ndof,bc_type,Phi,supp);
%% 5. Solve eigenproblem
eigMode = 20;
% opts.issym = 1; opts.isreal=1;
[V,D] = eigs(K,M,eigMode,'SM');
freq = sqrt(diag(D))/(2*pi);
[~, idx] = sort(freq);
freq = freq(idx);
V = V(:, idx);
% % 过滤非物理解
% valid_idx = imag(freq) == 0 & freq > 0 & freq < 1e10;
% freq = freq(valid_idx);
% V = V(:, valid_idx);
disp('Natural frequencies (Hz):'); disp(freq);
%% --- Functions ---
function Rb = radialRB(r, type, dc, q)
switch type
case 'MQ', Rb = (r^2 + dc^2)^q;
end
end
function Tm = tchebyBasis(xp,tp,mp,mq)
m=mp*mq;
if m==0
Tm=[];
else
Tm = zeros(length(xp),m);
for i=1:length(xp)
x=xp(i); t=tp(i);
Te = zeros(1,m);
% 计算二维切比雪夫基
for P=0:mp-1
for Q=0:mq-1
PQ=P*mq+Q+1;
Te(PQ)=cos(Q*acos(t))*cos(P*acos(x));
end
end
Tm(i,:)=Te;
end
end
end
function Phi = computeShape(suppCoord, evalCoord, Ginv, rbf_type, dc, q, mp,mq, L, phiSpan, R)
nn = size(suppCoord,1);
% distances & derivatives at eval
Re = zeros(1,nn); dRe_x=zeros(1,nn); dRe_t=zeros(1,nn);
dRe_xx=zeros(1,nn); dRe_xt=zeros(1,nn); dRe_tt=zeros(1,nn);
for i=1:nn
dx = evalCoord(1)-suppCoord(i,1);
% dt = evalCoord(2)-suppCoord(i,2);
% r = sqrt(dx^2 - R^2*(2*cos(dt)-2));
dt = wrapToPi((evalCoord(2)-suppCoord(i,2)));
r = sqrt(dx^2 + dt^2*R^2);
switch rbf_type
case 'MQ'
% Re(i) = (r^2 + dc^2)^q;
% dRe_x(i) = 2*q*(r^2 + dc^2)^(q-1)*dx;
% dRe_t(i) = 2*q*(r^2 + dc^2)^(q-1)*R^2*sin(dt);
% dRe_xx(i) = 2*q*(r^2 + dc^2)^(q-1) + 4*(q-1)*q*dx^2*(r^2 + dc^2)^(q-2);
% dRe_xt(i) = 4*(q-1)*q*dx*R^2*sin(dt)*(r^2 + dc^2)^(q-2);
% dRe_tt(i) = 2*q*(r^2 + dc^2)^(q-1)*R^2*cos(dt)+ 4*q*(q-1)*R^4*(sin(dt))^2*(r^2 + dc^2)^(q-2);
Re(i) = (r^2 + dc^2)^q;
dRe_x(i) = 2*q*(r^2 + dc^2)^(q-1)*dx;
dRe_t(i) = 2*q*(r^2 + dc^2)^(q-1)*R^2*dt;
dRe_xx(i) = 2*q*(r^2 + dc^2)^(q-1) + 4*(q-1)*q*dx^2*(r^2 + dc^2)^(q-2);
dRe_xt(i) = 4*(q-1)*q*dx*R^2*dt*(r^2 + dc^2)^(q-2);
dRe_tt(i) = 2*q*(r^2 + dc^2)^(q-1)*R^2+ 4*q*(q-1)*R^4*dt^2*(r^2 + dc^2)^(q-2);
end
end
% Tcheby at eval
x = 2*(evalCoord(1)-0)/L -1;
t = 2*(evalCoord(2)-0)/phiSpan -1;
m=mp*mq;
Te = zeros(1,m); dTe_x=zeros(1,m); dTe_t=zeros(1,m);
dTe_xx=zeros(1,m); dTe_xt=zeros(1,m); dTe_tt=zeros(1,m);
% 计算切比雪夫多项式
[Tx, dTx_norm, d2Tx_norm] = chebyshevT_derivatives(0:mp-1, x);
[Ty, dTy_norm, d2Ty_norm] = chebyshevT_derivatives(0:mq-1, t);
% 应用坐标变换的导数修正
dTx = dTx_norm * 2/L;
d2Tx = d2Tx_norm * (2/L)^2;
dTy = dTy_norm * 2/phiSpan;
d2Ty = d2Ty_norm * (2/phiSpan)^2;
% 计算二维切比雪夫基
idx=1;
for P=1:mp
for Q=1:mq
Te(idx)= Tx(P)*Ty(Q);
dTe_x(idx)= dTx(P)*Ty(Q);
dTe_t(idx)= Tx(P)*dTy(Q);
dTe_xx(idx)= d2Tx(P)*Ty(Q);
dTe_xt(idx)= dTx(P)*dTy(Q);
dTe_tt(idx)= Tx(P)*d2Ty(Q);
idx = idx + 1;
end
end
% assemble Phi at eval
% Phi_bar(1,:) = round([Re,Te]*Ginv,5);
Phi_bar(1,:) = [Re,Te]*Ginv;
Phi_bar(2,:) = [dRe_x,dTe_x]*Ginv;
Phi_bar(3,:) = [dRe_t,dTe_t]*Ginv;
Phi_bar(4,:) = [dRe_xx,dTe_xx]*Ginv;
Phi_bar(5,:) = [dRe_xt,dTe_xt]*Ginv;
Phi_bar(6,:) = [dRe_tt,dTe_tt]*Ginv;
Phi=Phi_bar(:,1:nn);
end
function [T, dT, d2T] = chebyshevT_derivatives(n, x)
% 计算切比雪夫多项式 T_n(x) 及其一阶、二阶导数
x = x(:); % 转换为列向量
T = zeros(length(x), length(n));
dT = zeros(length(x), length(n));
d2T = zeros(length(x), length(n));
for i = 1:length(n)
order = n(i);
if order == 0
T(:,i) = ones(size(x));
dT(:,i) = zeros(size(x));
d2T(:,i) = zeros(size(x));
elseif order == 1
T(:,i) = x;
dT(:,i) = ones(size(x));
d2T(:,i) = zeros(size(x));
else
% 初始化递推
T0 = ones(size(x));
T1 = x;
dT0 = zeros(size(x));
dT1 = ones(size(x));
d2T0 = zeros(size(x));
d2T1 = zeros(size(x));
% 递推计算
for k = 2:order
T2 = 2*x.*T1 - T0;
dT2 = 2*T1 + 2*x.*dT1 - dT0;
d2T2 = 4*dT1 + 2*x.*d2T1 - d2T0;
% 更新
T0 = T1;
T1 = T2;
dT0 = dT1;
dT1 = dT2;
d2T0 = d2T1;
d2T1 = d2T2;
end
T(:,i) = T1;
dT(:,i) = dT1;
d2T(:,i) = d2T1;
end
end
end
function [Ki,Mi] = localKM(Phi, A, B, D, As, R, I0, I1, I2)
ndof=5; nn=size(Phi,2);
Ki = zeros(ndof,ndof,nn); Mi = zeros(ndof,ndof,nn);
for j=1:nn
phi=Phi(1,j);
dphidx=Phi(2,j);
dphidt=Phi(3,j);
dphidxx=Phi(4,j);
dphidxt=Phi(5,j);
dphidtt=Phi(6,j);
% K
Tn11 = A(1,1) * dphidxx + (A(3,3)/R^2) * dphidtt + (2*A(1,3)/R) * dphidxt;
Tn12 = A(1,3) * dphidxx + (A(1,2)/R) * dphidxt + (A(3,3)/R) * dphidxt + (A(2,3)/R^2) * dphidtt;
Tn13 = (A(1,2)/R) * dphidx + (A(2,3)/R^2) * dphidt;
Tn14 = B(1,1) * dphidxx + (2*B(1,3)/R) * dphidxt + (B(3,3)/R^2) * dphidtt;
Tn15 = B(1,3) * dphidxx + (B(1,2)/R) * dphidxt + (B(3,3)/R) * dphidxt + (B(2,3)/R^2) * dphidtt;
Tn21 = Tn12;
Tn22 = A(3,3) * dphidxx + (A(2,2)/R^2) * dphidtt + (2*A(2,3)/R) * dphidxt - (As(1,1)/R^2) * phi;
Tn23 = (A(2,3)/R) * dphidx + (A(2,2)/R^2) * dphidt + (As(1,2)/R) * dphidx + (As(1,1)/R^2) * dphidt;
Tn24 = B(1,3) * dphidxx + (B(1,2)/R) * dphidxt + (B(3,3)/R) * dphidxt + (B(2,3)/R^2) * dphidtt + (As(1,2)/R) * phi;
Tn25 = B(3,3) * dphidxx + (2*B(2,3)/R) * dphidxt + (B(2,2)/R^2) * dphidtt + (As(1,1)/R) * phi;
Tn31 = -Tn13;
Tn32 = -Tn23;
Tn33 = As(2,2) * dphidxx - (A(2,2)/R^2) * phi + (As(1,1)/R^2) * dphidtt + (2*As(1,2)/R) * dphidxt;
Tn34 = As(2,2) * dphidx - (B(1,2)/R) * dphidx - (B(2,3)/R^2) * dphidt + (As(1,2)/R) * dphidt;
Tn35 = As(1,2) * dphidx - (B(2,3)/R) * dphidx - (B(2,2)/R^2) * dphidt + (As(1,1)/R) * dphidt;
Tn41 = Tn14;
Tn42 = Tn24;
Tn43 = -Tn34;
Tn44 = D(1,1) * dphidxx - As(2,2) * phi + (D(3,3)/R^2) * dphidtt + (2*D(1,3)/R) * dphidxt;
Tn45 = D(1,3) * dphidxx + (D(1,2)/R) * dphidxt + (D(3,3)/R) * dphidxt + (D(2,3)/R^2) * dphidtt - As(1,2) * phi;
Tn51 = Tn15;
Tn52 = Tn25;
Tn53 = -Tn35;
Tn54 = Tn45;
Tn55 = D(3,3) * dphidxx - As(1,1) * phi + (D(2,2)/R^2) * dphidtt + (2*D(2,3)/R) * dphidxt;
Kij= [Tn11, Tn12, Tn13, Tn14, Tn15;
Tn21, Tn22, Tn23, Tn24, Tn25;
Tn31, Tn32, Tn33, Tn34, Tn35;
Tn41, Tn42, Tn43, Tn44, Tn45;
Tn51, Tn52, Tn53, Tn54, Tn55];
Ki(:,:,j)=Kij;
% M
Mij= -[I0*phi,0,0,I1*phi,0;
0,I0*phi,0,0,I1*phi;
0,0,I0*phi,0,0;
I1*phi,0,0,I2*phi,0;
0,I1*phi,0,0,I2*phi];
Mi(:,:,j)=Mij;
end
end
function [K,M]=springBC(A,B,D,As,R,K,M,nodes,ndof,bc_type,allPhi,supp)
% Example: clamp x=0 and x=L, free elsewhere
switch bc_type
case 'F' %自由
ku=0;kv=0;kw=0;kx=0;kt=0;
case 'S' %简支
ku=10^14;kv=10^14;kw=10^14;kx=0;kt=10^14;
case 'SD' %弹性简支
ku=0;kv=10^14;kw=10^14;kx=0;kt=0;
case 'C' %固支
ku=10^14;kv=10^14;kw=10^14;kx=10^14;kt=10^14;
case 'E' %任意弹性
ku=10^7;kv=10^14;kw=10^14;kx=10^14;kt=10^14;
end
tol = 1e-8;
N = size(nodes,1);
for i=1:N
Phi=allPhi{i};
phi=Phi(1,:);
dphidx=Phi(2,:);
dphidt=Phi(3,:);
x = nodes(i,1);
t = nodes(i,2);
indices=supp{i};
nn=length(indices);
if abs(x)<tol
for j=1:nn
Cxn=calculate_Cxn(A,B,D,As,R,phi(j),dphidx(j),dphidt(j));
k_spring=blkdiag(ku,kv,kw,kx,kt)* phi(j);
dofl = (indices(j)-1)*ndof+1:indices(j)*ndof;
dofh = (i-1)*ndof+1:i*ndof;
K(dofh,dofl)=0;
M(dofh,dofl)=0;
K(dofh,dofl)=Cxn-k_spring;
end
elseif abs(x- max(nodes(:,1)))<tol
for j=1:nn
Cxn=calculate_Cxn(A,B,D,As,R,phi(j),dphidx(j),dphidt(j));
k_spring=blkdiag(ku,kv,kw,kx,kt)* phi(j);
dofl = (indices(j)-1)*ndof+1:indices(j)*ndof;
dofh = (i-1)*ndof+1:i*ndof;
K(dofh,dofl)=0;
M(dofh,dofl)=0;
K(dofh,dofl)=Cxn+k_spring;
end
elseif abs(t)<tol
for j=1:nn
Ctn=calculate_Ctn(A,B,D,As,R,phi(j),dphidx(j),dphidt(j));
k_spring=blkdiag(ku,kv,kw,kx,kt)* phi(j);
dofl = (indices(j)-1)*ndof+1:indices(j)*ndof;
dofh = (i-1)*ndof+1:i*ndof;
K(dofh,dofl)=0;
M(dofh,dofl)=0;
K(dofh,dofl)=Ctn-k_spring;
end
elseif abs(t- max(nodes(:,2)))<tol
for j=1:nn
Ctn=calculate_Ctn(A,B,D,As,R,phi(j),dphidx(j),dphidt(j));
k_spring=blkdiag(ku,kv,kw,kx,kt)* phi(j);
dofl = (indices(j)-1)*ndof+1:indices(j)*ndof;
dofh = (i-1)*ndof+1:i*ndof;
K(dofh,dofl)=0;
M(dofh,dofl)=0;
K(dofh,dofl)=Ctn+k_spring;
end
end
end
end
function Cxn=calculate_Cxn(A,B,D,As,R,phi,dphidx,dphidt)
Cxn11 = A(1,1) * dphidx + A(1,3)/R * dphidt;
Cxn12 = A(1,3) * dphidx + A(1,2)/R * dphidt;
Cxn13 = A(1,2)/R * phi;
Cxn14 = B(1,1) * dphidx + B(1,3)/R * dphidt;
Cxn15 = B(1,3) * dphidx + B(1,2)/R * dphidt;
Cxn21 = A(1,3) * dphidx + A(3,3)/R * dphidt;
Cxn22 = A(3,3) * dphidx + A(2,3)/R * dphidt;
Cxn23 = A(2,3)/R * phi;
Cxn24 = B(1,3) * dphidx + B(3,3)/R * dphidt;
Cxn25 = B(3,3) * dphidx + B(2,3)/R * dphidt;
Cxn31 = 0;
Cxn32 = -As(1,2)/R * phi;
Cxn33 = As(2,2) * dphidx + As(1,2)/R * dphidt;
Cxn34 = As(2,2) * phi;
Cxn35 = As(1,2) * phi;
Cxn41 = Cxn14;
Cxn42 = Cxn15;
Cxn43 = B(1,2)/R * phi;
Cxn44 = D(1,1) * dphidx + D(1,3)/R * dphidt;
Cxn45 = D(1,3) * dphidx + D(1,2)/R * dphidt;
Cxn51 = Cxn24;
Cxn52 = Cxn25;
Cxn53 = B(2,3)/R * phi;
Cxn54 = D(1,3) * dphidx + D(3,3)/R * dphidt;
Cxn55 = D(3,3) * dphidx + D(2,3)/R * dphidt;
Cxn= [Cxn11, Cxn12, Cxn13, Cxn14, Cxn15;
Cxn21, Cxn22, Cxn23, Cxn24, Cxn25;
Cxn31, Cxn32, Cxn33, Cxn34, Cxn35;
Cxn41, Cxn42, Cxn43, Cxn44, Cxn45;
Cxn51, Cxn52, Cxn53, Cxn54, Cxn55];
end
function Ctn=calculate_Ctn(A,B,D,As,R,phi,dphidx,dphidt)
Ctn11 = A(1,3) * dphidx + A(3,3)/R * dphidt;
Ctn12 = A(3,3) * dphidx + A(2,3)/R * dphidt;
Ctn13 = A(2,3)/R * phi;
Ctn14 = B(1,3) * dphidx + B(3,3)/R * dphidt;
Ctn15 = B(3,3) * dphidx + B(2,3)/R * dphidt;
Ctn21 = A(1,2) * dphidx + A(2,3)/R * dphidt;
Ctn22 = A(2,3) * dphidx + A(2,2)/R * dphidt;
Ctn23 = A(2,2)/R * phi;
Ctn24 = B(1,2) * dphidx + B(2,3)/R * dphidt;
Ctn25 = B(2,3) * dphidx + B(2,2)/R * dphidt;
Ctn31 = 0;
Ctn32 = -As(1,1)/R * phi;
Ctn33 = As(1,2) * dphidx + As(1,1)/R * dphidt;
Ctn34 = As(1,2) * phi;
Ctn35 = As(1,1) * phi;
Ctn41 = Ctn14;
Ctn42 = Ctn15;
Ctn43 = B(2,3)/R * phi;
Ctn44 = D(1,3) * dphidx + D(3,3)/R * dphidt;
Ctn45 = D(3,3) * dphidx + D(2,3)/R * dphidt;
Ctn51 = Ctn24;
Ctn52 = Ctn25;
Ctn53 = B(2,2)/R * phi;
Ctn54 = D(1,2) * dphidx + D(2,3)/R * dphidt;
Ctn55 = D(2,3) * dphidx + D(2,2)/R * dphidt;
Ctn= [Ctn11, Ctn12, Ctn13, Ctn14, Ctn15;
Ctn21, Ctn22, Ctn23, Ctn24, Ctn25;
Ctn31, Ctn32, Ctn33, Ctn34, Ctn35;
Ctn41, Ctn42, Ctn43, Ctn44, Ctn45;
Ctn51, Ctn52, Ctn53, Ctn54, Ctn55];
end
function [A,B,D,As,I] = laminateABD(layup,E1,E2,mu12,G12,G13,G23,h,rho,ks)
Nk = numel(layup); z = linspace(-h/2,h/2,Nk+1);
A=zeros(3); B=zeros(3); D=zeros(3); As = zeros(2); I = [0, 0, 0];
for k=1:Nk
th=layup(k)*pi/180; c=cos(th); s=sin(th);
[Q,As0] = orthotropicQ(E1,E2,mu12,G12,G13,G23,ks);
T = [c^2, s^2, -2*c*s; s^2, c^2, 2*c*s; c*s, -c*s, c^2-s^2];
Qbar = T*Q*T'; % transform
Ts = [c, s; -s, c];
Asbar = Ts * As0 * Ts';
A = A + Qbar*(z(k+1)-z(k));
B = B + Qbar*(z(k+1)^2-z(k)^2)/2;
D = D + Qbar*(z(k+1)^3-z(k)^3)/3;
As = As + Asbar * (z(k+1) - z(k));
I = I + rho * [(z(k+1) - z(k)), (z(k+1)^2 - z(k)^2)/2, (z(k+1)^3 - z(k)^3)/3];
% A=round(A,5);
% B=round(B,5);
% D=round(D,5);
% As=round(As,5);
end
end
function [Q,As] = orthotropicQ(E1,E2,mu12,G12,G13,G23,ks)
mu21 = E2*mu12/E1;
Q = [E1/(1-mu12*mu21), mu12*E2/(1-mu12*mu21), 0;
mu12*E2/(1-mu12*mu21), E2/(1-mu12*mu21), 0;
0,0,G12];
As = [ks*G23, 0;
0, ks*G13];
end
function f=draw_shell(nodes)
f=figure;
% 生成示例数据
x = nodes(:,1);
y = nodes(:,2);
z = nodes(:,3);
% 绘制三维点
scatter3(x, y, z, 'filled');
hold on;
% 添加标题和坐标轴标签
title('三维点图');
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
% 设置坐标轴比例一致
axis equal;
end
这是我的代码，你帮我分析分析
本文探讨了如何通过矩阵乘法来递推求解特定序列，并利用群论知识找到循环节，从而解决大指数问题。文中还讨论了几种不同的取模策略及其优劣。
扫一扫
5928
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
