杆系静力分析课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Subroutine,Elem_Stiff,(),说明,Stiff=0.0 !,单元刚度清零,Select Case (Type),Case (1),平面杆系结构单元,Case (2),空间杆系结构单元,Case Default,出错信息,End Select,End Subroutine,Elem_Stiff,3.4,杆系结构单元分析子程序,3.4.1,单元刚度总体设计,3.4.2,说明部分设计,Integer,Intent (in) : ,入口整型参数,Real(8),Intent(in) : ,入口实型参数,Real(8),Intent(out) : ,出口实型参数,Real(8) : Work1, ,Integer : i,j,k, ,实型和整型工作变量,3.4,杆系结构单元分析子程序,3.4.3,平面杆系结构设计,Select Case (Plane),Case (1),平面桁架元素赋值,Case (2),平面梁柱元素赋值,Case (3), ,Case Default,出错信息,End Select,3.4.4,空间杆系结构设计,Select Case (Space),Case (1),空间桁架元素赋值,Case (2),空间梁柱元素赋值,Case (3),交叉梁元素赋值,Case Default,出错信息,End Select,3.4,杆系结构单元分析子程序,3.4.5,有关单元等效结点荷载设计和进一步的考虑,1),单元等效结点荷载设计同仿单元刚度。,2),从各类单元刚度元素的计算，可看到要用到长度、,单元弹性特性、单元截面特性等数据。因此，要,确定存放它们的数据结构。要将它们作为出口。,3),为计算单元等效结点荷载元素，首先要建立各种,荷载情况等效荷载表达式，它们可由积分或载常,数表得到。然后要解决荷载信息的存放结构，也,要将它们作为出口量。,4),单元刚度矩阵、等效结点荷载矩阵都应先清零。,4.1,杆系结构整体分析,首先就全刚结点平面刚架进行讨论，然后推广。,4.1.1,总的思路,在单元特性搞清后，将单元拼装回去。在结点处,位移自动协调基础上，如果全部结点平衡，则求得,的结点位移将是实际结构的解。因此，整体分析就,是设法建立结点平衡方程。,4.1,.2,坐标转换,组成结构的杆件可以各个,方向，单元分析对局部坐标，,因此，必须将物理量转为统,一坐标,整体坐标。,1),力的转换关系,4.1,杆系结构整体分析,2),位移转换关系,3),转换矩阵,转换矩阵是正交矩阵。,4.1,杆系结构整体分析,4),杆端力转换,5),杆端位移转换,6),刚度方程的转换,如果记 称为整体单元刚度矩阵,则,这就是整体坐标下的单元刚度方程。,本节以后的讨论认为,都是对整体坐标的,4.1,杆系结构整体分析,4.1.3,结点平衡方程的建立,1),一简单例子（如图）,图中有两套编号，,红,的,是单元杆端编号，,黑,的是,结构整体编号。,1-1),结点示意,1,2,1,2,2,1,图中蓝色的表示结点荷载（已知），,红,色的表示,杆端力（未知的）， 、 分别,1,、,2,单元杆端力,子矩阵。对,1,、,4,结点,“荷载”,含有未知反力。,2,1-2),结点平衡,由示意图可见，结构,结点的平衡方程为,4.1,杆系结构整体分析,从例图可见，其全部结点,平衡方程为,1,2,1,2,2,1,若记,2,4.1,杆系结构整体分析,式中,I,、,0,分别为单位和零矩阵。,若引入矩阵记号,，则结点平衡方程可改写作,这一结论虽然是由一个例子得到的，但是显然,对一切结构都是成立的。问题在于不同结构,，,A,矩阵是不同的。,4.1,杆系结构整体分析,4.1.4,杆端位移用结点位移来表示,1,2,1,2,2,1,仍以简单例子来说明,若记,由结点、杆端位移的协调条件，可得,、,的对应关系为,式中,A,T,是前面力关系,A,的转置，因此,A,T,称,为,位移转换矩阵,。,4.1,杆系结构整体分析,4.1.5,整体刚度方程,结点平衡,1,2,1,2,2,1,若记,引入位移转换关系，则,这就是整体刚度方程，它的物理实质是结点平衡。,K,称作结构刚度矩阵（或整体刚度矩阵），,R,称作综合等效结点荷载矩阵，它由两部分组成。,单元个数,4.1,杆系结构整体分析,4.1.6,整体刚度矩阵的建立,1,2,1,2,2,1,若将,A,按单元分成图示,三个子矩阵,则,由此可见，整体刚度矩阵可由各单元整体刚度矩阵,装配累加,得到。为说明如何装配，先将单元刚度矩阵进行分割,整体结点码,则由矩阵乘法可证明，,A,i,k,i,A,i,T,的结果是，将,刚度矩阵子矩阵按整体结点码,r,、,s,送整体刚度矩阵,相应位置。,这一装配规则称为“对号入座”。,4.1,杆系结构整体分析,1),任意结构情况,上面结论是通过具体例子（全刚结点平面刚架）得到的，由虚位移原理或势能原理进行整体分析（见讲义），可得任意结构其结论同此例。,2),结点位移编号,如果按结点顺序，对结点非零位移进行依次编号，这一序号称作结点位移码。为便于计算机处理并减少结构刚度矩阵的阶次，将零位移的号码变为零。,对图示三铰刚架，当仅,用一种单元（梁柱自由是,单元）时结点位移编号如,图所示。,3),单元定位向量,按单元局部结点码顺序，,将结点位移码排成的向量，称作单元的,定位向量,。,4.1,杆系结构整体分析,对图示刚架各单元的定位向量为,(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),按单元定位向量集装刚度矩阵和综合荷载,前面说明的是分块子矩阵集装,下面说明如何按定,为向量来集装,.,如果如图所是采用各种不同的单元,（一端有铰,），,则定位向量为,如何获得,带铰的单元刚,度矩阵和等,效荷载矩阵,定位向量,4.1,杆系结构整体分析,4-1),刚度集装,以,3,单元为例来说明,定位向量,单元局部位移码,根据单元,局部位移码,和定位向量,的对应关系,用定位向量,位移码送元,素。,根据单元局部位移码和定位向量的对应关系用定,位向量位移码送元素，定位向量元素为零时不送。,4.1,杆系结构整体分析,4-2),荷载集装,以,4,单元为例来说明,定位向量,局,部,位,移,码,此结论同样适,用于刚度集装,4.1,杆系结构整体分析,4.1.7,整体分析总结,1),对局部坐标和整体坐标不一致的单元，要对刚,度、荷载进行坐标转换。,2),需对“结构”进行结点、位移的局部和整体编,号,。,4),集装所得整体刚度矩阵是对称、带状稀疏矩阵，,当支撑条件能限制刚体位移时，矩阵非奇异。,3),根据单元局部位移码和定位向量的对应关系用,定位向量位移码送元素，定位向量元素为零时,不送。据此可集装、累加得到整体刚度矩阵。,5),综合荷载由两部分组成，因此首先要将直接作,用结点的荷载按结点位移码送入，如果还有单,元等效荷载，再按定位向量集装、累加。,4.1,杆系结构整体分析,8),如果有某位移码方向弹性支撑，需进行将弹簧,刚度送入位移码对应的对角线元素位置累加。,9),如果有某位移码方向已知支撑位移，需进行将,“边界条件处理”。,具体做法以后介绍。,7),整体刚度方程实质是全部结点的平衡条件。,6),刚度矩阵带状稀疏，其带宽取决于结点、位移,编码。,最大半带宽,=,定位向量中最大元素差,+1,。,4.1,杆系结构整体分析,2.5.8,边界条件的处理,10),当用虚位移或势能原理作整体分析时（势能为,例），应变能为单元应变能之和，外力势能为,单元外力势能之和,+,结点外力势能。全部杆端力,的外力势能彼此抵消。,1),乘大数法,设第,i,个位移为已知值,a,，,N,=10,8,或更大的数。,乘大数法是,将刚度矩阵,K,ii,改为,N,K,ii,，将,R,i,改为,N,a,。,请考虑,为什麽这样做能使,边界条件得到满足？,2),置换法（划零置,1,）,设第,i,个位移为已知值,a,。,4.1,杆系结构整体分析,上述置换工作量大一些，显然可看出边界条件得,到精确满足。,4.1,杆系结构整体分析, 3,),关于斜边界的处理,如图示意的斜支座情况，有多种处理方案。,3-1),通过单元的坐标转换来处理,x,y,图示有斜支座单元，,r,结点处,以倾角,- ,来进行坐标转换，也,即在,r,结点处整体坐标为图示,xy,。,r,3-2),通过增加一个单元来处理,图示有斜支座单元，,r,结点处沿,y,方向增加一个,刚结的单元，,此单元,有,“无穷大”的抗拉刚度、,但,没有抗弯刚度。单元长度可任意。,3-3),对整体刚度矩阵进行处理（参见匡文起教材）,最大,半带宽,2.5,杆系结构整体分析,2.5.10 总刚度矩阵的存储与对应解法,因为总刚度矩阵对称、带状稀疏，利用这一特,点可减少存储、提高解算效率。,零元素,零元素,非零元素,最大,半带宽,主对角线元素,等,半,带,存,储,零元素,非,零,元,素,变,带,宽,一,维,存,储,带,宽,是,变,的,到P.55,2.5,杆系结构整体分析,目前一般都用变带宽存储，下面结合程序说明,存储和解法。首先介绍一些,F 90,的语法。,定义导出类型,导出类型结点,type :,typ,_Joint,real :,x,y,！,坐标,integer,: GDOF(3),！,整体位移码,end type,typ,_Joint,1),有关,F90,语法,导出类型,新特性,用结点导出类型作为成员导出单元类型:,type :,typ,_Element,integer,:,JointNo,(2),！,结点编号,type(,typ,_Joint) : Node(2),！,结点,integer :,GlbDOF,(6) !,定位向量,real,: EA,EI,end type,typ,_Element,type (,typ,_Element),:,Elem,(5),！,定义5个单元类型,! 对单元,i,的端点,j,的,x,y,GDOF(1:3),的赋值,Elem,(i)%Node(j) = Joint(,Elem,(i)%,JointNo,(j),由导出类型定义新类型,由导出类型定义变量,real : A(5),B(5,10),C(5),B=0.0,!,对,B,清零,A=1.0,!,对,A,赋1:,A(i)=1.0, i=1,5,C=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,语法,数组运算与赋值,：,matmul,(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=a0),对,a(1:10),中大于0的元素求和,（续）,有关,F90,语法,where,结构,新特性,例,where (C0),C = 0,A = B*D,end where,where (C0),A = B,end where,定义,where,(,数组关系表达式,),数组赋值语句,elsewhere,数组赋值语句,.,end where,规则:,1)同形数组; 2)不许嵌套; 3)最多二分叉,有关,F90,语法,cycle,和,exit,语句,新特性,用在,do,循环中,cycle,作下一个循环步,exit,跳出循环，执行,end do,后一条语句,等效例,do,.,if(.not.cond1)then,.,if(cond2),goto,5,.,end if,end do,5 .,用法,do,.,if (cond1) cycle,.,if (cond2) exit,.,end do,.,有关,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,2）,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,呢?,有关,F90,语法,指针,pointer,pointer,是变量的属性，可以指向相同类型的变量;,被指向的目标必须具有,target,属性或,pointer,属性,可以将指针变量理解为别名、称号,real,target : a,b,EDisp,(6),!,可被指针所指,real,pointer : p1,p2,!,称号：班长、课代表,！,p1,p2,是指针，可以指向实型数据,real,pointer : G(:),!,先进集体,！,G,是指针，可以指向一维实型数组,指针是一种,“,称号,”,，上述声明语句建立了,“,称号,”,，但并未,“,授予,”,某个变量这个称号，因此是指向,“,空,”,，并未占用内存,a = 3.0,p1 = a,! p1,指向,a,!,称号,p1,授予,a,，,a,的数据有两个名：固定名,a,和流动名,p1;,既可用,p1,也可用,a（p1,班长,，,a,张三,）,a = 4.0,! a,的值变为,4.0,p1,也变为,4.0,p1 = b,!,班长换人了,G =,EDisp,!,先进集体有了得主,！,EDisp,(:),的长度就是,G(:),的长度，用,G,和用,EDisp,同样效果,又如：,real,target : a,b,real,pointer: p,q,a = 3.14,b = 2.0,p = a,! p = a = 3.14,q = b,! q = b = 2.0,q = p,! (,q,指向的数据,b,) = (,p,指向的数据,a ),！,即： 所有,= 3.14,指针可以指向有名的数据区，也可以指向无名的数据区,real,pointer: b(:,:),integer : n,read,(*,*),n,allocate (b(n,n),!,指针指向了一个刚开辟的数组,!,以下可以当作常规数组用,b(1,1) = 1.1,b(1,2) = 1.2,.,d,eallocate,(b),有关,F90,语法,用指针建立动态数组,指针与,allocatable,数组的区别,具备,allocatable,数组的所有功能,还可以用在导出类型中，例如整体刚度矩阵的变带宽存储：,type :,typ,_,Kcol,!,整体刚度矩阵,K,的列,real(8),pointer : row(:),!,该列的行元素,end type,.,type (,typ,_,Kcol,),allocatable,:,Kcol,(:),allocate (,Kcol,(,NGlbDOF,),!,分配了,NGlbDOF,列,.,allocate (,Kcol,(5)%row(3:5),!,第5列只用3至5行,(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,类型,pos,1,均布 长度,2 集中 位置,3 .,2) 某程序输入数据说明,3,5,7,1,1,0,0, 0,0,0,0,4, 1,2,3,4,4, 4,5,6,4,4, 4,5,7,4,0, 0,0,0,1,2, 1.0E9, 16,2,3, 1.0E9, 24,5,4, 1.0E9, 12,2,1,10.0E3,!,结点,2,，自由度,1,，值为,10E3,2,1,4,-4.0E3,！,单元,2,，均布，长,4,米，值为,-4E3,2-1),数据文件例子：,(2),(1),(3),2,4,1,3,5,i,= 6,i,= 4,i,= 3,10,kN,4,kN/m,4m,4m,EA,= 10,9,N,(1),(2),坐标,位移码,(3),(4),(5),结点号,EA,EI,read(5,*),NElem,NJoint,NGlbDOF,NJLoad,NELoad,allocate (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) 程序读入计算所需数据：,2-3) 求始行码和分配带宽子程序,!=,subroutine,SetMatBand,(,Kcol,Elem,),!,接口简单,!=,type(,typ,_,Kcol,),intent(in out) :,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,NElem,NGlbDOF,= size(,Kcol,dim=1),！,使接口简单,NElem,= size(,Elem,dim=1),Row1 =,NGlbDOF,!,先设所有始行码同终行码,!,确定各列始行码向量,do,ie,=1,Nelem,!,对单元循环,! 确定定位向量,ELocVec,(:) =,Elem,(,ie,)%,GlbDOF,(:),!,寻找定位向量中大于零的最小值,minDOF,=,minval,(,ELocVec,mask=,ELocVec,0),!,屏蔽定位向量中小于等于零的,where (,ELocVec, 0),!,寻找,Row1(,ELocVec,),和,minDOF,中的最小值,Row1(,ELocVec,) = min(Row1(,ELocVec,), &,minDOF,),end where,end do,!,为各列的带宽分配空间,do j=1,NGlbDOF,!,对总自由度数循环,allocate (,Kcol,(j)%row(Row1(j):j),!,给,Kcol,(j)%row,从,Row1(j),到,j,个空间,Kcol,(j)%row = zero,!,总刚元素清零,end do,return,end subroutine,SetMatBand,3) 整体刚度矩阵的集成,do,ie,=1,NElem,!,计算单刚,EK(6,6),ELocVec,(6),do j=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),!,集成,end where,end do,end do,局部码,位移码,4),变带宽矩阵的分解求解,4-1),LDL,T,分解法求解 ,A,x,= ,b,若,对称正定，则可分解为,=,单位上三角阵,对角阵,原方程变为,求解步骤：,1. 分解：,2. 解,y,：,3. 解,x,：,LDL,T,分解法,Gauss,消去,前消去,处理右端项,回代,(不同的,b,只做一次分解),存放：,主对角,上三角 ,时不必求和,(,上三角:,i,j,),不动,存在 处,存在 处,(第,j,列系数),上三角:,i,1,(第,j,列系数),第,i,列中第1个非零元素的行码为：,4-3) 变列宽存贮的分解修正：,第,j,列中第1个非零元素的行码为：,则,分解顺序,4-4),F90,实现,!,三角分解,diag,(1:,ncol,) = (/ (,Kcol,(j)%row(j),j=1,ncol,) /),do j = 2,ncol,row1 =,lbound,(,Kcol,(j)%row,1),! i_1 j,do i = row1,j-1,row_1 = max(row1,lbound,(,Kcol,(i)%row,1),! i_1,s = 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),end do,！,第,j,列系数分解完毕,s = sum(,diag,(row1:j-1)*,Kcol,(j)%row(row1:j-1,)*,2),diag,(j) =,diag,(j) - s,!,第,j,列的主对角元素,end do,4-5),一般情况解,y,公式,可不动,可存在 处,4-6) 变列宽存储解,y,的修正,上三角的第,i,列,从第 行元素开始,!,回代步骤1: 自上而下,do i = 2,ncol,row1 =,lbound,(,Kcol,(i)%row,dim = 1),GP(i) = GP(i) - sum(,Kcol,(i)%row(row1:i-1) &,* GP(row1:i-1) ),end do,求出后， 不再用，可将 存在 处,4-7) 一般解,x,的,公式,自下而上：,这样做的缺欠：,自下而上逐行计算，行中遇到0元素需跳过，不方便！,4-8) 改为自右向左逐列计算！,记,解出,对,y,向量预处理,再解出,第,n,列乘上 后移到,右边去，修正,y,向量,!,回代步骤2: 自右向左,GP(:) = GP(:)/,diag,(:),do j =,ncol,2,-1,row1 =,lbound,(,Kcol,(j)%row,dim = 1),GP(row1:j-1) = GP(row1:j-1) - GP(j) * &,Kcol,(j)%row(row1:j-1),end do,!,现在,GP,就是解,4-9) 小结：,无须一维存贮，无须行列码转换定位,公式与程序直接对应翻译，直截了当,求和采用内部函数进行数组运算,动态内存，用多大、开多大,数据封装性好，接口简单:,subroutine,SetMatBand,(,Kcol,Elem,),subroutine,VarBandSolv,(,Disp,Kcol,GP),通用性强,2.6,杆系结构内力计算,2.6.1 杆端内力计算公式,解方程的结果可以得到结点位移，有了位移就,可以进一步求各单元的内力。,解算步骤：,从,根据定位向量取出单元杆端位移。,由单元倾角确定是算 还是 。,减去等效结点荷载或加上固端力矩阵。,2.6,杆系结构内力计算,2.6.2,杆中任意截面内力计算公式,注意事项：,由图示隔离体图，列杆段的平衡方程即可得到,任意截面的内力计算公式。,请自行写出。,按上述隔离体图所求得的内力是精确的。,求得单元结点位移后，代入单元位移场、求应,变、求内力，这样所得的内力一般不满足平衡条,件。只是近似解。,