1杆系矩阵分析1

WRITE (11,'(2(A,I5))') & " 总共有",neq," 个方程、一维存储的元素总数量为 ",kdiag(neq) !-----------------------整体刚度矩阵集装-----------------------------------kv=zero elements_2: DO iel=1,nels CALL rod_km(km,prop(1,etype(iel)),ell(iel)); g=g_g(:,iel) CALL fsparv(kv,km,g,kdiag) END DO elements_2 !-----------------------读入荷载 和/或 位移--------------------------------loads=zero !-----------------------读入荷载信息 --------------------------------------READ (10,*) loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes) !-----------------------读入给定位移自由度数 ------------------------------READ (10,*) fixed_freedoms
矩阵位移法的复习和扩展
矩阵位移法复习 简单杆件程序及应用 自学相关程序
矩阵位移法
有限元法除单元分析（刚度方程建立）外整 理论基础——结构力学； 个思想方法与矩阵位移法是一样的，特别是 数学基础——矩阵分析； 从定位向量集装开始完全一样！这也就是应 计算手段——计算机程序。 该复习的原因。 理论上没有新内容，从方法上强调采用 矩阵表达，建立结构的位移法基本方程， 以计算机为工具进行结构分析。 • 要点：离散化、单元分析、集成总装、 整体结构基于位移法基本方程求解
• • • •
矩阵位移法复习
• 位移法基本思想 • 局部坐标单元分析 • 坐标转换 • 定位集装 • 解方程与后处理
位移法基本思想
• 通过加约束使之能拆成三类杆件的集合 • 用形常数和载常数得到单位荷载内力 • 附加约束反力为零建立集合体平衡条件 • 解方程得结点位移基本未知量 • 单位放大对应未知量倍和荷载内力叠加 通过“拆整为零，集零归整”解决问题
!-----------------------读入单元特性参数，本程序为轴向刚度 EA -------------READ (10,*) prop; etype=1 !-----------------------如果单元属性多于一类，读入单元属性类型向量 --------IF (np_types > 1) READ (10,*) etype !-----------------------读入单元长度 --------------------------------------READ (10,*) ell; nf=1 !-----------------------读入约束信息 --------------------------------------READ (10,*) nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf) ALLOCATE (kdiag(neq),loads(0:neq)) !-----------------------对单元循环，寻找整体数组尺寸-----------------------kdiag=0 elements_1: DO iel=1,nels; num=(/iel,iel+1/) CALL num_to_g(num,nf,g); g_g(:,iel)=g; CALL fkdiag(kdiag,g) END DO elements_1 DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO ALLOCATE (kv(kdiag(neq)))
IMPLICIT NONE INTEGER,PARAMETER :: iwp=SELECTED_REAL_KIND(15) ! 返回 iwp = 8，也 ! 即双精度 INTEGER :: fixed_freedoms,i,iel,k,loaded_nodes,ndof=2,nels,neq,nod=2, & nodof=1,nn,nprops=1,np_types,nr REAL(iwp) :: penalty=1.0e20_iwp,zero=0.0_iwp ! 大数，零 !-----------------------dynamic arrays- 动态数组 --------------------------INTEGER,ALLOCATABLE :: etype(:),g(:),g_g(:,:),kdiag(:),nf(:,:),no(:), & node(:),num(:) REAL(iwp),ALLOCATABLE :: action(:),eld(:),ell(:),km(:,:),kv(:),loads(:), & prop(:,:),value(:) !-----------------------input and initialisation- 输入和初始化 ------------OPEN (10,FILE=‘p41_1.dat’) ! 打开输入数据文件 OPEN (11,FILE=‘p41_1.res’) ! 打开输出数据文件 !-----------------------读入单元数及不同的材料类型数目 --------------------READ (10,*) nels,np_types nn=nels+1 ALLOCATE (g(ndof),num(nod),nf(nodof,nn),etype(nels),ell(nels),eld(ndof), & km(ndof,ndof),action(ndof),g_g(ndof,nels),prop(nprops,np_types))
DO k=1,nn; WRITE (11,'(I5,2E12.4)') k,loads(nf(:,k)); END DO !-----------------------retrieve element end actions- 计算输出单元杆端作用 WRITE (11,'(/A)') " 单元 起始端 终止端 " elements_3: DO iel=1,nels CALL rod_km(km,prop(1,etype(iel)),ell(iel)); g=g_g(:,iel) eld=loads(g); action=MATMUL(km,eld) WRITE (11,'(I5,2E12.4)') iel,action END DO elements_3 STOP END PROGRAM p41
PROGRAM p41 USE main ! 模块中必须包含如下的一些子程序，以便主程序应用 ! formnf 形成结点自由度向量 nf ! num_to_g 从 num 和 nf 找得 g 向量 num 单元结点号向量，nf 结点自由 ! 度向量，g 定位向量 ! fkdiag 计算一维存储主对角线元素位置信息 g 单元定位向量，kdiag 一 ! 维存储对角线元素位置向量 ! rod_km 形成一维杆单元刚度矩阵 km 单元刚度矩阵，ea 弹性属性向 ! 量,length 单元长度向量 ! fsparv 将单元刚度集装到对称一维存储总刚度矩阵 g 定位向量，kdiag ! 一维存储对角线元素位置向量，km 单元刚度矩阵，kv 一维存储 ! 总刚度矩阵 ! sparin 对对称一维总刚矩阵执行 Cholesky 因式分解 kv 一维存储总刚度 ! 矩阵，kdiag 一维存储对角线元素位置向量 ! spabac 对对称一维总刚矩阵执行 Cholesky 前向替换和回代
IF (fixed_freedoms /= 0) THEN ALLOCATE (node(fixed_freedoms), & no(fixed_freedoms),value(fixed_freedoms)) !-----------------------读入给定位移信息 ----------------------------------READ (10,*) (node(i),value(i),i=1,fixed_freedoms) DO i=1,fixed_freedoms; no(i)=nf(1,node(i)); END DO !------------------------ 进行乘大数处理 -----------------------------------kv(kdiag(no))=kv(kdiag(no))*penalty ! 乘大数 loads(no)=kv(kdiag(no))*value END IF !-----------------------equation solution - 方程求解 ----------------------CALL sparin(kv,kdiag) CALL spabac(kv,loads,kdiag) loads(0)=zero WRITE (11,'(/A)') " 结点 位移"
局部坐标单元分析
总的原则是叠加法
例如
δ3
δ2 δ1 δ6
q(x)
δ4
δ5
还记得如何写出单元刚度方程吗？
坐标转换
建立同一量两座标下分量间的关系
例如
y
δ2
y x
α δ1 δ1 x
δ2
能写出两分量之间的关系吗？
定位集装
在单元、结点编号基础上，再进 行位移编号 例如
7 ⑨ (10,11,12) ⑥ 4 ④ (1,2,3) ① 1 8 ⑩ (13,14,15) ⑦ 5 ⑤ (4,5,6) ② 2
能说出如何将单元的元素进行集装吗？
解方程与后处理
• 根据所采用的集装策略(方阵、等带宽、 一维变带宽等)进行集装，完成后视情况 进行边界条件处理，再调用对应的解法 • 解得结构结点位移后，根据定位向量形 成单元整体位移向量，计算单元杆端力 • 求单元任意截面的位移和内力以备作图 • 绘制结构的变形和内力图，需要的话还 可进一步做设计、验算等工作
简单杆件程序及应用
• 主程序介绍 • 输入数据文件建立 • 程序运行结果
! 程序 4.1 用 2 结点杆单元轴向荷载作用弹性杆的一维分析 ! 整型变量说明： ! fixed_freedoms = 给定位移自由度数，loaded_nodes = 受荷结点数 ! ndof = 每单元的自由度数，nels = 单元数，neq = 总自由度数（方程个 ! 数）, nod = 单元结点数，nodof =每结点自由度数，nn = 结点数 ! nprops = 材料（弹性特性）参数个数，np_types = 不同的材料类型数 ! 目，nr = 约束结点数 ! 整型数组说明： ! etype = 单元属性类型向量，g = 单元定位向量，g_g = 总单元定位向量 ! kdiag = 一维存储对角线元素位置向量，nf = 结点自由度向量 ! no = 给定位移的自由度号向量，node = 给定位移的点号向量，num = ! 单元结点号向量 ! 实型数组说明： ! action = 单元结点的荷载响应向量，eld = 单元位移向量，ell = 单元长 ! 度向量, km = 单元刚度矩阵，kv = 整体刚度矩阵，loads = 总荷载向量 ! prop = 单元轴向刚度 EA 矩阵，value = 给定的位移向量
9(16，17，18)
⑧ 6 (7，8，9) ③ 3(Biblioteka Baidu，0，0)
能写出各单元的定位向量吗？
定位集装
①的定位向量为（0,0,0,1,2,3） ④的定位向量为（1,2,3,4,5,6）
7 ⑥ 4 ① 1 ④ 5 ② 2 ⑨ 8 ⑦ ⑤ ⑩ 9(16，17，18) ⑧
例如
6 (7，8，9)
③ 3(0，0，0)