一阶四面体刚度矩阵计算

目录概述杆单元为什么要引入局部坐标局部坐标系下的刚度矩阵全局坐标系下的刚度矩阵程序编写&算例所有的计算代码概述   杆系结构是工程中应用较为广泛的结构体系，包含了平面或空间形式的梁、钢架和行桁架等等，这些结构的每个单独的杆件都可以视为一个单独的单元，单元与单元之间用节点连接，根据结构体系的不同，节点可以传递轴向力、切向力或者弯矩等等。   帖子以二维平面杆单元为例，采用材料力学更直观的方法建立杆单元的刚度矩阵，详细介绍了具体的推导过程，最后以一个桁架为例讲解matlab程序并采用abaqus做对比计算，文末给出了详细的程序。杆单元为什么要引入局部坐标系  我们在做一件事情的时候，一定要知道为什么做。  杆单元为什么不直接在全局坐标系中求解，非得人为的引入一个局部坐标系？我想第一次学习杆单元的人都会有这个疑问，然后很遗憾，单单用杆单元部分的知识并不能很好的回答这个问题。一个简单的答案是：为了统一单元、最终是为了简化计算。   我们知道，科学家最喜欢将特殊的东西泛化，然后给出一个简洁美丽的公式来解释广泛的事物。在科学计算中也是如此，就拿这次的杆单元来说，我们脑子里面浮现的是一个结构体系中横七竖八的段杆，就像是体育场上上空的架子，乱糟糟的，用数学的语言描述就是他们的空间位置都不相同，这在直觉上意味着极大的工作量。   这时候，天马行空的想法随着工作量的增加会自然而然的从脑子里面冒出来（很多发明一开始都是天马行空的偷懒想法）：如果我们能只求一次某个量，然后乘以一个系数把这个量映射到所有的杆单元上面就好了，这个想法就是引入局部坐标系的火苗，有了这个思想的指导，后面的数学家就开始干活了，最终的杆单元计算过程就是上面那句加粗的文字那样：首先在局部坐标系下计算刚度矩阵，然后乘以一个坐标转换矩阵，将据局部坐标系下的刚度矩阵映射到全局坐标系中参与组装。下面列出一些引入局部坐标系具体方便了哪些简化计算    局部坐标系可以简化杆单元的几何和材料属性的描述，使得在进行刚度矩阵和质量矩阵的计算时更为方便。通过将杆单元的方向与局部坐标系对齐，可以使得计算过程中的数学表达更为直观。方便载荷施加    在局部坐标系中，可以更方便地定义和施加载荷。比如，集中载荷或分布载荷可以直接在杆单元的局部坐标系中施加，避免了转换成全球坐标系的复杂性。更好地处理变形和应力分析    杆单元在变形时会有不同方向的应力和应变，局部坐标系可以帮助明确杆单元的拉伸、压缩及剪切行为，进而提高计算精度。支持多方向的连接    在复杂的结构中，杆单元可能以不同的方向连接。在局部坐标系中，每个单元的方向都是独立定义的，这使得处理复杂结构的连接关系更为灵活。易于与其他单元类型结合    在实际工程应用中，结构通常由多种单元组成。引入局部坐标系可以更方便地将杆单元与其他类型的单元（如面单元或体单元）结合，确保它们之间的力学性质和边界条件协调一致。   我在这里只是从定性的角度给出了一些解释，要想彻底的从底层理解这一套工作思路，还需要学习等参元，等参元的引入把简化计算的思想发挥的淋漓尽致，在等参元中，科学家通过坐标变换解决了积分难题，局部坐标系下的刚度矩阵   材料力学给出了二维杆单元的受力特点每个杆单元的节点只有两个平动自由度二维杆单元的节点不传递弯矩外荷载均为作用在节点上的集中力    基于此，我们从一个结构体系中取出一个杆件，并引入与杆轴向重合的局部坐标系，设杆单元的截面积为 ，长度为 ，两端的节点分别用 和 表示，并且从节点 到节点 的方向为局部有坐标轴 ，垂直于杆轴向为 ，具体的示意图为   根据材料力学的知识，对两端节点列方程   整理成矩阵形式   以上就是二维杆单元在局部坐标系中的刚度矩阵，下面推导局部坐标系与全局坐标系的转化矩阵，即全局坐标系的刚度矩阵。全局坐标系下的刚度矩阵  推导全局坐标系下的刚度矩阵的重点就是推导转换矩阵，即轴向量与全局坐标系之间的转换关系，如下图展示了轴向数据与全局坐标系数据   通过上面图片可知，整体坐标系下，杆单元两端节点的位移矩阵为   局部坐标系下，杆单元两端节点的位移矩阵为    根据三角函数关系，写出全局坐标系位移与局部坐标系位移之间的关系为    整理成矩阵形式   其中，转换矩阵 为   其中   上面就是局部坐标系向全局坐标系转换数据的矩阵，局部和全局量尊遵循上述转换关系。对于全局坐标系下的刚度矩阵，则有   结合上文中的局部坐标系下的刚度矩阵，代入转换矩阵有   好了，现在有了全局坐标系下的杆单元刚度矩阵，下面就可以写程序了。程序编写&算例  这里直接给出一个算例，以算例的形式详细讲解程序。二维杆单元的弹性模量 ，截面面积 ，具体的网格和边界条件示意图为  上面的桁架一共15个杆单元，单元编号用斜体字表示；9个节点，节点用有背景色的方块表示，在节点1和6上施加竖直向的集中力，力的大小 ，节点9和4为的两个平动自由度被限制。   下面开始介绍关键的代码，全部的代码在文末给出。   首先是杆单元的材料属性，杆单元只需要弹性模量和截面积。E=210e3; % 弹性模量 (GPa=1000 N/mm^2)A=500; % 杆单元截面积 (mm^2)    然后给出模型的信息。node=9; % 模型节点总数element=15; % 杆单元的数目nodecord=[ 1.5, 0.866025388; 2.5, 0.866025388; 2., 0.; 3., 0.; 1., 0.; 0.5, 0.866025388; 0., 0.; -0.5, 0.866025388; -1., 0.];xcord=nodecord(:,1); % 模型的x坐标ycord=nodecord(:,2); % 模型的y坐标elementcon=[1, 2 ; 2, 3 ;3, 4 ;4, 2 ;1, 5 ;6, 1 ;6, 7 ;7, 5 ; 3, 1 ; 5, 3 ; 5, 6 ; 8, 6 ; 7, 8 ; 9, 7 ; 8, 9 ];Wd =2*node; % 模型的整体自由度    下面计算杆单元刚度矩阵和施加节点力。计算杆单元的全部代码见文末。% 计算杆单元刚度矩阵[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);% 施加荷载 (units: N)F=zeros(Wd,1);% 给第二个自由度施加45e3N的力，即给节点1的y方向施加45e3N的力F(2)=-45e3; % 给第十二个自由度施加45e3N的力，即给节点6的y方向施加45e3N的力F(12)=-45e3;    下面给节点4和9施加固定边界条件。%% 施加固定位移条件gdof=1: Wd;cdof=[7 8 17 18];    求解方程，这里是静力计算，实质上的计算量就是矩阵求逆。[U]=solve(Wd,gdof,cdof,S,F);displacements=[gdof' U]    下面是一些后处理的程序，以注释的形式给出解释说明。%% 统计变形后的节点坐标ux=1:2:2*node-1;uy=2:2:2*node;UX=U(ux);UY=U(uy);newcordx=xcord+UX;newcordy=ycord+UY;newcord=[newcordx,newcordy]%% 计算杆单元的力[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);A=1:size(p,1);Elemental_Force=[A' p]%% 绘制杆单元SC=100; %缩放因子，便以观察变形Snewcordx=xcord+SC*UX;Snewcordy=ycord+SC*UY;Snewcord=[Snewcordx,Snewcordy];drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);   运行上面的matlab程序，并将变形前后的桁架绘制在同一个窗口中观察变形，图形为   可以看到，变形的形状符合直觉，为了进一步验证程序的计算效果，采用abaqus创建同样的模型并对比计算结果，abaqus的模型为   abaqus计算完毕后，统计所有节点的竖直向位移与matlab程序对比，绘制表格为matlababaqus-0.003214286-0.00321429-0.001071429-0.00107143-0.002285714-0.002285710-4.50E-32-0.003428571-0.00342857-0.003214286-0.00321429-0.002285714-0.00228571-0.001071429-0.001071430-4.50E-32   可以看到，上面的计算结果一模一样，跟造数据了一样，事实上就算是我造数据了，你们也不知道，但是按我的习惯，我必然会把所的计算文件放在文末 。所有的计算代码   下面是matlab的计算主程序，主程序调用的函数我就不放在这里了，有需要的可以私聊我，私聊必发！clc ; clear ; % 清空变量%%% 杆单元材料参数E=210e3; % 弹性模量 (GPa=1000 N/mm^2)A=500; %杆单元截面积 (mm^2)node=9; % 模型节点总数element=15; % 模型单元总数nodecord=[ 1.5, 0.866025388; 2.5, 0.866025388; 2., 0.; 3., 0.; 1., 0.; 0.5, 0.866025388; 0., 0.; -0.5, 0.866025388; -1., 0.];xcord=nodecord(:,1); ycord=nodecord(:,2);elementcon=[1, 2 ; 2, 3 ;3, 4 ;4, 2 ;1, 5 ;6, 1 ;6, 7 ;7, 5 ; 3, 1 ; 5, 3 ; 5, 6 ; 8, 6 ; 7, 8 ; 9, 7 ; 8, 9 ];Wd =2*node; % Whole DOFs% 计算杆单元刚度矩阵[S]=planetruss(E,A, Wd,element,elementcon,xcord,ycord);% 施加节点力 (units: N)F=zeros(Wd,1);F(2)=-45e3;F(12)=-45e3;%% 施加固定位移边界条件gdof=1: Wd;cdof=[7 8 17 18];%% 求解位移[U]=solve(Wd,gdof,cdof,S,F);displacements=[gdof' U]%% 统计变形后的坐标ux=1:2:2*node-1;uy=2:2:2*node;UX=U(ux);UY=U(uy);newcordx=xcord+UX;newcordy=ycord+UY;newcord=[newcordx,newcordy]%% 计算单元力[p]=trussforce(E,A,element,elementcon,xcord,ycord,U);A=1:size(p,1);Elemental_Force=[A' p]%% 绘制变形前后的对比图SC=100; % 缩放因子，便于观察Snewcordx=xcord+SC*UX;Snewcordy=ycord+SC*UY;Snewcord=[Snewcordx,Snewcordy];drawplanetruss(element,elementcon,nodecord,newcord,Snewcord);    下面是abaqus的计算inp文件，这个inp文件拿过来就能直接算，算完的结果跟上面的matlab一样，曲线吻合的跟假的一样，非常能提升自己搞科研的自信心，建议保存起来，每当其他程序出bug的时候运行一下这个程序，以此欺骗和麻痹自己。*Heading** Job name: Job-1 Model name: Model-1** Generated by: Abaqus/CAE 6.14-2*Preprint, echo=NO, model=NO, history=NO, contact=NO**** PARTS***Part, name=Part-1*Node 1, 1.5, 0.866025388 2, 2.5, 0.866025388 3, 2., 0. 4, 3., 0. 5, 1., 0. 6, 0.5, 0.866025388 7, 0., 0. 8, -0.5, 0.866025388 9, -1., 0.*Element, type=T2D21, 1, 22, 2, 33, 3, 44, 4, 25, 1, 56, 6, 17, 6, 78, 7, 5 9, 3, 110, 5, 311, 5, 612, 8, 613, 7, 814, 9, 715, 8, 9*Nset, nset=Set-1 1, 5, 6, 7, 8*Elset, elset=Set-1 6, 7, 8, 11, 13*Nset, nset=Set-6, generate 1, 9, 1*Elset, elset=Set-6, generate 1, 15, 1*Nset, nset=fixed 4, 9*Nset, nset=middle 1, 6** Section: Section-1*Solid Section, elset=Set-6, material=Material-1500.,*End Part** **** ASSEMBLY***Assembly, name=Assembly** *Instance, name=Part-1-1, part=Part-1*End Instance** *End Assembly** ** MATERIALS** *Material, name=Material-1*Elastic210000.,0.** ** BOUNDARY CONDITIONS** ** Name: fixed Type: Displacement/Rotation*BoundaryPart-1-1.fixed, 1, 1Part-1-1.fixed, 2, 2** ----------------------------------------------------------------** ** STEP: Step-1** *Step, name=Step-1, nlgeom=NO*Static1., 1., 1e-05, 1.** ** LOADS** ** Name: middle Type: Concentrated force*CloadPart-1-1.middle, 2, -45000.** ** OUTPUT REQUESTS** *Restart, write, frequency=0** ** FIELD OUTPUT: F-Output-1** *Output, field, variable=PRESELECT** ** HISTORY OUTPUT: H-Output-1** *Output, history, variable=PRESELECT*End Step 来源：有限元先生