newmark.f95

实现Newmark法的fortran 控制台程序

弹性地基梁计算（纽玛克法）NDJL程序

采用纽玛克（Newmark）法计算弹性地基梁具有速度快、计算简单、占机内存少（占机内存与节点数无关）等特点。NDJL程序可用于变截面、变刚度、变地基系数的任意截面形状的弹性地基梁或基桩的内力和位移计算，梁（桩）端边界条件任意。使用时应尽量利用结构对称性，并合理地选用边界条件，以减少划分的单元数，提高计算速度和精度。

该软件 采用DATA语句输入，其输入参数依次如下：

N——梁（桩）所划分的单元数；

X——梁（桩）所划分的单元长度（m）各单元长度不等时为0；

E——梁（桩）的抗弯刚度（KM m2），变截面、变刚度时为0;

K——每个单元的地基刚度（KM/m），其值为地基系数乘以单元面积，变地基系数时为0；

JJ——输出控制参数，JJ为输出时所隔的单元数；

Qi——作用在第i个点上的垂直荷载（KM）,向下为正;

Mi——作用在第i个点上的力矩值（KM m），顺时针方向为正；

Ei——第i节点梁（桩）的抗弯刚度（KM m2），只当变截面、变刚度时输入;

Ki——第i个单元的地基刚度（KM/m），只有变地基系数时输入；

Xi——第i个单元的单元长度（m）;

A$B$——已知边界条件名称，其值为“V”、“M”、“Y”、“Phi”中的任意两个字符串，其中“V”表示剪力，“M”表示弯矩，“Y”表示位移，“Phi”表示转角。

程序中所给算例为《岩土工程学报》1982年1期第13页算例。其源程序如下：

10 READ N,X,E,K,JJ:P$=”    ”:Z$=”       ”

20 DIM K(N),A(N+1),E(N),Y(1,N),Q(N),T(N),X(N),P(N),I(1),J(1)

30 PRINT TAB(25);”---原始计算数据---”:PRINT:PRINT TAB(30);”外荷载”

40 PRINT P$;”节点号  距离（m）  垂直力（KM）  弯矩（KM m）”

50 FOR I=0 TO N: READ Q(I),T(I)

60 IF E=0 THEN READ E(I):  GOTO 80

70 E(I)=E

80 IF K=0 THEN READ K(I)

90 IF X>0 THEN X(I)=X: GOTO 110

100 IF I<N READ X(I)

110 IF Q(I) OR T(I)THEN PRINT USING”&######&######&######&######”;

120 D=D+X(I): NEXT I

130 IF K THEN FOR I=1 TO N-1; K(I)=K*(X(I-1)+X(I))/2: NEXT I

140 K(0)=K*X(0)/2: K(N)=K*X(N-1)/2

150 READ A$,B$:PRINT:S=0:PRINT TAB(27);”各单元的物理参数”

160 PRINT P$;”节点号  抗弯刚度”（KM m2）  刚度系数（KM/m）  单元长度（m）”

170 FOR I=0 TO N: IF JJ>1 THEN IF I<>INT(I/JJ)*JJ THEN 220

180 PRINT USIMG “&#####&#######.##&#######.##”;P$;I;Z$;E(I);

Z$;K(I);IF I<N THEN PRINT TAB(54);USING”#####.#####”; X(I)

190 GOTO 220

220 NEXT I

230 PRINT: PRINT: PRINT TAB(28);”边界条件：”;A$;”,” B$: PRINT

240 M=-T(N):V=-Q(N): P(N)=0:P=0:A(N)=T(N)/E(N)

250 FOR K=N-1 TO 0 STEP -1: M=M+V*X(K)-T(K): A(K)=-M/E(K)

260 A=A(K+1)+T(K+1)/E(K+1):IF K=N-1 THEN B=(2*A(N)+A(K))*X(K)/6:GOTO 280

270 B=(A(K+2)*X(K+1)+(A+A(K+1))*(X(K+1)+X(K))+X(K)*A(K))/6

280 P=P+B: P(K)=P(K+1)+P*X(K):V=K(K)*P(K)+V-Q(K)

290 NENT K: I1=(ASC(A$)-77)/3: I2=(ASC(B$)-77)/3: A(4)=P(0)

300 A(1)=P+(A(1)+2*A(0)*X(0)/6:A(3)=V:A(0)=M:Q=A(I1):R=A(I2)

310 FOR J=0 TO 1:M=0: Y(J,N)=1-J:P=J:V=K(N*Y(J,N):A(N)=0

320 FOR K=N-1 TO 0 STEP -1:M=M+V*X(K):A(K)=-M/E(K)

330 IF K=N-1 THEN B=A(K)*X(K)/6: GOTO 350

340 B=(A(K+2)*X(K+1)+2*A(K+1)*(X(K+1)+X(K))+X(K)*A(K))/6

350 P=P+B:Y(J,K)=Y(J,K+1)+P*X(K): V=K(K)*Y(J,K)+V

360 NEXT K:A(4)=Y(J,0):A(1)=P+(A(1)+2*A(0))*X(0)/6:A(3)=V:A(0)=M

370 I(J)=A(I1):J(J)=A(I2):NEXT J

380 G=I(1)*J(0)-I(0)*J(1):YA=(Q*J(1)-R*(1))/G:PA=(R*I(0)-Q*J(0))/G

390 P(N)=YA:A(N)=YA*K(N)-Q(N):E(N)=-T(N):M=E(N):K=M

400 FOR I=N-1 TO 0 STEP -1

410 P(I)=Y(0,I)*YA+Y(1,I)*PA+P(I): A(I)=K(I)*P(I)+A(I+1)-Q(I)

420 E(I)=E(I+1)+A(I+1)*X(I): IF M<E(I) THEN M=E(N): G=I

430 IF K>E(I) THEN K=E(I): D=I

440 NEXT I:PRINT TAB(25);”---计算结果---”:PRINT:S=0

450 PRINT P$;”距离（m）  位移（mm）  弯矩（KM m）  剪力（KM）”

460 X=0: FOR I=0 TO N: A(I)=A(I)+Q(I):S=P(I)*1000

470 IF JJ>1 THEN IF I<>TNT(I/JJ)*JJ THEN 490

480 PRINT USING”&####.###&####.##&####.##&####.##”;P$;Z$;E(I);Z$;A(I)

490 X=X+X(I): NEXT I

500 END

510 DATA 5,6096,76352.4,2345.249,1,0,0,0,0,0,0,0,0,140.685,0,0,0,V,Phi