桿系靜力分析課件

SubroutineElem_Stiff(······)

說明Stiff=0.0!單元剛度清零SelectCase(Type)Case(1)

平面桿系結(jié)構(gòu)單元Case(2)

空間桿系結(jié)構(gòu)單元CaseDefault

出錯信息EndSelectEndSubroutineElem_Stiff3.4桿系結(jié)構(gòu)單元分析子程序3.4.1單元剛度總體設(shè)計3.4.2說明部分設(shè)計Integer,Intent(in)::···

入口整型參數(shù)Real(8),Intent(in)::···

入口實(shí)型參數(shù)Real(8),Intent(out)::···

出口實(shí)型參數(shù)Real(8)::Work1,···Integer::i,j,k,···實(shí)型和整型工作變量SubroutineElem_Stiff(······)313.4桿系結(jié)構(gòu)單元分析子程序3.4.3平面桿系結(jié)構(gòu)設(shè)計SelectCase(Plane)Case(1)

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

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

出錯信息EndSelect3.4.4空間桿系結(jié)構(gòu)設(shè)計SelectCase(Space)Case(1)

空間桁架元素賦值Case(2)

空間梁柱元素賦值Case(3)

交叉梁元素賦值CaseDefault

出錯信息EndSelect3.4桿系結(jié)構(gòu)單元分析子程序3.4.3平面桿系結(jié)構(gòu)設(shè)計23.4桿系結(jié)構(gòu)單元分析子程序3.4.5有關(guān)單元等效結(jié)點(diǎn)荷載設(shè)計和進(jìn)一步的考慮1)單元等效結(jié)點(diǎn)荷載設(shè)計同仿單元剛度。2)從各類單元剛度元素的計算，可看到要用到長度、單元彈性特性、單元截面特性等數(shù)據(jù)。因此，要確定存放它們的數(shù)據(jù)結(jié)構(gòu)。要將它們作為出口。3)為計算單元等效結(jié)點(diǎn)荷載元素，首先要建立各種荷載情況等效荷載表達(dá)式，它們可由積分或載常數(shù)表得到。然后要解決荷載信息的存放結(jié)構(gòu)，也要將它們作為出口量。4)單元剛度矩陣、等效結(jié)點(diǎn)荷載矩陣都應(yīng)先清零。3.4桿系結(jié)構(gòu)單元分析子程序3.4.5有關(guān)單元等效結(jié)點(diǎn)荷34.1桿系結(jié)構(gòu)整體分析

首先就全剛結(jié)點(diǎn)平面剛架進(jìn)行討論，然后推廣。4.1.1總的思路在單元特性搞清后，將單元拼裝回去。在結(jié)點(diǎn)處位移自動協(xié)調(diào)基礎(chǔ)上，如果全部結(jié)點(diǎn)平衡，則求得的結(jié)點(diǎn)位移將是實(shí)際結(jié)構(gòu)的解。因此，整體分析就是設(shè)法建立結(jié)點(diǎn)平衡方程。4.1.2坐標(biāo)轉(zhuǎn)換組成結(jié)構(gòu)的桿件可以各個方向，單元分析對局部坐標(biāo)，因此，必須將物理量轉(zhuǎn)為統(tǒng)一坐標(biāo)——整體坐標(biāo)。1)力的轉(zhuǎn)換關(guān)系4.1桿系結(jié)構(gòu)整體分析首先就全剛結(jié)點(diǎn)平面剛架進(jìn)行44.1桿系結(jié)構(gòu)整體分析2)位移轉(zhuǎn)換關(guān)系3)轉(zhuǎn)換矩陣轉(zhuǎn)換矩陣是正交矩陣。4.1桿系結(jié)構(gòu)整體分析2)位移轉(zhuǎn)換關(guān)系3)轉(zhuǎn)換矩陣54.1桿系結(jié)構(gòu)整體分析4)桿端力轉(zhuǎn)換5)桿端位移轉(zhuǎn)換6)剛度方程的轉(zhuǎn)換如果記稱為整體單元剛度矩陣則這就是整體坐標(biāo)下的單元剛度方程。本節(jié)以后的討論認(rèn)為都是對整體坐標(biāo)的4.1桿系結(jié)構(gòu)整體分析4)桿端力轉(zhuǎn)換5)桿端位移轉(zhuǎn)換664.1桿系結(jié)構(gòu)整體分析4.1.3結(jié)點(diǎn)平衡方程的建立1)一簡單例子（如圖）圖中有兩套編號，紅的是單元桿端編號，黑的是結(jié)構(gòu)整體編號。1-1)結(jié)點(diǎn)示意121221①②③圖中藍(lán)色的表示結(jié)點(diǎn)荷載（已知），紅色的表示桿端力（未知的），、分別1、2單元桿端力子矩陣。對1、4結(jié)點(diǎn)“荷載”含有未知反力。21-2)結(jié)點(diǎn)平衡由示意圖可見，結(jié)構(gòu)結(jié)點(diǎn)的平衡方程為4.1桿系結(jié)構(gòu)整體分析4.1.3結(jié)點(diǎn)平衡方程的建立1)74.1桿系結(jié)構(gòu)整體分析從例圖可見，其全部結(jié)點(diǎn)平衡方程為121221①②③若記24.1桿系結(jié)構(gòu)整體分析從例圖可見，其全部結(jié)點(diǎn)1284.1桿系結(jié)構(gòu)整體分析式中[I]、[0]分別為單位和零矩陣。若引入矩陣記號，則結(jié)點(diǎn)平衡方程可改寫作這一結(jié)論雖然是由一個例子得到的，但是顯然對一切結(jié)構(gòu)都是成立的。問題在于不同結(jié)構(gòu)，[A]矩陣是不同的。4.1桿系結(jié)構(gòu)整體分析式中[I]、[0]分別為單位和零矩陣94.1桿系結(jié)構(gòu)整體分析4.1.4桿端位移用結(jié)點(diǎn)位移來表示121221①②③仍以簡單例子來說明若記由結(jié)點(diǎn)、桿端位移的協(xié)調(diào)條件，可得[]、[]的對應(yīng)關(guān)系為式中[A]T是前面力關(guān)系[A]的轉(zhuǎn)置，因此[A]T稱為位移轉(zhuǎn)換矩陣。4.1桿系結(jié)構(gòu)整體分析4.1.4桿端位移用結(jié)點(diǎn)位移來表示104.1桿系結(jié)構(gòu)整體分析4.1.5整體剛度方程——結(jié)點(diǎn)平衡121221①②③若記引入位移轉(zhuǎn)換關(guān)系，則這就是整體剛度方程，它的物理實(shí)質(zhì)是結(jié)點(diǎn)平衡。[K]稱作結(jié)構(gòu)剛度矩陣（或整體剛度矩陣），[R]稱作綜合等效結(jié)點(diǎn)荷載矩陣，它由兩部分組成。4.1桿系結(jié)構(gòu)整體分析4.1.5整體剛度方程——結(jié)點(diǎn)平衡11單元個數(shù)4.1桿系結(jié)構(gòu)整體分析4.1.6整體剛度矩陣的建立121221①②③若將[A]按單元分成圖示三個子矩陣則由此可見，整體剛度矩陣可由各單元整體剛度矩陣裝配累加得到。為說明如何裝配，先將單元剛度矩陣進(jìn)行分割整體結(jié)點(diǎn)碼則由矩陣乘法可證明，[A]i[k]i[A]iT的結(jié)果是，將剛度矩陣子矩陣按整體結(jié)點(diǎn)碼r、s送整體剛度矩陣相應(yīng)位置。這一裝配規(guī)則稱為“對號入座”。單元個數(shù)4.1桿系結(jié)構(gòu)整體分析4.1.6整體剛度矩陣的建124.1桿系結(jié)構(gòu)整體分析1)任意結(jié)構(gòu)情況上面結(jié)論是通過具體例子（全剛結(jié)點(diǎn)平面剛架）得到的，由虛位移原理或勢能原理進(jìn)行整體分析（見講義），可得任意結(jié)構(gòu)其結(jié)論同此例。2)結(jié)點(diǎn)位移編號如果按結(jié)點(diǎn)順序，對結(jié)點(diǎn)非零位移進(jìn)行依次編號，這一序號稱作結(jié)點(diǎn)位移碼。為便于計算機(jī)處理并減少結(jié)構(gòu)剛度矩陣的階次，將零位移的號碼變?yōu)榱?。對圖示三鉸剛架，當(dāng)僅用一種單元（梁柱自由是單元）時結(jié)點(diǎn)位移編號如圖所示。3)單元定位向量按單元局部結(jié)點(diǎn)碼順序，將結(jié)點(diǎn)位移碼排成的向量，稱作單元的定位向量。①②③④4.1桿系結(jié)構(gòu)整體分析1)任意結(jié)構(gòu)情況上面結(jié)論13①②③④4.1桿系結(jié)構(gòu)整體分析對圖示剛架各單元的定位向量為①(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桿系結(jié)構(gòu)整體分析對圖示剛架各單元的14定位向量①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-1)剛度集裝以3單元為例來說明定位向量單元局部位移碼根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素。定位向量①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-1)剛15根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素，定位向量元素為零時不送。①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-2)荷載集裝以4單元為例來說明定位向量局部位移碼此結(jié)論同樣適用于剛度集裝根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定①②③④①164.1桿系結(jié)構(gòu)整體分析4.1.7整體分析總結(jié)1)對局部坐標(biāo)和整體坐標(biāo)不一致的單元，要對剛度、荷載進(jìn)行坐標(biāo)轉(zhuǎn)換。2)需對“結(jié)構(gòu)”進(jìn)行結(jié)點(diǎn)、位移的局部和整體編號。4)集裝所得整體剛度矩陣是對稱、帶狀稀疏矩陣，當(dāng)支撐條件能限制剛體位移時，矩陣非奇異。3)根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素，定位向量元素為零時不送。據(jù)此可集裝、累加得到整體剛度矩陣。4.1桿系結(jié)構(gòu)整體分析4.1.7整體分析總結(jié)1)對局部175)綜合荷載由兩部分組成，因此首先要將直接作用結(jié)點(diǎn)的荷載按結(jié)點(diǎn)位移碼送入，如果還有單元等效荷載，再按定位向量集裝、累加。4.1桿系結(jié)構(gòu)整體分析8)如果有某位移碼方向彈性支撐，需進(jìn)行將彈簧剛度送入位移碼對應(yīng)的對角線元素位置累加。9)如果有某位移碼方向已知支撐位移，需進(jìn)行將“邊界條件處理”。具體做法以后介紹。7)整體剛度方程實(shí)質(zhì)是全部結(jié)點(diǎn)的平衡條件。6)剛度矩陣帶狀稀疏，其帶寬取決于結(jié)點(diǎn)、位移編碼。最大半帶寬=定位向量中最大元素差+1。5)綜合荷載由兩部分組成，因此首先要將直接作4.1桿系184.1桿系結(jié)構(gòu)整體分析2.5.8邊界條件的處理10)當(dāng)用虛位移或勢能原理作整體分析時（勢能為例），應(yīng)變能為單元應(yīng)變能之和，外力勢能為單元外力勢能之和+結(jié)點(diǎn)外力勢能。全部桿端力的外力勢能彼此抵消。1)乘大數(shù)法設(shè)第i個位移為已知值a，N=108或更大的數(shù)。乘大數(shù)法是將剛度矩陣Kii改為NKii，將Ri改為

Na。請考慮為什麼這樣做能使邊界條件得到滿足？2)置換法（劃零置1）設(shè)第i個位移為已知值a。4.1桿系結(jié)構(gòu)整體分析2.5.8邊界條件的處理10)當(dāng)194.1桿系結(jié)構(gòu)整體分析上述置換工作量大一些，顯然可看出邊界條件得到精確滿足。4.1桿系結(jié)構(gòu)整體分析上述置換工作量大一些，顯然可204.1桿系結(jié)構(gòu)整體分析★3)關(guān)于斜邊界的處理如圖示意的斜支座情況，有多種處理方案。3-1)通過單元的坐標(biāo)轉(zhuǎn)換來處理xy圖示有斜支座單元，r結(jié)點(diǎn)處以傾角-來進(jìn)行坐標(biāo)轉(zhuǎn)換，也即在r結(jié)點(diǎn)處整體坐標(biāo)為圖示xy。r3-2)通過增加一個單元來處理圖示有斜支座單元，r結(jié)點(diǎn)處沿y方向增加一個剛結(jié)的單元，此單元有“無窮大”的抗拉剛度、但沒有抗彎剛度。單元長度可任意。3-3)對整體剛度矩陣進(jìn)行處理（參見匡文起教材）4.1桿系結(jié)構(gòu)整體分析★3)關(guān)于斜邊界的處理21最大半帶寬2.5桿系結(jié)構(gòu)整體分析2.5.10總剛度矩陣的存儲與對應(yīng)解法因?yàn)榭倓偠染仃噷ΨQ、帶狀稀疏，利用這一特點(diǎn)可減少存儲、提高解算效率。零元素零元素非零元素最大半帶寬主對角線元素等半帶存儲零元素非零元素變帶寬一維存儲帶寬是變的到P.55最大2.5桿系結(jié)構(gòu)整體分析2.5.10總剛度矩陣的存儲與222.5桿系結(jié)構(gòu)整體分析目前一般都用變帶寬存儲，下面結(jié)合程序說明存儲和解法。首先介紹一些F90的語法。

定義導(dǎo)出類型導(dǎo)出類型——結(jié)點(diǎn)type::typ_Jointreal::x,y！坐標(biāo)integer::GDOF(3)！整體位移碼endtypetyp_Joint1)有關(guān)F90語法導(dǎo)出類型

新特性

2.5桿系結(jié)構(gòu)整體分析目前一般都用變帶寬存儲，下23用結(jié)點(diǎn)導(dǎo)出類型作為成員導(dǎo)出單元類型:type::typ_Element

integer::JointNo(2)！結(jié)點(diǎn)編號

type(typ_Joint)::Node(2)！結(jié)點(diǎn)

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

real::EA,EIendtypetyp_Elementtype(typ_Element)::Elem(5)

！定義5個單元類型…!對單元i的端點(diǎn)j的x,y,GDOF(1:3)的賦值Elem(i)%Node(j)=Joint(Elem(i)%JointNo(j))…由導(dǎo)出類型定義新類型由導(dǎo)出類型定義變量用結(jié)點(diǎn)導(dǎo)出類型作為成員導(dǎo)出單元類型: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

!數(shù)組與標(biāo)量運(yùn)算:A(1:5)+(/2,2,2,2,2/)A=C+A

!數(shù)組與數(shù)組運(yùn)算（同形）C=sqrt(A)!數(shù)組的函數(shù)運(yùn)算:C(i)=sqrt(A(i),i=1,5數(shù)組內(nèi)部函數(shù)：dot_product(vector_a,vector_b)

!點(diǎn)積如：dot_product((/1,2,3/),(/2,3,4/))的值為20（待續(xù)）有關(guān)F90語法數(shù)組運(yùn)算與賦值：real::A(5),B(5,10),C(5)有關(guān)F925matmul(matrix_a,matrix_b)

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

矩陣轉(zhuǎn)置如：glbEDisp=matmul(transpose(T),locEDisp)size(array,dim)

求數(shù)組第dim維的長度dim為可選變元:size(a,dim=2)若array為一維時，可不用dim。sum(array,dim,mask)

數(shù)組元素求和dim,mask為可選變元；mask=條件表達(dá)式sum(a(1:10))對a的1到10元素求和sum(a(1:10),mask=a>0)

對a(1:10)中大于0的元素求和（續(xù)）matmul(matrix_a,matrix_b)矩26有關(guān)F90語法where結(jié)構(gòu)

新特性例where(C>0)C=0A=B*Dendwherewhere(C>0)A=Bendwhere定義where(數(shù)組關(guān)系表達(dá)式)

數(shù)組賦值語句

…elsewhere

數(shù)組賦值語句

...endwhere規(guī)則:

1)同形數(shù)組;2)不許嵌套;3)最多二分叉有關(guān)F90語法例定義規(guī)則:1)同形數(shù)組;2)不許嵌套27有關(guān)F90語法cycle和exit語句

新特性

用在do循環(huán)中

cycle——作下一個循環(huán)步

exit——跳出循環(huán)，執(zhí)行enddo后一條語句等效例

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

...if(cond1)cycle ...

if(cond2)exit...enddo...有關(guān)F90語法等效例用法28有關(guān)F90語法數(shù)組構(gòu)造函數(shù)spread語法spread(數(shù)組名,dim,ncopies)

將數(shù)組沿dim維方向復(fù)制ncopies形成新數(shù)組

dim,ncopies

—整型、位置變元、關(guān)鍵字變元（若按位置引用，可略關(guān)鍵字）例:（僅限一維數(shù)組）1）

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

[1,1,1]或[1,1,1]T有關(guān)F90語法語法例:（僅限一維數(shù)組）[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有關(guān)F90語法指針pointerpointer是變量的屬性，可以指向相同類型的變量;被指向的目標(biāo)必須具有target屬性或pointer屬性可以將指針變量理解為別名、稱號real,target::a,b,EDisp(6)

!可被指針?biāo)?/p>

real,pointer::p1,p2!稱號：班長、課代表！p1,p2是指針，可以指向?qū)嵭蛿?shù)據(jù)real,pointer::G(:)!先進(jìn)集體！G是指針，可以指向一維實(shí)型數(shù)組有關(guān)F90語法pointer是變量的屬性，可以指向相同類31指針是一種“稱號”，上述聲明語句建立了“稱號”，但并未“授予”某個變量這個稱號，因此是指向“空”，并未占用內(nèi)存a=3.0p1=>a!p1指向a!稱號p1授予a，a的數(shù)據(jù)有兩個名：固定名a和流動名p1;既可用p1也可用a（p1—班長，a—張三）a=4.0!a的值變?yōu)?.0,p1也變?yōu)?.0p1=>b!班長換人了G=>EDisp!先進(jìn)集體有了得主！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指向的數(shù)據(jù)b)=(p指向的數(shù)據(jù)a)！即：所有=3.14又如：33指針可以指向有名的數(shù)據(jù)區(qū)，也可以指向無名的數(shù)據(jù)區(qū)real,pointer::b(:,:)integer::nread(*,*)nallocate(b(n,n))

!指針指向了一個剛開辟的數(shù)組!以下可以當(dāng)作常規(guī)數(shù)組用b(1,1)=1.1b(1,2)=1.2...deallocate(b)有關(guān)F90語法用指針建立動態(tài)數(shù)組指針可以指向有名的數(shù)據(jù)區(qū)，也可以指向無名的數(shù)據(jù)區(qū)有關(guān)F9034指針與allocatable數(shù)組的區(qū)別具備allocatable數(shù)組的所有功能還可以用在導(dǎo)出類型中，例如整體剛度矩陣的變帶寬存儲：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數(shù)組的區(qū)別35(1)NElem,NJoint,NGlbDOF,NJLoad,NELoad單元數(shù)結(jié)點(diǎn)數(shù)總自由度數(shù)結(jié)點(diǎn)荷載數(shù)單元荷載數(shù)[Joint-結(jié)點(diǎn)]…NJoint行(2)Joint%X,Joint%Y,GDOF(1:3)X坐標(biāo)Y坐標(biāo)結(jié)點(diǎn)位移碼[Elem-單元]…NElem行(3)JointNo1,JointNo2,EA,EI起點(diǎn)號終點(diǎn)號剛度[JLoad-結(jié)點(diǎn)荷載]…NJLoad行(4)JointNo,LodDOF,LodVal作用點(diǎn)號局部位移碼荷載值[ELoad-單元荷載]…NELoad行(5)ElemNo,Indx,Pos,LodVal單元號類型號位置荷載值Indx

類型

pos1均布長度2集中位置3...2)某程序輸入數(shù)據(jù)說明(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

!結(jié)點(diǎn)2，自由度1，值為10E32,1,4,-4.0E3

！單元2，均布，長4米，值為-4E32-1)數(shù)據(jù)文件例子：(2)(1)(3)24135i=6i=4i=310kN4kN/m4m4mEA=109N(1)(2)坐標(biāo)位移碼(3)(4)(5)結(jié)點(diǎn)號EA,EI3,5,7,1,12-1)數(shù)據(jù)文件例子：(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)程序讀入計算所需數(shù)據(jù)：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為自動數(shù)組，子程序結(jié)束后自動釋放。！這樣做可使接口簡單，減少了數(shù)組的控制變量。

integer(ikind)::ie,j,NGlbDOF,NElemNGlbDOF=size(Kcol,dim=1)！使接口簡單2-3)求始行碼和分配帶寬子程序!===========39NElem=size(Elem,dim=1)Row1=NGlbDOF

!先設(shè)所有始行碼同終行碼!確定各列始行碼向量

doie=1,Nelem!對單元循環(huán)!確定定位向量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

!對總自由度數(shù)循環(huán)allocate(Kcol(j)%row(Row1(j):j)

!給Kcol(j)%row從Row1(j)到j(luò)個空間Kcol(j)%row=zero

!總剛元素清零

enddoreturnendsubroutineSetMatBand!為各列的帶寬分配空間413)整體剛度矩陣的集成

doie=1,NElem…!計算單剛

EK(6,6),ELocVec(6)doj=1,6!對單元逐列集成JGDOF=ELocVec(j)!取出位移碼if(JGDOF=0)cycle!作下一循環(huán)步where(ELocVec>0.And.ELocVec<=JGDOF)

!位移碼非零同時小于第j個局部碼對應(yīng)的位移碼Kcol(JGDOF)%row(ELocVec)=&Kcol(JGDOF)%row(ELocVec)+EK(:,j)!集成endwhereenddoenddo局部碼位移碼3)整體剛度矩陣的集成doie=1,NElem局部424)變帶寬矩陣的分解求解4-1)LDLT分解法求解[A]{x}

=若對稱正定，則可分解為=單位上三角陣對角陣原方程變?yōu)榍蠼獠襟E：1.分解：2.解y：3.解x：LDLT分解法Gauss消去前消去處理右端項(xiàng)回代(不同的b只做一次分解)4)變帶寬矩陣的分解求解若對稱正定，則可分解為=單位上三角43存放：主對角—

上三角—

存放：主對角—上三角—44時不必求和(上三角:i<j)不動存在處存在處時不必求和(上三角:i<j)不動存在處存在45(第j列系數(shù))上三角:i<j

4-2)分解一般情況：逐列分解對角線:i=

j>1(第j列系數(shù))上三角:i<j4-2)分解一般情46(第j列系數(shù))第

i列中第1個非零元素的行碼為：4-3)變列寬存貯的分解修正：

第

j

列中第1個非零元素的行碼為：

則分解順序(第j列系數(shù))第i列中第1個非零元素的行碼為：4-3474-4)F90實(shí)現(xiàn)!三角分解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列系數(shù)分解完畢

s=sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)diag(j)=diag(j)-s

!第j列的主對角元素enddo4-4)F90實(shí)現(xiàn)!三角分解484-5)一般情況解y公式可不動可存在處4-5)一般情況解y公式可不動可存在處494-6)變列寬存儲解y的修正上三角的第

i列從第行元素開始!回代步驟1:自上而下doi=2,ncolrow1=lbound(Kcol(i)%row,dim=1)GP(i)=GP(i)-sum(Kcol(i)%row(row1:i-1)&*GP(row1:i-1))enddo

求出后，不再用，可將存在處4-6)變列寬存儲解y的修正上三角的第i列從第行504-7)一般解x的公式自下而上：4-7)一般解x的公式自下而上：51這樣做的缺欠：自下而上逐行計算，行中遇到0元素需跳過，不方便！4-8)改為自右向左逐列計算！記解出對y向量預(yù)處理再解出第n列乘上后移到右邊去，修正y向量這樣做的缺欠：自下而上逐行計算，行中遇到0元素需跳過，不方便52!回代步驟2:自右向左GP(:)=GP(:)/diag(:)doj=ncol,2,-1row1=lbound(Kcol(j)%row,dim=1)GP(row1:j-1)=GP(row1:j-1)-GP(j)*&Kcol(j)%row(row1:j-1)enddo!現(xiàn)在GP就是解!回代步驟2:自右向左534-9)小結(jié)：

無須一維存貯，無須行列碼轉(zhuǎn)換定位公式與程序直接對應(yīng)翻譯，直截了當(dāng)求和采用內(nèi)部函數(shù)進(jìn)行數(shù)組運(yùn)算動態(tài)內(nèi)存，用多大、開多大數(shù)據(jù)封裝性好，接口簡單:

subroutineSetMatBand(Kcol,Elem)subroutineVarBandSolv(Disp,Kcol,GP)通用性強(qiáng)4-9)小結(jié)：542.6桿系結(jié)構(gòu)內(nèi)力計算2.6.1桿端內(nèi)力計算公式解方程的結(jié)果可以得到結(jié)點(diǎn)位移，有了位移就可以進(jìn)一步求各單元的內(nèi)力。

解算步驟：從[]根據(jù)定位向量取出單元桿端位移。由單元傾角確定是算還是。減去等效結(jié)點(diǎn)荷載或加上固端力矩陣。2.6桿系結(jié)構(gòu)內(nèi)力計算2.6.1桿端內(nèi)力計算公式552.6桿系結(jié)構(gòu)內(nèi)力計算2.6.2桿中任意截面內(nèi)力計算公式

注意事項(xiàng)：由圖示隔離體圖，列桿段的平衡方程即可得到任意截面的內(nèi)力計算公式。請自行寫出。按上述隔離體圖所求得的內(nèi)力是精確的。求得單元結(jié)點(diǎn)位移后，代入單元位移場、求應(yīng)變、求內(nèi)力，這樣所得的內(nèi)力一般不滿足平衡條件。只是近似解。2.6桿系結(jié)構(gòu)內(nèi)力計算2.6.2桿中任意截面內(nèi)力計算公式56SubroutineElem_Stiff(······)

說明Stiff=0.0!單元剛度清零SelectCase(Type)Case(1)

平面桿系結(jié)構(gòu)單元Case(2)

空間桿系結(jié)構(gòu)單元CaseDefault

出錯信息EndSelectEndSubroutineElem_Stiff3.4桿系結(jié)構(gòu)單元分析子程序3.4.1單元剛度總體設(shè)計3.4.2說明部分設(shè)計Integer,Intent(in)::···

入口整型參數(shù)Real(8),Intent(in)::···

入口實(shí)型參數(shù)Real(8),Intent(out)::···

出口實(shí)型參數(shù)Real(8)::Work1,···Integer::i,j,k,···實(shí)型和整型工作變量SubroutineElem_Stiff(······)3573.4桿系結(jié)構(gòu)單元分析子程序3.4.3平面桿系結(jié)構(gòu)設(shè)計SelectCase(Plane)Case(1)

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

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

出錯信息EndSelect3.4.4空間桿系結(jié)構(gòu)設(shè)計SelectCase(Space)Case(1)

空間桁架元素賦值Case(2)

空間梁柱元素賦值Case(3)

交叉梁元素賦值CaseDefault

出錯信息EndSelect3.4桿系結(jié)構(gòu)單元分析子程序3.4.3平面桿系結(jié)構(gòu)設(shè)計583.4桿系結(jié)構(gòu)單元分析子程序3.4.5有關(guān)單元等效結(jié)點(diǎn)荷載設(shè)計和進(jìn)一步的考慮1)單元等效結(jié)點(diǎn)荷載設(shè)計同仿單元剛度。2)從各類單元剛度元素的計算，可看到要用到長度、單元彈性特性、單元截面特性等數(shù)據(jù)。因此，要確定存放它們的數(shù)據(jù)結(jié)構(gòu)。要將它們作為出口。3)為計算單元等效結(jié)點(diǎn)荷載元素，首先要建立各種荷載情況等效荷載表達(dá)式，它們可由積分或載常數(shù)表得到。然后要解決荷載信息的存放結(jié)構(gòu)，也要將它們作為出口量。4)單元剛度矩陣、等效結(jié)點(diǎn)荷載矩陣都應(yīng)先清零。3.4桿系結(jié)構(gòu)單元分析子程序3.4.5有關(guān)單元等效結(jié)點(diǎn)荷594.1桿系結(jié)構(gòu)整體分析

首先就全剛結(jié)點(diǎn)平面剛架進(jìn)行討論，然后推廣。4.1.1總的思路在單元特性搞清后，將單元拼裝回去。在結(jié)點(diǎn)處位移自動協(xié)調(diào)基礎(chǔ)上，如果全部結(jié)點(diǎn)平衡，則求得的結(jié)點(diǎn)位移將是實(shí)際結(jié)構(gòu)的解。因此，整體分析就是設(shè)法建立結(jié)點(diǎn)平衡方程。4.1.2坐標(biāo)轉(zhuǎn)換組成結(jié)構(gòu)的桿件可以各個方向，單元分析對局部坐標(biāo)，因此，必須將物理量轉(zhuǎn)為統(tǒng)一坐標(biāo)——整體坐標(biāo)。1)力的轉(zhuǎn)換關(guān)系4.1桿系結(jié)構(gòu)整體分析首先就全剛結(jié)點(diǎn)平面剛架進(jìn)行604.1桿系結(jié)構(gòu)整體分析2)位移轉(zhuǎn)換關(guān)系3)轉(zhuǎn)換矩陣轉(zhuǎn)換矩陣是正交矩陣。4.1桿系結(jié)構(gòu)整體分析2)位移轉(zhuǎn)換關(guān)系3)轉(zhuǎn)換矩陣614.1桿系結(jié)構(gòu)整體分析4)桿端力轉(zhuǎn)換5)桿端位移轉(zhuǎn)換6)剛度方程的轉(zhuǎn)換如果記稱為整體單元剛度矩陣則這就是整體坐標(biāo)下的單元剛度方程。本節(jié)以后的討論認(rèn)為都是對整體坐標(biāo)的4.1桿系結(jié)構(gòu)整體分析4)桿端力轉(zhuǎn)換5)桿端位移轉(zhuǎn)換6624.1桿系結(jié)構(gòu)整體分析4.1.3結(jié)點(diǎn)平衡方程的建立1)一簡單例子（如圖）圖中有兩套編號，紅的是單元桿端編號，黑的是結(jié)構(gòu)整體編號。1-1)結(jié)點(diǎn)示意121221①②③圖中藍(lán)色的表示結(jié)點(diǎn)荷載（已知），紅色的表示桿端力（未知的），、分別1、2單元桿端力子矩陣。對1、4結(jié)點(diǎn)“荷載”含有未知反力。21-2)結(jié)點(diǎn)平衡由示意圖可見，結(jié)構(gòu)結(jié)點(diǎn)的平衡方程為4.1桿系結(jié)構(gòu)整體分析4.1.3結(jié)點(diǎn)平衡方程的建立1)634.1桿系結(jié)構(gòu)整體分析從例圖可見，其全部結(jié)點(diǎn)平衡方程為121221①②③若記24.1桿系結(jié)構(gòu)整體分析從例圖可見，其全部結(jié)點(diǎn)12644.1桿系結(jié)構(gòu)整體分析式中[I]、[0]分別為單位和零矩陣。若引入矩陣記號，則結(jié)點(diǎn)平衡方程可改寫作這一結(jié)論雖然是由一個例子得到的，但是顯然對一切結(jié)構(gòu)都是成立的。問題在于不同結(jié)構(gòu)，[A]矩陣是不同的。4.1桿系結(jié)構(gòu)整體分析式中[I]、[0]分別為單位和零矩陣654.1桿系結(jié)構(gòu)整體分析4.1.4桿端位移用結(jié)點(diǎn)位移來表示121221①②③仍以簡單例子來說明若記由結(jié)點(diǎn)、桿端位移的協(xié)調(diào)條件，可得[]、[]的對應(yīng)關(guān)系為式中[A]T是前面力關(guān)系[A]的轉(zhuǎn)置，因此[A]T稱為位移轉(zhuǎn)換矩陣。4.1桿系結(jié)構(gòu)整體分析4.1.4桿端位移用結(jié)點(diǎn)位移來表示664.1桿系結(jié)構(gòu)整體分析4.1.5整體剛度方程——結(jié)點(diǎn)平衡121221①②③若記引入位移轉(zhuǎn)換關(guān)系，則這就是整體剛度方程，它的物理實(shí)質(zhì)是結(jié)點(diǎn)平衡。[K]稱作結(jié)構(gòu)剛度矩陣（或整體剛度矩陣），[R]稱作綜合等效結(jié)點(diǎn)荷載矩陣，它由兩部分組成。4.1桿系結(jié)構(gòu)整體分析4.1.5整體剛度方程——結(jié)點(diǎn)平衡67單元個數(shù)4.1桿系結(jié)構(gòu)整體分析4.1.6整體剛度矩陣的建立121221①②③若將[A]按單元分成圖示三個子矩陣則由此可見，整體剛度矩陣可由各單元整體剛度矩陣裝配累加得到。為說明如何裝配，先將單元剛度矩陣進(jìn)行分割整體結(jié)點(diǎn)碼則由矩陣乘法可證明，[A]i[k]i[A]iT的結(jié)果是，將剛度矩陣子矩陣按整體結(jié)點(diǎn)碼r、s送整體剛度矩陣相應(yīng)位置。這一裝配規(guī)則稱為“對號入座”。單元個數(shù)4.1桿系結(jié)構(gòu)整體分析4.1.6整體剛度矩陣的建684.1桿系結(jié)構(gòu)整體分析1)任意結(jié)構(gòu)情況上面結(jié)論是通過具體例子（全剛結(jié)點(diǎn)平面剛架）得到的，由虛位移原理或勢能原理進(jìn)行整體分析（見講義），可得任意結(jié)構(gòu)其結(jié)論同此例。2)結(jié)點(diǎn)位移編號如果按結(jié)點(diǎn)順序，對結(jié)點(diǎn)非零位移進(jìn)行依次編號，這一序號稱作結(jié)點(diǎn)位移碼。為便于計算機(jī)處理并減少結(jié)構(gòu)剛度矩陣的階次，將零位移的號碼變?yōu)榱?。對圖示三鉸剛架，當(dāng)僅用一種單元（梁柱自由是單元）時結(jié)點(diǎn)位移編號如圖所示。3)單元定位向量按單元局部結(jié)點(diǎn)碼順序，將結(jié)點(diǎn)位移碼排成的向量，稱作單元的定位向量。①②③④4.1桿系結(jié)構(gòu)整體分析1)任意結(jié)構(gòu)情況上面結(jié)論69①②③④4.1桿系結(jié)構(gòu)整體分析對圖示剛架各單元的定位向量為①(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桿系結(jié)構(gòu)整體分析對圖示剛架各單元的70定位向量①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-1)剛度集裝以3單元為例來說明定位向量單元局部位移碼根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素。定位向量①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-1)剛71根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素，定位向量元素為零時不送。①②③④①②③④4.1桿系結(jié)構(gòu)整體分析4-2)荷載集裝以4單元為例來說明定位向量局部位移碼此結(jié)論同樣適用于剛度集裝根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定①②③④①724.1桿系結(jié)構(gòu)整體分析4.1.7整體分析總結(jié)1)對局部坐標(biāo)和整體坐標(biāo)不一致的單元，要對剛度、荷載進(jìn)行坐標(biāo)轉(zhuǎn)換。2)需對“結(jié)構(gòu)”進(jìn)行結(jié)點(diǎn)、位移的局部和整體編號。4)集裝所得整體剛度矩陣是對稱、帶狀稀疏矩陣，當(dāng)支撐條件能限制剛體位移時，矩陣非奇異。3)根據(jù)單元局部位移碼和定位向量的對應(yīng)關(guān)系用定位向量位移碼送元素，定位向量元素為零時不送。據(jù)此可集裝、累加得到整體剛度矩陣。4.1桿系結(jié)構(gòu)整體分析4.1.7整體分析總結(jié)1)對局部735)綜合荷載由兩部分組成，因此首先要將直接作用結(jié)點(diǎn)的荷載按結(jié)點(diǎn)位移碼送入，如果還有單元等效荷載，再按定位向量集裝、累加。4.1桿系結(jié)構(gòu)整體分析8)如果有某位移碼方向彈性支撐，需進(jìn)行將彈簧剛度送入位移碼對應(yīng)的對角線元素位置累加。9)如果有某位移碼方向已知支撐位移，需進(jìn)行將“邊界條件處理”。具體做法以后介紹。7)整體剛度方程實(shí)質(zhì)是全部結(jié)點(diǎn)的平衡條件。6)剛度矩陣帶狀稀疏，其帶寬取決于結(jié)點(diǎn)、位移編碼。最大半帶寬=定位向量中最大元素差+1。5)綜合荷載由兩部分組成，因此首先要將直接作4.1桿系744.1桿系結(jié)構(gòu)整體分析2.5.8邊界條件的處理10)當(dāng)用虛位移或勢能原理作整體分析時（勢能為例），應(yīng)變能為單元應(yīng)變能之和，外力勢能為單元外力勢能之和+結(jié)點(diǎn)外力勢能。全部桿端力的外力勢能彼此抵消。1)乘大數(shù)法設(shè)第i個位移為已知值a，N=108或更大的數(shù)。乘大數(shù)法是將剛度矩陣Kii改為NKii，將Ri改為

Na。請考慮為什麼這樣做能使邊界條件得到滿足？2)置換法（劃零置1）設(shè)第i個位移為已知值a。4.1桿系結(jié)構(gòu)整體分析2.5.8邊界條件的處理10)當(dāng)754.1桿系結(jié)構(gòu)整體分析上述置換工作量大一些，顯然可看出邊界條件得到精確滿足。4.1桿系結(jié)構(gòu)整體分析上述置換工作量大一些，顯然可764.1桿系結(jié)構(gòu)整體分析★3)關(guān)于斜邊界的處理如圖示意的斜支座情況，有多種處理方案。3-1)通過單元的坐標(biāo)轉(zhuǎn)換來處理xy圖示有斜支座單元，r結(jié)點(diǎn)處以傾角-來進(jìn)行坐標(biāo)轉(zhuǎn)換，也即在r結(jié)點(diǎn)處整體坐標(biāo)為圖示xy。r3-2)通過增加一個單元來處理圖示有斜支座單元，r結(jié)點(diǎn)處沿y方向增加一個剛結(jié)的單元，此單元有“無窮大”的抗拉剛度、但沒有抗彎剛度。單元長度可任意。3-3)對整體剛度矩陣進(jìn)行處理（參見匡文起教材）4.1桿系結(jié)構(gòu)整體分析★3)關(guān)于斜邊界的處理77最大半帶寬2.5桿系結(jié)構(gòu)整體分析2.5.10總剛度矩陣的存儲與對應(yīng)解法因?yàn)榭倓偠染仃噷ΨQ、帶狀稀疏，利用這一特點(diǎn)可減少存儲、提高解算效率。零元素零元素非零元素最大半帶寬主對角線元素等半帶存儲零元素非零元素變帶寬一維存儲帶寬是變的到P.55最大2.5桿系結(jié)構(gòu)整體分析2.5.10總剛度矩陣的存儲與782.5桿系結(jié)構(gòu)整體分析目前一般都用變帶寬存儲，下面結(jié)合程序說明存儲和解法。首先介紹一些F90的語法。

定義導(dǎo)出類型導(dǎo)出類型——結(jié)點(diǎn)type::typ_Jointreal::x,y！坐標(biāo)integer::GDOF(3)！整體位移碼endtypetyp_Joint1)有關(guān)F90語法導(dǎo)出類型

新特性

2.5桿系結(jié)構(gòu)整體分析目前一般都用變帶寬存儲，下79用結(jié)點(diǎn)導(dǎo)出類型作為成員導(dǎo)出單元類型:type::typ_Element

integer::JointNo(2)！結(jié)點(diǎn)編號

type(typ_Joint)::Node(2)！結(jié)點(diǎn)

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

real::EA,EIendtypetyp_Elementtype(typ_Element)::Elem(5)

！定義5個單元類型…!對單元i的端點(diǎn)j的x,y,GDOF(1:3)的賦值Elem(i)%Node(j)=Joint(Elem(i)%JointNo(j))…由導(dǎo)出類型定義新類型由導(dǎo)出類型定義變量用結(jié)點(diǎn)導(dǎo)出類型作為成員導(dǎo)出單元類型:type(typ_El80real::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

!數(shù)組與標(biāo)量運(yùn)算:A(1:5)+(/2,2,2,2,2/)A=C+A

!數(shù)組與數(shù)組運(yùn)算（同形）C=sqrt(A)!數(shù)組的函數(shù)運(yùn)算:C(i)=sqrt(A(i),i=1,5數(shù)組內(nèi)部函數(shù)：dot_product(vector_a,vector_b)

!點(diǎn)積如：dot_product((/1,2,3/),(/2,3,4/))的值為20（待續(xù)）有關(guān)F90語法數(shù)組運(yùn)算與賦值：real::A(5),B(5,10),C(5)有關(guān)F981matmul(matrix_a,matrix_b)

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

矩陣轉(zhuǎn)置如：glbEDisp=matmul(transpose(T),locEDisp)size(array,dim)

求數(shù)組第dim維的長度dim為可選變元:size(a,dim=2)若array為一維時，可不用dim。sum(array,dim,mask)

數(shù)組元素求和dim,mask為可選變元；mask=條件表達(dá)式sum(a(1:10))對a的1到10元素求和sum(a(1:10),mask=a>0)

對a(1:10)中大于0的元素求和（續(xù)）matmul(matrix_a,matrix_b)矩82有關(guān)F90語法where結(jié)構(gòu)

新特性例where(C>0)C=0A=B*Dendwherewhere(C>0)A=Bendwhere定義where(數(shù)組關(guān)系表達(dá)式)

數(shù)組賦值語句

…elsewhere

數(shù)組賦值語句

...endwhere規(guī)則:

1)同形數(shù)組;2)不許嵌套;3)最多二分叉有關(guān)F90語法例定義規(guī)則:1)同形數(shù)組;2)不許嵌套83有關(guān)F90語法cycle和exit語句

新特性

用在do循環(huán)中

cycle——作下一個循環(huán)步

exit——跳出循環(huán)，執(zhí)行enddo后一條語句等效例

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

...if(cond1)cycle ...

if(cond2)exit...enddo...有關(guān)F90語法等效例用法84有關(guān)F90語法數(shù)組構(gòu)造函數(shù)spread語法spread(數(shù)組名,dim,ncopies)

將數(shù)組沿dim維方向復(fù)制ncopies形成新數(shù)組

dim,ncopies

—整型、位置變元、關(guān)鍵字變元（若按位置引用，可略關(guān)鍵字）例:（僅限一維數(shù)組）1）

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

[1,1,1]或[1,1,1]T有關(guān)F90語法語法例:（僅限一維數(shù)組）[1,1,1852）

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,086有關(guān)F90語法指針pointerpointer是變量的屬性，可以指向相同類型的變量;被指向的目標(biāo)必須具有target屬性或pointer屬性可以將指針變量理解為別名、稱號real,target::a,b,EDisp(6)

!可被指針?biāo)?/p>

real,pointer::p1,p2!稱號：班長、課代表！p1,p2是指針，可以指向?qū)嵭蛿?shù)據(jù)real,pointer::G(:)!先進(jìn)集體！G是指針，可以指向一維實(shí)型數(shù)組有關(guān)F90語法pointer是變量的屬性，可以指向相同類87指針是一種“稱號”，上述聲明語句建立了“稱號”，但并未“授予”某個變量這個稱號，因此是指向“空”，并未占用內(nèi)存a=3.0p1=>a!p1指向a!稱號p1授予a，a的數(shù)據(jù)有兩個名：固定名a和流動名p1;既可用p1也可用a（p1—班長，a—張三）a=4.0!a的值變?yōu)?.0,p1也變?yōu)?.0p1=>b!班長換人了G=>EDisp!先進(jìn)集體有了得主！EDisp(:)的長度就是G(:)的長度，用G和用EDisp同樣效果指針是一種“稱號”，上述聲明語句建立了“稱號”，但并未“授予88又如：real,target::a,breal,pointer::p,qa=3.14b=2.0p=>a!p=a=3.14q=>b!q=b=2.0q=p!(q指向的數(shù)據(jù)b)=(p指向的數(shù)據(jù)a)！即：所有=3.14又如：89指針可以指向有名的數(shù)據(jù)區(qū)，也可以指向無名的數(shù)據(jù)區(qū)real,pointer::b(:,:)integer::nread(*,*)nallocate(b(n,n))

!指針指向了一個剛開辟的數(shù)組!以下可以當(dāng)作常規(guī)數(shù)組用b(1,1)=1.1b(1,2)=1.2...deallocate(b)有關(guān)F90語法用指針建立動態(tài)數(shù)組指針可以指向有名的數(shù)據(jù)區(qū)，也可以指向無名的數(shù)據(jù)區(qū)有關(guān)F9090指針與allocatable數(shù)組的區(qū)別具備allocatable數(shù)組的所有功能還可以用在導(dǎo)出類型中，例如整體剛度矩陣的變帶寬存儲：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數(shù)組的區(qū)別91(1)NElem,NJoint,NGlbDOF,NJLoad,NELoad單元數(shù)結(jié)點(diǎn)數(shù)總自由度數(shù)結(jié)點(diǎn)荷載數(shù)單元荷載數(shù)[Joint-結(jié)點(diǎn)]…NJoint行(2)Joint%X,Joint%Y,GDOF(1:3)X坐標(biāo)Y坐標(biāo)結(jié)點(diǎn)位移碼[Elem-單元]…NElem行(3)JointNo1,JointNo2,EA,EI起點(diǎn)號終點(diǎn)號剛度[JLoad-結(jié)點(diǎn)荷載]…NJLoad行(4)JointNo,LodDOF,LodVal作用點(diǎn)號局部位移碼荷載值[ELoad-單元荷載]…NELoad行(5)ElemNo,Indx,Pos,LodVal單元號類型號位置荷載值Indx

類型

pos1均布長度2集中位置3...2)某程序輸入數(shù)據(jù)說明(1)NElem,NJoint,NGlbDOF,923,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

!結(jié)點(diǎn)2，自由度1，值為10E32,1,4,-4.0E3

！單元2，均布，長4米，值為-4E32-1)數(shù)據(jù)文件例子：(2)(1)(3)24135i=6i=4i=310kN4kN/m4m4mEA=109N(1)(2)坐標(biāo)位移碼(3)(4)(5)結(jié)點(diǎn)號EA,EI3,5,7,1,12-1)數(shù)據(jù)文件例子：(2)(1)(3)93read(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)程序讀入計算所需數(shù)據(jù)：read(5,*)NElem,NJoint,NGlbDOF942-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為自動數(shù)組，子程序結(jié)束后自動釋放。！這樣做可使接口簡單，減少了數(shù)組的控制變量。

integer(ikind)::ie,j,NGlbDOF,NElemNGlbDOF=size(Kcol,dim=1)！使接口簡單2-3)求始行碼和分配帶寬子程序!===========95NElem=size(Elem,dim=1)Row1=NGlbDOF

!先設(shè)所有始行碼同終行碼!確定各列始行碼向量

doie=1,Nelem!對單元循環(huán)!確定定位向量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)96!為各列的帶寬分配空間

doj=1,NGlbDOF

!對總自由度數(shù)循環(huán)allocate(Kcol(j)%row(Row1(j):j)

!給Kcol(j)%row從Row1(j)到j(luò)個空間Kcol(j)%row=zero