FE-Ch09-2杆系静力分析的有限元课件

SubroutineElem_Stiff(······)

说明Stiff=0.0!单元刚度清零SelectCase(Type)Case(1)

平面杆系结构单元Case(2)

空间杆系结构单元CaseDefault

出错信息EndSelectEndSubroutineElem_Stiff杆系结构单元分析子程序1.单元刚度总体设计2.说明部分设计Integer,Intent(in)::···

入口整型参数Real(8),Intent(in)::···

入口实型参数Real(8),Intent(out)::···

出口实型参数Real(8)::Work1,···Integer::i,j,k,···实型和整型工作变量SubroutineElem_Stiff(······)杆13.平面杆系结构设计SelectCase(Plane)Case(1)

平面桁架元素赋值Case(2)

平面梁柱元素赋值Case(3)······CaseDefault

出错信息EndSelect4.空间杆系结构设计SelectCase(Space)Case(1)

空间桁架元素赋值Case(2)

空间梁柱元素赋值Case(3)

交叉梁元素赋值CaseDefault

出错信息EndSelect3.平面杆系结构设计SelectCase(Pl25.有关单元等效结点荷载设计和进一步的考虑1)单元等效结点荷载设计同仿单元刚度。2)从各类单元刚度元素的计算，可看到要用到长度、单元弹性特性、单元截面特性等数据。因此，要确定存放它们的数据结构。要将它们作为出口。3)为计算单元等效结点荷载元素，首先要建立各种荷载情况等效荷载表达式，它们可由积分或载常数表得到。然后要解决荷载信息的存放结构，也要将它们作为出口量。4)单元刚度矩阵、等效结点荷载矩阵都应先清零。5.有关单元等效结点荷载设计和进一步的考虑1)单元等效结3杆系结构整体分析

首先就全刚结点平面刚架进行讨论，然后推广。1.总的思路在单元特性搞清后，将单元拼装回去。在结点处位移自动协调基础上，如果全部结点平衡，则求得的结点位移将是实际结构的解。因此，整体分析就是设法建立结点平衡方程。4.1.2坐标转换组成结构的杆件可以各个方向，单元分析对局部坐标，因此，必须将物理量转为统一坐标——整体坐标。1)力的转换关系杆系结构整体分析首先就全刚结点平面刚架进行讨论，然42)位移转换关系3)转换矩阵转换矩阵是正交矩阵。2)位移转换关系3)转换矩阵转换矩阵是正交矩阵54)杆端力转换5)杆端位移转换6)刚度方程的转换如果记称为整体单元刚度矩阵则这就是整体坐标下的单元刚度方程。本节以后的讨论认为都是对整体坐标的4)杆端力转换5)杆端位移转换6)刚度方程的转换63.结点平衡方程的建立1)一简单例子（如图）图中有两套编号，红的是单元杆端编号，黑的是结构整体编号。1-1)结点示意121221①②③图中蓝色的表示结点荷载（已知），红色的表示杆端力（未知的），、分别1、2单元杆端力子矩阵。对1、4结点“荷载”含有未知反力。21-2)结点平衡由示意图可见，结构结点的平衡方程为3.结点平衡方程的建立1)一简单例子（如图）图7从例图可见，其全部结点平衡方程为121221①②③若记2从例图可见，其全部结点121221①②③若记28式中[I]、[0]分别为单位和零矩阵。若引入矩阵记号，则结点平衡方程可改写作这一结论虽然是由一个例子得到的，但是显然对一切结构都是成立的。问题在于不同结构，[A]矩阵是不同的。式中[I]、[0]分别为单位和零矩阵。若引入矩阵记号，则结点94.杆端位移用结点位移来表示121221①②③仍以简单例子来说明若记由结点、杆端位移的协调条件，可得[]、[]的对应关系为式中[A]T是前面力关系[A]的转置，因此[A]T称为位移转换矩阵。4.杆端位移用结点位移来表示121221①②③仍以简单例子105.整体刚度方程——结点平衡121221①②③若记引入位移转换关系，则这就是整体刚度方程，它的物理实质是结点平衡。[K]称作结构刚度矩阵（或整体刚度矩阵），[R]称作综合等效结点荷载矩阵，它由两部分组成。5.整体刚度方程——结点平衡121221①②③若记11单元个数6.整体刚度矩阵的建立121221①②③若将[A]按单元分成图示三个子矩阵则由此可见，整体刚度矩阵可由各单元整体刚度矩阵装配累加得到。为说明如何装配，先将单元刚度矩阵进行分割整体结点码则由矩阵乘法可证明，[A]i[k]i[A]iT的结果是，将刚度矩阵子矩阵按整体结点码r、s送整体刚度矩阵相应位置。这一装配规则称为“对号入座”。单元个数6.整体刚度矩阵的建立121221①②③121)任意结构情况上面结论是通过具体例子（全刚结点平面刚架）得到的，由虚位移原理或势能原理进行整体分析（见讲义），可得任意结构其结论同此例。2)结点位移编号如果按结点顺序，对结点非零位移进行依次编号，这一序号称作结点位移码。为便于计算机处理并减少结构刚度矩阵的阶次，将零位移的号码变为零。对图示三铰刚架，当仅用一种单元（梁柱自由是单元）时结点位移编号如图所示。3)单元定位向量按单元局部结点码顺序，将结点位移码排成的向量，称作单元的定位向量。①②③④1)任意结构情况上面结论是通过具体例子（全刚结点13①②③④对图示刚架各单元的定位向量为①(0,0,1,3,4,5)②(0,0,2,10,11,12)③(3,4,5,6,7,8)④(6,7,9,10,11,12)①(0,0,1,2,3)②(0,0,,6,7,8)③(1,2,3,4,5)④(4,5,6,7,8)4)按单元定位向量集装刚度矩阵和综合荷载前面说明的是分块子矩阵集装,下面说明如何按定为向量来集装.如果如图所是采用各种不同的单元（一端有铰），则定位向量为①②③④①②③④如何获得带铰的单元刚度矩阵和等效荷载矩阵①②③④对图示刚架各单元的定位向量为①(014定位向量①②③④①②③④1)刚度集装以3单元为例来说明定位向量单元局部位移码根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素。定位向量①②③④①②③④1)刚度集装以3单元15根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。①②③④①②③④2)荷载集装以4单元为例来说明定位向量局部位移码此结论同样适用于刚度集装根据单元局部位移码和定位向量的对应关系用定①②③④①167整体分析总结1)对局部坐标和整体坐标不一致的单元，要对刚度、荷载进行坐标转换。2)需对“结构”进行结点、位移的局部和整体编号。4)集装所得整体刚度矩阵是对称、带状稀疏矩阵，当支撑条件能限制刚体位移时，矩阵非奇异。3)根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。据此可集装、累加得到整体刚度矩阵。7整体分析总结1)对局部坐标和整体坐标不一致的单元，要对175)综合荷载由两部分组成，因此首先要将直接作用结点的荷载按结点位移码送入，如果还有单元等效荷载，再按定位向量集装、累加。8)如果有某位移码方向弹性支撑，需进行将弹簧刚度送入位移码对应的对角线元素位置累加。9)如果有某位移码方向已知支撑位移，需进行将“边界条件处理”。具体做法以后介绍。7)整体刚度方程实质是全部结点的平衡条件。6)刚度矩阵带状稀疏，其带宽取决于结点、位移编码。最大半带宽=定位向量中最大元素差+1。5)综合荷载由两部分组成，因此首先要将直接作8)如果有182.5.8边界条件的处理10)当用虚位移或势能原理作整体分析时（势能为例），应变能为单元应变能之和，外力势能为单元外力势能之和+结点外力势能。全部杆端力的外力势能彼此抵消。1)乘大数法设第i个位移为已知值a，N=108或更大的数。乘大数法是将刚度矩阵Kii改为NKii，将Ri改为

Na。请考虑为什麽这样做能使边界条件得到满足？2)置换法（划零置1）设第i个位移为已知值a。2.5.8边界条件的处理10)当用虚位移或势能原理作整体19上述置换工作量大一些，显然可看出边界条件得到精确满足。上述置换工作量大一些，显然可看出边界条件得20★3)关于斜边界的处理如图示意的斜支座情况，有多种处理方案。3-1)通过单元的坐标转换来处理xy图示有斜支座单元，r结点处以倾角-来进行坐标转换，也即在r结点处整体坐标为图示xy。r3-2)通过增加一个单元来处理图示有斜支座单元，r结点处沿y方向增加一个刚结的单元，此单元有“无穷大”的抗拉刚度、但没有抗弯刚度。单元长度可任意。3-3)对整体刚度矩阵进行处理（参见教材）★3)关于斜边界的处理如图示意的斜支座情况，21最大半带宽10.总刚度矩阵的存储与对应解法因为总刚度矩阵对称、带状稀疏，利用这一特点可减少存储、提高解算效率。零元素零元素非零元素最大半带宽主对角线元素等半带存储零元素非零元素变带宽一维存储带宽是变的前面章节已经讨论最大10.总刚度矩阵的存储与对应解法因为总刚度矩22目前一般都用变带宽存储，下面结合程序说明存储和解法。首先介绍一些F90的语法。

定义导出类型导出类型——结点type::typ_Jointreal::x,y！坐标integer::GDOF(3)！整体位移码endtypetyp_Joint1)有关F90语法导出类型

新特性

目前一般都用变带宽存储，下面结合程序说明23用结点导出类型作为成员导出单元类型:type::typ_Element

integer::JointNo(2)！结点编号

type(typ_Joint)::Node(2)！结点

integer::GlbDOF(6)!定位向量

real::EA,EIendtypetyp_Elementtype(typ_Element)::Elem(5)！定义5个单元类型…!对单元i的端点j的x,y,GDOF(1:3)的赋值Elem(i)%Node(j)=Joint(Elem(i)%JointNo(j))…由导出类型定义新类型由导出类型定义变量用结点导出类型作为成员导出单元类型:type(typ_El24real::A(5),B(5,10),C(5)B=0.0!对B清零A=1.0!对A赋1:A(i)=1.0,i=1,5C=A+2!数组与标量运算:A(1:5)+(/2,2,2,2,2/)A=C+A!数组与数组运算（同形）C=sqrt(A)!数组的函数运算:C(i)=sqrt(A(i),i=1,5数组内部函数：dot_product(vector_a,vector_b)

!点积如：dot_product((/1,2,3/),(/2,3,4/))的值为20（待续）有关F90语法数组运算与赋值：real::A(5),B(5,10),C(5)有关F925matmul(matrix_a,matrix_b)

矩阵相乘如：locEDisp=matmul(T,glbEDisp)transpose(matrix)

矩阵转置如：glbEDisp=matmul(transpose(T),locEDisp)size(array,dim)

求数组第dim维的长度dim为可选变元:size(a,dim=2)若array为一维时，可不用dim。sum(array,dim,mask)数组元素求和dim,mask为可选变元；mask=条件表达式sum(a(1:10))对a的1到10元素求和sum(a(1:10),mask=a>0)对a(1:10)中大于0的元素求和（续）matmul(matrix_a,matrix_b)矩26有关F90语法where结构

新特性例where(C>0)C=0A=B*Dendwherewhere(C>0)A=Bendwhere定义where(数组关系表达式)

数组赋值语句

…elsewhere

数组赋值语句

...endwhere规则:

1)同形数组;2)不许嵌套;3)最多二分叉有关F90语法例定义规则:1)同形数组;2)不许嵌套27有关F90语法cycle和exit语句

新特性

用在do循环中

cycle——作下一个循环步

exit——跳出循环，执行enddo后一条语句等效例

do...if(.not.cond1)then...if(cond2)goto5...endifenddo5...用法do

...if(cond1)cycle ...

if(cond2)exit...enddo...有关F90语法等效例用法28有关F90语法数组构造函数spread语法spread(数组名,dim,ncopies)

将数组沿dim维方向复制ncopies形成新数组

dim,ncopies

—整型、位置变元、关键字变元（若按位置引用，可略关键字）例:（仅限一维数组）1）

spread(one,dim=1,ncopies=3)spread(one,1,3)spread(one,ncopies=3,dim=1)

[1,1,1]或[1,1,1]T有关F90语法语法例:（仅限一维数组）[1,1,1292）

ELocVec(1:6)=(/1,0,3,4,5,0/)spread(ELocVec,dim=1,ncopies=3)3）

spread(A(2,2:),dim=1,ncopies=2)如果dim=2呢?2）ELocVec(1:6)=(/1,0,3,4,5,030有关F90语法指针pointerpointer是变量的属性，可以指向相同类型的变量;被指向的目标必须具有target属性或pointer属性可以将指针变量理解为别名、称号real,target::a,b,EDisp(6)!可被指针所指

real,pointer::p1,p2!称号：班长、课代表！p1,p2是指针，可以指向实型数据real,pointer::G(:)!先进集体！G是指针，可以指向一维实型数组有关F90语法pointer是变量的属性，可以指向相同类31指针是一种“称号”，上述声明语句建立了“称号”，但并未“授予”某个变量这个称号，因此是指向“空”，并未占用内存a=3.0p1=>a!p1指向a!称号p1授予a，a的数据有两个名：固定名a和流动名p1;既可用p1也可用a（p1—班长，a—张三）a=4.0!a的值变为4.0,p1也变为4.0p1=>b!班长换人了G=>EDisp!先进集体有了得主！EDisp(:)的长度就是G(:)的长度，用G和用EDisp同样效果指针是一种“称号”，上述声明语句建立了“称号”，但并未“授予32又如：real,target::a,breal,pointer::p,qa=3.14b=2.0p=>a!p=a=3.14q=>b!q=b=2.0q=p!(q指向的数据b)=(p指向的数据a)！即：所有=3.14又如：33指针可以指向有名的数据区，也可以指向无名的数据区real,pointer::b(:,:)integer::nread(*,*)nallocate(b(n,n))

!指针指向了一个刚开辟的数组!以下可以当作常规数组用b(1,1)=1.1b(1,2)=1.2...deallocate(b)有关F90语法用指针建立动态数组指针可以指向有名的数据区，也可以指向无名的数据区有关F9034指针与allocatable数组的区别具备allocatable数组的所有功能还可以用在导出类型中，例如整体刚度矩阵的变带宽存储：type::typ_Kcol

!整体刚度矩阵K的列

real(8),pointer::row(:)

!该列的行元素endtype...type(typ_Kcol),allocatable::Kcol(:)…allocate(Kcol(NGlbDOF))!分配了NGlbDOF列...allocate(Kcol(5)%row(3:5))!第5列只用3至5行指针与allocatable数组的区别35(1)NElem,NJoint,NGlbDOF,NJLoad,NELoad单元数结点数总自由度数结点荷载数单元荷载数[Joint-结点]…NJoint行(2)Joint%X,Joint%Y,GDOF(1:3)X坐标Y坐标结点位移码[Elem-单元]…NElem行(3)JointNo1,JointNo2,EA,EI起点号终点号刚度[JLoad-结点荷载]…NJLoad行(4)JointNo,LodDOF,LodVal作用点号局部位移码荷载值[ELoad-单元荷载]…NELoad行(5)ElemNo,Indx,Pos,LodVal单元号类型号位置荷载值Indx

类型

pos1均布长度2集中位置3...2)某程序输入数据说明(1)NElem,NJoint,NGlbDOF,363,5,7,1,10,0,0,0,00,4,1,2,34,4,4,5,64,4,4,5,74,0,0,0,01,2,1.0E9,162,3,1.0E9,245,4,1.0E9,122,1,10.0E3

!结点2，自由度1，值为10E32,1,4,-4.0E3

！单元2，均布，长4米，值为-4E32-1)数据文件例子：(2)(1)(3)24135i=6i=4i=310kN4kN/m4m4mEA=109N(1)(2)坐标位移码(3)(4)(5)结点号EA,EI3,5,7,1,12-1)数据文件例子：(2)(1)(3)37read(5,*)NElem,NJoint,NGlbDOF,NJLoad,NELoadallocate(Joint(NJoint))allocate(Elem(NElem))allocate(JLoad(NJLoad))allocate(ELoad(NELoad))...read(5,*)(Joint(i),i=1,NJoint)read(5,*)(Elem(ie)%JointNo,Elem(ie)%EA,&ELem(ie)%EI,ie=1,NElem)if(NJLoad>0)read(5,*)(JLoad(i),i=1,NJLoad)if(NELoad>0)read(5,*)(ELoad(i),i=1,NELoad)2-2)程序读入计算所需数据：read(5,*)NElem,NJoint,NGlbDOF382-3)求始行码和分配带宽子程序!==================================subroutineSetMatBand(Kcol,Elem)!接口简单!==================================

type(typ_Kcol),intent(inout)::Kcol(:)！总刚列type(typ_Element),intent(in))::Elem(:)！单元integer(ikind)::minDOF,ELocVec(6)integer(ikind)::Row1(size(Kcol,dim=1))

！Row1为自动数组，子程序结束后自动释放。！这样做可使接口简单，减少了数组的控制变量。

integer(ikind)::ie,j,NGlbDOF,NElemNGlbDOF=size(Kcol,dim=1)！使接口简单2-3)求始行码和分配带宽子程序!===========39NElem=size(Elem,dim=1)Row1=NGlbDOF

!先设所有始行码同终行码!确定各列始行码向量

doie=1,Nelem!对单元循环!确定定位向量ELocVec(:)=Elem(ie)%GlbDOF(:)!寻找定位向量中大于零的最小值minDOF=minval(ELocVec,mask=ELocVec>0)

!屏蔽定位向量中小于等于零的where(ELocVec>0)

!寻找Row1(ELocVec)和minDOF中的最小值Row1(ELocVec)=min(Row1(ELocVec),&minDOF)endwhereenddoNElem=size(Elem,dim=1)40!为各列的带宽分配空间

doj=1,NGlbDOF!对总自由度数循环allocate(Kcol(j)%row(Row1(j):j)

!给Kcol(j)%row从Row1(j)到j个空间Kcol(j)%row=zero

!总刚元素清零

enddoreturnendsubroutineSetMatBand!为各列的带宽分配空间413)整体刚度矩阵的集成

doie=1,NElem…!计算单刚

EK(6,6),ELocVec(6)doj=1,6!对单元逐列集成JGDOF=ELocVec(j)!取出位移码if(JGDOF=0)cycle!作下一循环步where(ELocVec>0.And.ELocVec<=JGDOF)

!位移码非零同时小于第j个局部码对应的位移码Kcol(JGDOF)%row(ELocVec)=&Kcol(JGDOF)%row(ELocVec)+EK(:,j)!集成endwhereenddoenddo局部码位移码3)整体刚度矩阵的集成doie=1,NElem局部42变带宽矩阵的分解求解1)LDLT分解法求解[A]{x}

={b}若对称正定，则可分解为=单位上三角阵对角阵原方程变为求解步骤：1.分解：2.解y：3.解x：LDLT分解法Gauss消去前消去处理右端项回代(不同的b只做一次分解)变带宽矩阵的分解求解若对称正定，则可分解为=单位上三角阵对角43存放：主对角—

上三角—

存放：主对角—上三角—44时不必求和(上三角:i<j)不动存在处存在处时不必求和(上三角:i<j)不动存在处存在45(第j列系数)上三角:i<j2)分解一般情况：逐列分解对角线:i=

j>1(第j列系数)上三角:i<j2)分解一般情况：46(第j列系数)第

i列中第1个非零元素的行码为：3)变列宽存贮的分解修正：

第

j

列中第1个非零元素的行码为：

则分解顺序(第j列系数)第i列中第1个非零元素的行码为：3)474)F90实现!三角分解diag(1:ncol)=(/(Kcol(j)%row(j),j=1,ncol)/)doj=2,ncolrow1=lbound(Kcol(j)%row,1)!i_1jdoi=row1,j-1row_1=max(row1,lbound(Kcol(i)%row,1))!i_1s=sum(diag(row_1:i-1)*Kcol(i)%row(row_1:i-1)*&Kcol(j)%row(row_1:i-1))

!求和部分

Kcol(j)%row(i)=(Kcol(j)%row(i)-s)/diag(i)enddo

！第j列系数分解完