这份代码不能正常折叠折纸的原因：clc
clear
close all
%% 结构参数
a = 1; % 四边形边长
ad = 1; % 虚拟点高度系数
th = 60; % 边夹角
h = ad * a; % 虚拟点高度
dt = 0.01; % 时间增量
steps = 500; % 步数
%% 节点定义 (13个节点：9实体+4虚拟)
nodes = zeros(13, 3);
% 实体节点 (z=0)
nodes(1,:) = [0, 0, 0]; % 节点0
nodes(2,:) = [a*cosd(th), a*sind(th), 0]; % 节点1
nodes(3,:) = [2*a*cosd(th), 2*a*sind(th), 0]; % 节点2
nodes(4,:) = [a, 0, 0]; % 节点3
nodes(5,:) = [a + a*cosd(th), a*sind(th), 0]; % 节点4
nodes(6,:) = [a + 2*a*cosd(th), 2*a*sind(th), 0]; % 节点5
nodes(7,:) = [a + a*cosd(th), -a*sind(th), 0]; % 节点6
nodes(8,:) = [2*a, 0, 0]; % 节点7
nodes(9,:) = [2*a + a*cosd(th), a*sind(th), 0]; % 节点8
% 虚拟节点 (z=h)
nodes(10,:) = [(a + a*cosd(th))/2, a*sind(th)/2, h]; % 节点9
nodes(11,:) = [(a + a*cosd(th))/2 + a*cosd(th), a*sind(th)/2 + a*sind(th), h]; % 节点10
nodes(12,:) = [3*a/2, 0, h]; % 节点11
nodes(13,:) = [3*a/2 + a*cosd(th), a*sind(th), h]; % 节点12
%% 定义杆件 (32条边)
edges = [
% 实体边 (16)
0 1; 1 2; 3 4; 4 5; 6 7; 7 8;
0 3; 3 6; 1 4; 4 7; 2 5; 5 8;
1 3; 3 7; 2 4; 4 8;
% 虚拟连接 (16)
0 9; 1 9; 3 9; 4 9;
1 10; 2 10; 4 10; 5 10;
3 11; 4 11; 6 11; 7 11;
4 12; 5 12; 7 12; 8 12] + 1; % 节点索引从1开始，所以加1
%% 定义折痕线（旋转弹簧） - 每条折痕线由4个节点定义 [a,c,b,d]
folds = [2,5,1,3;
4,5,1,7;
5,6,3,9;
5,8,4,6;
];
% 设置折叠刚度 (山折/谷折不同)
k_f = 0; % 谷折刚度 (小)
k_b = 0; % 山折刚度 (大)
fold_stiffness = [k_f, k_b, k_b, k_f]; % 每条折痕的刚度
%% 物理参数
EA = 1e8; % 杆件轴向刚度 (大值模拟刚性)
k_rot = fold_stiffness; % 旋转弹簧刚度数组
theta0 = pi; % 初始二面角 (平铺状态)
tol = 1e-6; % 牛顿迭代容差
max_iter = 100; % 最大迭代次数
%% 边界条件 (固定部分节点)
% 固定节点方向
fix_dofs = [3,21,27,9,13,14];
% 自由自由度
all_dofs = 1:3*size(nodes,1);
free_dofs = setdiff(all_dofs, fix_dofs);
%% 施加外力 (在节点4的z方向施加力)
F_ext = zeros(3*size(nodes,1), 1);
% 在节点5（索引5）的z方向施加负向力（向下）
node_force = 5; % 节点索引
F_ext(3*node_force) = 10;
%% 主循环 (伪时间步，用于逐步施加载荷)
noden = nodes; % 当前节点坐标
U_his = zeros(3*size(nodes,1), steps); % 记录位移历史
for step = 1:steps
% 逐步增加载荷
lambda = step / steps;
F_current = lambda * F_ext;
% 牛顿迭代
residual_norm = inf;
iter = 0;
noden_prev = noden; % 保存上一步的节点坐标
while residual_norm > tol && iter < max_iter
% 组装全局刚度矩阵和残差向量
[K_global, R] = assemble_system(noden, edges, folds, EA, k_rot, theta0);
% 计算残差（包括外部力）
R_total = R - F_current;
% 处理边界条件
K_red = K_global(free_dofs, free_dofs);
R_red = R_total(free_dofs);
% 求解线性系统
du_red = K_red \ R_red;
% 更新自由自由度位移
du_full = zeros(3*size(nodes,1), 1);
du_full(free_dofs) = du_red;
% 更新节点坐标
for i = 1:size(noden, 1)
idx = 3*(i-1)+1:3*i;
noden(i,:) = noden(i,:) + du_full(idx)';
end
% 检查收敛
residual_norm = norm(R_red);
iter = iter + 1;
end
% 记录位移
U_his(:, step) = noden(:) - nodes(:);
% 可视化
if mod(step, 10) == 0
visualize_structure(noden, edges, folds, step);
end
end
%% 绘制折叠角度变化
figure(2);
hold on;
for i = 1:size(folds,1)
angles = [];
for step = 1:steps
n = folds(i,:);
a = noden(n(1),:);
c = noden(n(2),:);
b = noden(n(3),:);
d = noden(n(4),:);
theta = calculate_dihedral_angle(a, c, b, d);
angles(step) = theta;
end
plot(1:steps, rad2deg(angles), 'LineWidth', 1.5);
end
title('折痕角度变化');
xlabel('时间步');
ylabel('角度 (度)');
legend('折痕线1', '折痕线2', '折痕线3', '折痕线4', '折痕线5');
grid on;
%% 系统组装函数
function [K_global, R] = assemble_system(nodes, edges, folds, EA, k_rot, theta0)
n_nodes = size(nodes, 1);
n_dofs = 3 * n_nodes;
% 初始化全局矩阵
K_global = zeros(n_dofs);
R = zeros(n_dofs, 1);
% 处理杆件单元
for i = 1:size(edges, 1)
n1 = edges(i, 1);
n2 = edges(i, 2);
[ke, re] = bar_element(nodes(n1,:), nodes(n2,:), EA);
% 组装到全局系统
dofs = [3*(n1-1)+1:3*n1, 3*(n2-1)+1:3*n2];
K_global(dofs, dofs) = K_global(dofs, dofs) + ke;
R(dofs) = R(dofs) + re;
end
% 处理旋转弹簧单元
for i = 1:size(folds, 1)
n = folds(i, :); % [a,c,b,d]
[ke_rot, re_rot] = rotational_spring_element(...
nodes(n(1),:), nodes(n(2),:), nodes(n(3),:), nodes(n(4),:), ...
k_rot(i), theta0);
% 组装到全局系统
dofs = [];
for j = 1:4
dofs = [dofs, 3*(n(j)-1)+1:3*n(j)];
end
K_global(dofs, dofs) = K_global(dofs, dofs) + ke_rot;
R(dofs) = R(dofs) + re_rot;
end
end
%% 杆件单元计算 (线性)
function [ke, re] = bar_element(n1, n2, EA)
vec = n2 - n1;
L0 = norm(vec);
dir = vec / L0;
% 当前长度
L = norm(vec);
% 应变
strain = (L - L0) / L0;
% 力
force = EA * strain;
% 局部刚度矩阵 (线性)
ke_local = EA/L0 * [1, -1; -1, 1];
% 转换矩阵 (从局部坐标到全局坐标)
T = [dir; null(dir)']'; % 局部坐标系的基
T_blk = blkdiag(T, T);
ke = T_blk' * kron(ke_local, eye(3)) * T_blk;
% 残差向量 (杆端力)
re = force * [-dir, dir]';
end
%% 旋转弹簧单元计算 - 完整实现
function [ke_rot, re_rot] = rotational_spring_element(a, c, b, d, k_rot, theta0)
% 计算二面角
theta = calculate_dihedral_angle(a, c, b, d);
% 计算力矩（线性模型）
M = k_rot * (theta - theta0);
% 计算梯度 (∂θ/∂u)
grad_theta = zeros(1, 12); % 12个自由度（4个节点*3）
epsilon = 1e-6; % 数值微分步长
% 数值计算梯度
for i = 1:4
for j = 1:3
% 复制节点坐标
nodes_copy = {a, c, b, d};
% 扰动第i个节点的第j个坐标
nodes_copy{i}(j) = nodes_copy{i}(j) + epsilon;
% 重新计算二面角
theta_plus = calculate_dihedral_angle(nodes_copy{1}, nodes_copy{2}, ...
nodes_copy{3}, nodes_copy{4});
% 计算偏导数
grad_index = (i-1)*3 + j;
grad_theta(grad_index) = (theta_plus - theta) / epsilon;
end
end
% 残差向量 (力矩对应的广义力)
re_rot = M * grad_theta';
% 切线刚度矩阵
ke_rot = k_rot * (grad_theta' * grad_theta);
end
%% 计算二面角函数
function theta = calculate_dihedral_angle(a, c, b, d)
% 向量ac
vec_ac = c - a;
% 向量ab
vec_ab = b - a;
% 向量ad
vec_ad = d - a;
% 两个面的法向量
n1 = cross(vec_ac, vec_ab);
n2 = cross(vec_ac, vec_ad);
n1 = n1 / norm(n1);
n2 = n2 / norm(n2);
% 计算二面角（两个法向量的夹角）
cos_theta = dot(n1, n2);
% 确保角度在有效范围内
cos_theta = max(min(cos_theta, 1), -1);
theta = acos(cos_theta);
% 调整角度范围到0到π之间
if theta > pi
theta = 2*pi - theta;
end
end
%% 可视化函数
function visualize_structure(nodes, edges, folds, step)
figure(1);
clf;
% 绘制节点
plot3(nodes(:,1), nodes(:,2), nodes(:,3), 'bo', 'MarkerSize', 6, 'MarkerFaceColor', 'b');
hold on;
% 绘制杆件
for i = 1:size(edges, 1)
n1 = edges(i, 1);
n2 = edges(i, 2);
plot3([nodes(n1,1), nodes(n2,1)], [nodes(n1,2), nodes(n2,2)], [nodes(n1,3), nodes(n2,3)], 'k-', 'LineWidth', 1.5);
end
% 绘制折痕线 (红色高亮)
for i = 1:size(folds, 1)
n = folds(i, [1,2]); % a-c轴
plot3([nodes(n(1),1), nodes(n(2),1)], [nodes(n(1),2), nodes(n(2),2)], [nodes(n(1),3), nodes(n(2),3)], 'r-', 'LineWidth', 2.5);
end
% 绘制面板
panels = {
[1, 2, 5, 4], % 左下
[2, 3, 6, 5], % 左上
[5, 6, 9, 8], % 右上
[4, 5, 8, 7] % 右下
};
colors = [0.8, 0.8, 1.0; % 浅蓝
1.0, 0.8, 0.8; % 浅红
0.8, 1.0, 0.8; % 浅绿
1.0, 1.0, 0.8]; % 浅黄
for i = 1:length(panels)
verts = nodes(panels{i}, :);
patch('Vertices', verts, 'Faces', [1,2,3,4], ...
'FaceColor', colors(i,:), 'FaceAlpha', 0.6, ...
'EdgeColor', 'k', 'LineWidth', 1.2);
end
% 标注
title(sprintf('非均匀折叠模拟\n时间步: %d', step), 'FontSize', 14);
xlabel('X'); ylabel('Y'); zlabel('Z');
axis equal;
grid on;
view(3);
rotate3d on;
set(gca, 'FontSize', 12);
drawnow;
end
此篇博客介绍了如何通过`get_normals`函数获取边缘向量的法向量，利用向量点积和弧度转换计算两个面之间的二面角。重点讲解了axis=0和clip(-1,1)的作用，以及np.expand_dims的使用。适合理解三维空间中几何计算的读者。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
