📄 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 + -