FEM-Course-Matlab/16.几何非线性有限元matlab编程/几何非线性有限元-Williams frame/intForcesUL.m

55 lines
1.7 KiB
Matlab
Raw Permalink Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

function [F,Elem] = intForcesUL(Model,Elem,U,D_U,update_angle)
% Initialize global vector of internal forces
F = zeros(Model.neq,1);
for i = 1:Model.nel
% Lengths: Beginning of step, current, and step increment
L_1 = Elem(i).L_1;%每个增量步开始时进行更新
L_c = elemLength(Elem(i),U);
D_L = L_c - L_1;
% 刚体位移角度增量(相对于增量步初始状态)
rbr = elemAngleIncr(Elem(i),U,D_U);%rbr带有符号增加正值减少负值
%rbr为相对荷载步初始转角的单元刚性转角即当前构型与荷载步初始时的构型的整体转角差梁单元轴向角度
angle = Elem(i).angle_1 + rbr;
% 更新单元刚体转角
if (update_angle)
Elem(i).angle = angle;
end
% 局部坐标系-全局坐标系的坐标转换矩阵
c = cos(angle);
s = sin(angle);
rot = [ c -s 0 0 0 0;
s c 0 0 0 0;
0 0 1 0 0 0;
0 0 0 c -s 0;
0 0 0 s c 0;
0 0 0 0 0 1 ];
% 形变转角,相对增量步初始构型
r1 = D_U(Elem(i).n1.dof(3)) - rbr;
r2 = D_U(Elem(i).n2.dof(3)) - rbr;
% 局部坐标下的位移向量增量,相对增量步初始构型
dl = [0; 0 ; r1; D_L; 0; r2];%没有y向变形只有轴向和转动变形因为y向的变形其实等同于转角局部坐标系下不存在y向变形
% 局部坐标系下内力向量增量
D_fl = Elem(i).ke * dl;
% 局部坐标系下总内力
fl = Elem(i).fi_1 + D_fl;
% 更新Elem中局部坐标系总内力
Elem(i).fi = fl;
% 内力从局部到全局的转换
fg = rot * fl;
% 单元内力组装到整体内力矩阵中
gle = Elem(i).gle;
F(gle) = F(gle) + fg;
end
end