1、一個牛頓拉弗遜法求解非線性問題迭代過程的簡單程序例子! Example source code for Newton-Raphson Method in Non-linear Problems ! 可以進行荷載步長自動調整! Load step can be automatically adjusted! 讀者可參考該模塊加入相應的子程序編制自己的非線性計算程序module Main_Iterationuse Lxz_Tools ! 通用子程序模塊 use TypeDef ! 數據結構定義模塊 use Data_Input ! 數據輸入模塊 use Data_Output ! 數據輸出模塊 u

2、se Elem_Prop ! 單元屬性模塊 use MatSolve ! 矩陣求解模塊 implicit none contains subroutine RC2D_Main() type(typ_GValue) : GValue ! 總體控制變量數據結構數組節點數據結構數組 混凝土單元數據結構數組 鋼筋單元數據結構數組荷載數據結構數組 初始荷載 材料數據結構數組 支座數據結構數組type(typ_Node),pointer : Node(:) ! type(typ_Elem),pointer : Elem(:) ! type(typ_Rebar),pointer : Rebar(:) ! t

3、ype(typ_Load),pointer : Load(:) ! type(typ_Load),pointer : Initial_Load(:) ! type(typ_Material),pointer : Material(:) ! type(typ_Support),pointer : Support(:) !call DataInput(GValue,Node,Elem,Rebar,Material,Load,Support,& Initial_Load) ! 讀入數據! 迭代輸出數call Iteration(GValue,Node,Elem,Rebar,Load,Supp

4、ort,Initial_Load) 核心程序call DataOutput(GValue,Node,Elem,Rebar,Material,Load,Support) ! 據write(*,*) "計算成功結束 "read(*,*) deallocate(node) deallocate(Elem) deallocate(Rebar) deallocate(Load) deallocate(Support)return end subroutine subroutine Iteration(GValue0,Node0,Elem0,Rebar0,Load0,Support0,

5、Initial_Load0) type(typ_GValue) : GValue0 type(typ_Node) : Node0(:) type(typ_Elem) : Elem0(:) type(typ_Rebar) : Rebar0(:) type(typ_Load): Load0(:) type(typ_Load): Initial_Load0(:) type(typ_Support) : Support0(:) type(typ_GValue) : GValue type(typ_Node) : Node(GValue0%NNode) type(typ_Elem) : Elem(GVa

6、lue0%NElem) type(typ_Rebar) : Rebar(GValue0%NRebar) type(typ_Load): Load(GValue0%NLoad) type(typ_Load): Initial_Load(GValue0%NInitial_Load) type(typ_Support) : Support(GValue0%NSupport) Integer(ikind) : I,J,II,IIIGValue=GValue0 GValue%FinishedStep=0II=1 ! 開始計算加載步數LOOP_A: dowrite(*,*) "= write(*

7、,*) 'Step',II," 總加載步數 ",GValue%NStepwrite(*,*) ' 完成 ',(GValue%FinishedStep)/real(GValue%NStep)*100.0,'%'! 對數據進行備份，在荷載步調整操作時恢復上一步的參數 do I=1,GValue%NNode; Node(I)=Node0(I); end do do I=1,GValue%NElem; Elem(I)=Elem0(I); end do do I=1,GValue%NRebar; Rebar(I)=Rebar0(I);

8、end do do I=1,Gvalue%NLoad; Load(I)=Load0(I); end do do I=1,Gvalue%NInitial_Load; Initial_Load(I)=Initial_Load0(I); end do do I=1,GValue%NSupport; Support(I)=Support0(I); end doGValue%GError=1.0 ! 初始誤差設定為 1.0 != GValue%ForError=0.0d0 ! 前一次迭代誤差清零III=1LOOP_B: do write(*,*) 'Iteration Number:',

9、 III call Get_Elem_B(GValue,Elem,Node) ! 計算形函數 call Get_Elem_D(GValue,Elem) ! 計算混凝土單元材料本構矩陣 call Get_Elem_K(GValue,Elem) ! 計算混凝土單元剛度矩陣 call Get_Rebar_K(GValue,Rebar,Node) ! 計算鋼筋單元剛度矩陣 call Get_Elem_F(GValue,Elem) ! 計算混凝土單元荷載向量 call Get_Rebar_F(GValue,Rebar) ! 計算鋼筋單元荷載向量 ! 計算節點不平衡力call Cal_Node_DForc

10、e(GValue,Node,Elem,Rebar,Load,Initial_Load)! 計算荷載誤差call Get_Error_F(GValue, Node, Load, Initial_Load,support)if(III.ne.1) then !開始迭代! call Get_Error_F(GValue, Node) ! write(*,*) " 當前迭代誤差 :",GValue%GError if(GValue%GError<=GValue%ErrorT) then ! 誤差小于誤差容限，進行下一步計 算!do I=1,GValue%NElem! 更新相應

11、的參數 do J=1,4Elem(I)%Concrete(J)%EPS=Elem(I)%Concrete(J)%STRAIN Elem(I)%Concrete(J)%SIG=Elem(I)%Concrete(J)%Stress Elem(I)%SRebar(J)%EPS=Elem(I)%SRebar(J)%STRAIN Elem(I)%SRebar(J)%SIG=Elem(I)%SRebar(J)%Stress end do end dodo I=1, GValue%NRebarRebar(I)%SIG=Rebar(I)%Stress Rebar(I)%EPS=Rebar(I)%Strain

12、end do! 把新的參數賦給保留參數do I=1,GValue%NNode; Node0(I)=Node(I); end dodo I=1,GValue%NElem; Elem0(I)=Elem(I); end dodo I=1,GValue%NRebar; Rebar0(I)=Rebar(I); end dodo I=1,Gvalue%NLoad; Load0(I)=Load(I); end dodo I=1,GValue%NSupport; Support0(I)=Support(I); end do write(*,*) " 計算收斂 , 進行下一步計算 " II=

13、II+1GValue%FinishedStep=GValue%FinishedStep+1 ! 完成一步計算exit LOOP_Bend ifend if! 得到前 3 此誤差修正if(III=2) GValue%ForError(1)=GValue%GError if(III=3) GValue%ForError(2)=GValue%GError if(III=4) GValue%ForError(3)=GValue%GError if(III>4) thenGValue%ForError(1)=GValue%ForError(2)GValue%ForError(2)=GValue%F

14、orError(3)GValue%ForError(3)=GValue%GError end ifif(III>=30 .and. GValue%ForError(1)<GValue%ForError(2) .and. & GValue%ForError(2)<GValue%ForError(3).or.& GValue%GError>1.0d6) then ! 三次誤差持續增大! 荷載折半GValue%NStep=GValue%NStep*2if(GValue%NStep>GValue%MaxStep) then write(*,*) "最大分割步數 , 計算不收斂 "stop end ifGValue%FinishedStep=GValue%FinishedStep*2 write(*,*) "計算收斂困難 ,步長減半 "exit LOOP_B end if! 得到整體剛度矩陣并求解 ! c