?? pre_statically_load.f90
字號:
PROGRAM Pre_Statically_Load_Grads !面內沖擊載荷加筋板非線性屈曲,預靜載,雙樣條力法,(Galerkin),梯度法;
DOUBLE PRECISION,PARAMETER::PI=3.14159265 !定義π
INTEGER,PARAMETER::MK=2,MMK=1 !MK筋的根數;MMK筋間樣條點的個數,為奇數;
INTEGER,PARAMETER::N=(MMK+1)*MK+MMK,M=5 !N為y方向內點數,M為x方向內點數,且均為奇數;
DOUBLE PRECISION P(N,N),P0(0:N+1,N),P1(N,N),P10(0:N+1,N),P2(N,N),P20(0:N+1,N),&
P4(N,N),P40(0:N+1,N)
DOUBLE PRECISION F(M,M),F0(0:M+1,M),F1(M,M),F10(0:M+1,M),F2(M,M),F20(0:M+1,M),&
F4(M,M),F40(0:M+1,M)
DOUBLE PRECISION GM(N,N),GM0(0:N+1,N),GM1(N,N),GM10(0:N+1,N),GM2(N,N),GM20(0:N+1,N),&
GM4(N,N),GM40(0:N+1,N)
DOUBLE PRECISION CI(M,M),CI0(0:M+1,M),CI1(M,M),CI10(0:M+1,M),CI2(M,M),CI20(0:M+1,M),&
CI4(M,M),CI40(0:M+1,M)
DOUBLE PRECISION PB(MK,N),PB2(MK,N),GMB(MK,N),GMB2(MK,N)
DOUBLE PRECISION S1(N*M,N*M),S2(N*M,N*M),S3(N*M,N*M),S4(N*M,N*M),S5(N*M,N*M),&
S6(N*M,N*M),S7(N*M,N*M),S8(N*M,N*M)
DOUBLE PRECISION SS1(N*M,N*M),SS2(N*M,N*M),SS3(N*M,N*M),SS4(N*M,N*M),SS5(N*M,N*M),SS6(N*M,N*M)
DOUBLE PRECISION G1(N*M,1),G2(N*M,1),G3(N*M,1),G4(N*M,1),G5(N*M,1),G6(N*M,1),G7(N*M,1)
DOUBLE PRECISION GG1(N*M,1),GG2(N*M,1),GG3(N*M,1),GG4(N*M,1),GG5(N*M,1),GG6(N*M,1),&
GG7(N*M,1),GG8(N*M,1),GG9(N*M,1),GG10(N*M,1)
DOUBLE PRECISION GGG1(N*M,1),GGG2(N*M,1),GGG3(N*M,1)
DOUBLE PRECISION QQ0(N*M,1)
DOUBLE PRECISION U(2*N*M+1),UC(2*N*M+1),DU(2*N*M+1),UU(2*N*M+1),UW(N*M),FS(N*M,1),FT(N*M,1)
DOUBLE PRECISION M1(N*M,N*M),M1N(N*M,N*M),M2(N*M,N*M),M2N(N*M,N*M),M3(N*M,N*M)
DOUBLE PRECISION L1(N*M,1),L2(N*M,1),L3(N*M,1),L4(N*M,1),L5(N*M,1),L6(N*M,1),L7(N*M,1),&
L8(N*M,1),L9(N*M,1),L10(N*M,1)
DOUBLE PRECISION LL1(N*M,1),LL2(N*M,1),LL3(N*M,1),LL4(N*M,1),LL5(N*M,1),LL6(N*M,1),&
LL7(N*M,1),LL8(N*M,1)
DOUBLE PRECISION LLL1(N*M,1),LLL2(N*M,1)
DOUBLE PRECISION PF(N*M,N*M),PFD(N*M),PFN(N*M,N*M),CN(N*M,1),AS(N*M,1),A0(N*M,1)
DOUBLE PRECISION ASX(N*M),ASCX(N*M),ASDF(N*M)
INTEGER ASM
DOUBLE PRECISION ASC,EPS
DOUBLE PRECISION R1,R2,R3,R4
DOUBLE PRECISION A1,A2,B1,B2,C1,C2,D1,D2,A3,A4,B3,B4,C3,C4,D3,D4
DOUBLE PRECISION D,HX,HY
DOUBLE PRECISION AA,BB,HP
DOUBLE PRECISION E,V,RHO
DOUBLE PRECISION TK,HK,EK,VK,RHOK,AK,IXK,GK,JK,EEK
DOUBLE PRECISION ETA1,ETA2,TD,DT,T,PP0,Q0,PMAX,PT !PP0為面內預靜載,Q0為面內預靜載,PT為動力載荷;
DOUBLE PRECISION W,W0
DOUBLE PRECISION F5,F51,F52,F54
INTEGER MM,NN,I,J,IS(N*M),JS(N*M)
OPEN(6,FILE="INPUT.TXT",STATUS="OLD")
OPEN(8,FILE="OUTPUT.TXT",STATUS="OLD")
READ(6,*)R1,R2,R3,R4 !R1,R2為x方向彈性約束,R3,R4為y方向彈性約束;
READ(6,*)AA,BB,HP !分別為板的尺寸參數(長、寬、厚);
READ(6,*)E,V,RHO !分別為板的材料參數(楊氏模量、泊松比、密度);
READ(6,*)TD,DT,T !TD為動力載荷作用時間,DT為時間離散的步長,T為計算時間;
READ(6,*)PP0,Q0,PMAX,W0 !PP0為面內預靜載,Q0為面內預靜載,PMAX為載荷峰值,W0為初缺陷;
READ(6,*)TK,HK !筋的尺寸參數(寬、高);
READ(6,*)EK,VK,RHOK !筋的材料參數(楊氏模量、泊松比、密度);
MM=M+1 !MM為板沿X方向的等分數m;
NN=N+1 !NN為板沿Y方向的等分數n;
D=E*HP**3/12.0/(1-V**2)
HX=AA/MM !板長步長;
HY=BB/NN !板寬步長;
AK=HK*TK !筋的截面積;
GK=EK/2.0/(1+VK)
JK=HK*TK**3/3.0 !筋的扭轉慣性矩;
EEK=(HP+HK)/2.0 !筋的偏心距;
IXK=TK*HK**3/12.0+TK*HK*EEK**2 !筋的截面慣性矩;
A1=-65*R1/144/NN/(55*R1/96/NN+1)
A2=(55*R1/96/NN-1)/(55*R1/96/NN+1)
B1=-R1/144/NN/(11*R1/48/NN+1)
B2=(11*R1/48/NN-1)/(11*R1/48/NN+1)
C1=-R2/144/NN/(11*R2/48/NN+1)
C2=(11*R2/48/NN-1)/(11*R2/48/NN+1)
D1=-65*R2/144/NN/(55*R2/96/NN+1)
D2=(55*R2/96/NN-1)/(55*R2/96/NN+1)
A3=-65*R3/144/MM/(55*R3/96/MM+1)
A4=(55*R3/96/MM-1)/(55*R3/96/MM+1)
B3=-R3/144/MM/(11*R3/48/MM+1)
B4=(11*R3/48/MM-1)/(11*R3/48/MM+1)
C3=-R4/144/MM/(11*R4/48/MM+1)
C4=(11*R4/48/MM-1)/(11*R4/48/MM+1)
D3=-65*R4/144/MM/(55*R4/96/MM+1)
D4=(55*R4/96/MM-1)/(55*R4/96/MM+1)
ETA1=-4.0*PMAX/TD**2 !η1;
ETA2=4.0*PMAX/TD !η2;
!組集Ψj(yi)矩陣yi=[0,b];
DO I=0,N+1
P0(I,1)=F5(I-1)+A1*F5(I)+A2*F5(I+1)
P0(I,2)=F5(I-2)+B1*F5(I)+B2*F5(I+2)
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
P0(I,J)=F5(I-J)
ELSE
P0(I,J)=0
ENDIF
ENDDO
P0(I,N-1)=F5(I-NN+2)+C1*F5(I-NN)+C2*F5(I-NN-2)
P0(I,N)=F5(I-NN+1)+D1*F5(I-NN)+D2*F5(I-NN-1)
ENDDO
!取Ψj(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
P(I,J)=P0(I,J)
ENDDO
ENDDO
!組集Φj(xi)矩陣xi=[0,a];
DO I=0,M+1
F0(I,1)=F5(I-1)+A3*F5(I)+A4*F5(I+1)
F0(I,2)=F5(I-2)+B3*F5(I)+B4*F5(I+2)
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
F0(I,J)=F5(I-J)
ELSE
F0(I,J)=0
ENDIF
ENDDO
F0(I,M-1)=F5(I-MM+2)+C3*F5(I-MM)+C4*F5(I-MM-2)
F0(I,M)=F5(I-MM+1)+D3*F5(I-MM)+D4*F5(I-MM-1)
ENDDO
!取Φj(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
F(I,J)=F0(I,J)
ENDDO
ENDDO
!組集Γj(yi)矩陣yi=[0,b];
DO I=0,N+1
GM0(I,1)=F5(I-1)-26.0/33.0*F5(I)+F5(I+1)
GM0(I,2)=F5(I-2)-1.0/33.0*F5(I)+F5(I+2)
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
GM0(I,J)=F5(I-J)
ELSE
GM0(I,J)=0
ENDIF
ENDDO
GM0(I,N-1)=F5(I-NN+2)-1.0/33.0*F5(I-NN)+F5(I-NN-2)
GM0(I,N)=F5(I-NN+1)-26.0/33.0*F5(I-NN)+F5(I-NN-1)
ENDDO
!取Γj(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
GM(I,J)=GM0(I,J)
ENDDO
ENDDO
!組集Χj(xi)矩陣xi=[0,a];
DO I=0,M+1
CI0(I,1)=F5(I-1)-26.0/33.0*F5(I)+F5(I+1)
CI0(I,2)=F5(I-2)-1.0/33.0*F5(I)+F5(I+2)
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
CI0(I,J)=F5(I-J)
ELSE
CI0(I,J)=0
ENDIF
ENDDO
CI0(I,M-1)=F5(I-MM+2)-1.0/33.0*F5(I-MM)+F5(I-MM-2)
CI0(I,M)=F5(I-MM+1)-26.0/33.0*F5(I-MM)+F5(I-MM-1)
ENDDO
!取Χj(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
CI(I,J)=CI0(I,J)
ENDDO
ENDDO
!組集Ψj′(yi)矩陣yi=[0,b];
DO I=0,N+1
P10(I,1)=(F51(I-1)+A1*F51(I)+A2*F51(I+1))/HY
P10(I,2)=(F51(I-2)+B1*F51(I)+B2*F51(I+2))/HY
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
P10(I,J)=(F51(I-J))/HY
ELSE
P10(I,J)=0
ENDIF
ENDDO
P10(I,N-1)=(F51(I-NN+2)+C1*F51(I-NN)+C2*F51(I-NN-2))/HY
P10(I,N)=(F51(I-NN+1)+D1*F51(I-NN)+D2*F51(I-NN-1))/HY
ENDDO
!取Ψj′(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
P1(I,J)=P10(I,J)
ENDDO
ENDDO
!組集Φj′(xi)矩陣xi=[0,a];
DO I=0,M+1
F10(I,1)=(F51(I-1)+A3*F51(I)+A4*F51(I+1))/HX
F10(I,2)=(F51(I-2)+B3*F51(I)+B4*F51(I+2))/HX
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
F10(I,J)=(F51(I-J))/HX
ELSE
F10(I,J)=0
ENDIF
ENDDO
F10(I,M-1)=(F51(I-MM+2)+C3*F51(I-MM)+C4*F51(I-MM-2))/HX
F10(I,M)=(F51(I-MM+1)+D3*F51(I-MM)+D4*F51(I-MM-1))/HX
ENDDO
!取Φj′(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
F1(I,J)=F10(I,J)
ENDDO
ENDDO
!組集Γj′(yi)矩陣yi=[0,b];
DO I=0,N+1
GM10(I,1)=(F51(I-1)-26.0/33.0*F51(I)+F51(I+1))/HY
GM10(I,2)=(F51(I-2)-1.0/33.0*F51(I)+F51(I+2))/HY
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
GM10(I,J)=(F51(I-J))/HY
ELSE
GM10(I,J)=0
ENDIF
ENDDO
GM10(I,N-1)=(F51(I-NN+2)-1.0/33.0*F51(I-NN)+F51(I-NN-2))/HY
GM10(I,N)=(F51(I-NN+1)-26.0/33.0*F51(I-NN)+F51(I-NN-1))/HY
ENDDO
!取Γj′(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
GM1(I,J)=GM10(I,J)
ENDDO
ENDDO
!組集Χj′(xi)矩陣xi=[0,a];
DO I=0,M+1
CI10(I,1)=(F51(I-1)-26.0/33.0*F51(I)+F51(I+1))/HX
CI10(I,2)=(F51(I-2)-1.0/33.0*F51(I)+F51(I+2))/HX
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
CI10(I,J)=(F51(I-J))/HX
ELSE
CI10(I,J)=0
ENDIF
ENDDO
CI10(I,M-1)=(F51(I-MM+2)-1.0/33.0*F51(I-MM)+F51(I-MM-2))/HX
CI10(I,M)=(F51(I-MM+1)-26.0/33.0*F51(I-MM)+F51(I-MM-1))/HX
ENDDO
!取Χj′(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
CI1(I,J)=CI10(I,J)
ENDDO
ENDDO
!組集Ψj″(yi)矩陣yi=[0,b];
DO I=0,N+1
P20(I,1)=(F52(I-1)+A1*F52(I)+A2*F52(I+1))/HY**2
P20(I,2)=(F52(I-2)+B1*F52(I)+B2*F52(I+2))/HY**2
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
P20(I,J)=(F52(I-J))/HY**2
ELSE
P20(I,J)=0
ENDIF
ENDDO
P20(I,N-1)=(F52(I-NN+2)+C1*F52(I-NN)+C2*F52(I-NN-2))/HY**2
P20(I,N)=(F52(I-NN+1)+D1*F52(I-NN)+D2*F52(I-NN-1))/HY**2
ENDDO
!取Ψj″(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
P2(I,J)=P20(I,J)
ENDDO
ENDDO
!組集Φj″(xi)矩陣xi=[0,a];
DO I=0,M+1
F20(I,1)=(F52(I-1)+A3*F52(I)+A4*F52(I+1))/HX**2
F20(I,2)=(F52(I-2)+B3*F52(I)+B4*F52(I+2))/HX**2
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
F20(I,J)=(F52(I-J))/HX**2
ELSE
F20(I,J)=0
ENDIF
ENDDO
F20(I,M-1)=(F52(I-MM+2)+C3*F52(I-MM)+C4*F52(I-MM-2))/HX**2
F20(I,M)=(F52(I-MM+1)+D3*F52(I-MM)+D4*F52(I-MM-1))/HX**2
ENDDO
!取Φj″(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
F2(I,J)=F20(I,J)
ENDDO
ENDDO
!組集Γj″(yi)矩陣yi=[0,b];
DO I=0,N+1
GM20(I,1)=(F52(I-1)-26.0/33.0*F52(I)+F52(I+1))/HY**2
GM20(I,2)=(F52(I-2)-1.0/33.0*F52(I)+F52(I+2))/HY**2
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
GM20(I,J)=(F52(I-J))/HY**2
ELSE
GM20(I,J)=0
ENDIF
ENDDO
GM20(I,N-1)=(F52(I-NN+2)-1.0/33.0*F52(I-NN)+F52(I-NN-2))/HY**2
GM20(I,N)=(F52(I-NN+1)-26.0/33.0*F52(I-NN)+F52(I-NN-1))/HY**2
ENDDO
!取Γj″(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
GM2(I,J)=GM20(I,J)
ENDDO
ENDDO
!組集Χj″(xi)矩陣xi=[0,a];
DO I=0,M+1
CI20(I,1)=(F52(I-1)-26.0/33.0*F52(I)+F52(I+1))/HX**2
CI20(I,2)=(F52(I-2)-1.0/33.0*F52(I)+F52(I+2))/HX**2
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
CI20(I,J)=(F52(I-J))/HX**2
ELSE
CI20(I,J)=0
ENDIF
ENDDO
CI20(I,M-1)=(F52(I-MM+2)-1.0/33.0*F52(I-MM)+F52(I-MM-2))/HX**2
CI20(I,M)=(F52(I-MM+1)-26.0/33.0*F52(I-MM)+F52(I-MM-1))/HX**2
ENDDO
!取Χj″(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
CI2(I,J)=CI20(I,J)
ENDDO
ENDDO
!組集Ψj⑷(yi)矩陣yi=[0,b];
DO I=0,N+1
P40(I,1)=(F54(I-1)+A1*F54(I)+A2*F54(I+1))/HY**4
P40(I,2)=(F54(I-2)+B1*F54(I)+B2*F54(I+2))/HY**4
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
P40(I,J)=(F54(I-J))/HY**4
ELSE
P40(I,J)=0
ENDIF
ENDDO
P40(I,N-1)=(F54(I-NN+2)+C1*F54(I-NN)+C2*F54(I-NN-2))/HY**4
P40(I,N)=(F54(I-NN+1)+D1*F54(I-NN)+D2*F54(I-NN-1))/HY**4
ENDDO
!取Ψj⑷(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
P4(I,J)=P40(I,J)
ENDDO
ENDDO
!組集Φj⑷(xi)矩陣xi=[0,a];
DO I=0,M+1
F40(I,1)=(F54(I-1)+A3*F54(I)+A4*F54(I+1))/HX**4
F40(I,2)=(F54(I-2)+B3*F54(I)+B4*F54(I+2))/HX**4
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
F40(I,J)=(F54(I-J))/HX**4
ELSE
F40(I,J)=0
ENDIF
ENDDO
F40(I,M-1)=(F54(I-MM+2)+C3*F54(I-MM)+C4*F54(I-MM-2))/HX**4
F40(I,M)=(F54(I-MM+1)+D3*F54(I-MM)+D4*F54(I-MM-1))/HX**4
ENDDO
!取Φj⑷(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
F4(I,J)=F40(I,J)
ENDDO
ENDDO
!組集Γj⑷(yi)矩陣yi=[0,b];
DO I=0,N+1
GM40(I,1)=(F54(I-1)-26.0/33.0*F54(I)+F54(I+1))/HY**4
GM40(I,2)=(F54(I-2)-1.0/33.0*F54(I)+F54(I+2))/HY**4
DO J=3,N-2
IF(ABS(I-J)<=2)THEN
GM40(I,J)=(F54(I-J))/HY**4
ELSE
GM40(I,J)=0
ENDIF
ENDDO
GM40(I,N-1)=(F54(I-NN+2)-1.0/33.0*F54(I-NN)+F54(I-NN-2))/HY**4
GM40(I,N)=(F54(I-NN+1)-26.0/33.0*F54(I-NN)+F54(I-NN-1))/HY**4
ENDDO
!取Γj⑷(yi)矩陣在yi=(0,b);
DO I=1,N
DO J=1,N
GM4(I,J)=GM40(I,J)
ENDDO
ENDDO
!組集Χj⑷(xi)矩陣xi=[0,a];
DO I=0,M+1
CI40(I,1)=(F54(I-1)-26.0/33.0*F54(I)+F54(I+1))/HX**4
CI40(I,2)=(F54(I-2)-1.0/33.0*F54(I)+F54(I+2))/HX**4
DO J=3,M-2
IF(ABS(I-J)<=2)THEN
CI40(I,J)=(F54(I-J))/HX**4
ELSE
CI40(I,J)=0
ENDIF
ENDDO
CI40(I,M-1)=(F54(I-MM+2)-1.0/33.0*F54(I-MM)+F54(I-MM-2))/HX**4
CI40(I,M)=(F54(I-MM+1)-26.0/33.0*F54(I-MM)+F54(I-MM-1))/HX**4
ENDDO
!取Χj⑷(xi)矩陣在xi=(0,a);
DO I=1,M
DO J=1,M
CI4(I,J)=CI40(I,J)
ENDDO
ENDDO
!求Ψj(bk)及Ψj″(bk);求Γj(bk)及Γj″(bk);
DO K=1,MK
DO J=1,N
PB(K,J)=P((MMK+1)*K,J)
PB2(K,J)=P2((MMK+1)*K,J)
GMB(K,J)=GM((MMK+1)*K,J)
GMB2(K,J)=GM2((MMK+1)*K,J)
ENDDO
ENDDO
CALL KRONECKER1(S1,P0,P0,N,HY,F0,F0,M,HX)
CALL KRONECKER1(S2,P0,P0,N,HY,F0,F40,M,HX)
CALL KRONECKER1(S3,P0,P20,N,HY,F0,F20,M,HX)
CALL KRONECKER1(S4,P0,P40,N,HY,F0,F0,M,HX)
CALL KRONECKER1(S5,P0,P0,N,HY,F0,F20,M,HX)
CALL KRONECKER1(S6,GM0,GM0,N,HY,CI0,CI40,M,HX)
CALL KRONECKER1(S7,GM0,GM20,N,HY,CI0,CI20,M,HX)
CALL KRONECKER1(S8,GM0,GM40,N,HY,CI0,CI0,M,HX)
CALL KRONECKER2(SS1,PB,PB,N,MK,F0,F0,M,HX)
CALL KRONECKER2(SS2,PB,PB2,N,MK,F0,F20,M,HX)
CALL KRONECKER2(SS3,PB,PB,N,MK,F0,F40,M,HX)
CALL KRONECKER2(SS4,PB,PB,N,MK,F0,F20,M,HX)
CALL KRONECKER2(SS5,PB,GMB2,N,MK,F0,CI20,M,HX)
CALL KRONECKER2(SS6,PB,GMB,N,MK,F0,CI40,M,HX)
CALL KRONECKER5(QQ0,P0,N,HY,F0,M,HX)
DO J=1,N*M
DO I=1,N*M
M1(I,J)=RHO*HP*S1(I,J)+RHOK*AK*SS1(I,J)
M2(I,J)=S6(I,J)+2.0*S7(I,J)+S8(I,J)
M3(I,J)=D*(S2(I,J)+2.0*S3(I,J)+S4(I,J))
ENDDO
ENDDO
DO J=1,N*M
DO I=1,N*M
M1N(I,J)=M1(I,J)
M2N(I,J)=M2(I,J)
ENDDO
ENDDO
CALL INV(M1N,N*M,L,IS,JS)
CALL INV(M2N,N*M,L,IS,JS)
CALL KRONECKER6(PF,P,N,F,M)
DO I=1,N*M
DO J=1,N*M
PFN(I,J)=PF(I,J)
ENDDO
ENDDO
CALL INV(PFN,N*M,L,IS,JS)
!求初缺陷對應的樣條點;
DO J=1,N
DO I=(J-1)*M+1,J*M
CN(I,1)=W0*SIN(PI*J/(N+1))*SIN(PI*(I-(J-1)*M)/(M+1))
ENDDO
ENDDO
A0=MATMUL(PFN,CN)
!求定樣條點的雙樣條函數值
DO J=1,N*M
PFD(J)=PF((N+1)/2*M-(M-1)/2,J)
ENDDO
DO 11 I=1,N*M
11 WRITE(8,21)I,A0(I,1)
21 FORMAT(1X,'A0(',I3,',1)=',E15.6)
!求解預靜載對應的撓度AS;
DO I=1,N*M
ASX(I)=A0(I,1)
ENDDO
ASC=1E-6
EPS=5E-6
ASM=0
CALL NLGRAD(N*M,ASX,ASC,ASCX,ASDF,EPS,ASM)
DO I=1,N*M
AS(I,1)=ASX(I)
ENDDO
WRITE(8,20)(I,ASX(I),I=1,N*M)
20 FORMAT(1X,'ASX(',I3,')=',E15.6)
WRITE(*,200)ASM
200 FORMAT(1X,'ASM=',I6)
!求解方程主程序
U(1)=0.0
DO I=2,N*M+1
U(I)=AS(I-1,1)
ENDDO
DO I=N*M+2,2*N*M+1
U(I)=0.0
ENDDO
WRITE(8,40)
40 FORMAT(10X,'T=',13X,'W=')
CALL DIF(2*N*M+1,U,DU,N,M,MK)
DO J=1,N*M
UW(J)=U(J+1)
ENDDO
W=DOT_PRODUCT(PFD,UW)
WRITE(8,50)U(1),W
60 CALL RGKT2(2*N*M+1,DT,U,DU,UC,UU,N,M,MK)
DO J=1,N*M
UW(J)=U(J+1)
ENDDO
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -