121 DO 27 К 2,NYP

122 DO 23 J 2,NXP

123 JM J - 1

124 JP J + 1

125 B(2,JM) = EMXd) - CCXA - CXA*U(K,JH) - 0.5*EMX{1)*DBT*RE*

126 1{U(K,JM)*U(K,JM) + V(K,JM)*V{K,JM))

127 B(3,JM) = EMX(2) + 2.*CCXA - 0.5*EMX{2)*DBT*RE*(U(K,J)*U(K,J>

128 1 + V(K,J)*V(K,J))

129 B(4,JM) = EMX{3) - CCXA + CXA*U(K,JP) - 0.5*EMX{3)*DBT*RE*

130 1(U(K,JP)*U(K,JP) + V{K,JP) V(K,JP))

131 RRU{JH) = RU(K,J)

132 RRV(JM) RV(K,J)

133 IF(ME .NE. 2)G0T0 23

134 LST = 1

135 ir(J .EQ. 2)LST = 2 13в Bd,JM) 0.

137 DO 22 L=LST,4

138 LJ = J + L - 3

139 22 B(L,JM) = B(L,JM) + CXQ(L)*U{K,LJ)*DBT

140 23 CONTINUE

141 B(2,l) 0.

142 B(4,JM) 0.

143 IF(ME .NE. 2)G0TO 25

144 С

145 С REDUCE TO TRIDIAGONAL FORM

146 С

147 DO 24 JM 3,NXP

148 JMM JM - 1

149 DUM = Bd,JM)/B(2,JMM)

150 B(2,JM) = B(2,JM) - B(3,JMM)*DUM

151 B(3,JM) B(3,JM) - B(4,JMM)*DUM

152 Bd,JM) 0.

153 24 RRU(JM) RRU(JM) - RRiJ(JMM) *DUM

154 B{1,2) 0.

155 С

156 25 CALL BANFAC(B,NXPP,1)

157 CALL BANSOL(RRU,DDU,B,NXPPД)

158 CALL BANSOL(RRV,DDV,B,NXPP,1)

159 С

160 DO 26 J = 2,NXP

161 JM = J - 1

162 RU{K,J) = DDUWM)

163 26 RV(K,J) = DDV(JM)

164 27 CONTINUE

165 С

166 С TRIDIAGONAL SYSTEMS IN THE Y-DIRECTION

167 С

168 DO 33 J = 2,NXP

169 DO 29 К = 2,NYP

170 KM = К - 1

171 KP = К + 1

172 B{2,KM) = EMYd) - CCYA - CYA*V(KM,J) - 0.5*ЕМУ (1) *DBT*RE*

173 1(U(KM,J)*U(KM,J) + V(KM,J)*V(KM,J))

174 В(ЗДМ) = EMY{2) + 2.*CCYA - 0.5*EMY(2) *DBT*RE* (U(K, J) *U(K, J,

175 1 + V(K,J)*V(K,J))

176 B(4,KM) EMY(3) - CCYA + CYA*V(KP,J) - 0.5*ЕМУ(3)*DBT*RE*

177 1(U(KP,J)*U(KP,J) + V(KP,J)*V(KP,J))

178 RRU(KM) = RU(K,J)

179 RRV(KM) = RV(K,J)

180 IF(ME .NE. 2)GOTO 29

181 LST = 1

182 IF(K .EQ. 2)LST = 2

183 В(1ДМ) = 0.

184 DO 28 L=LST,4

185 LK = К + L - 3

186 28 B(L,KM) = B(L,KM) + CYQ(L)*V(LK,J)*DBT

187 29 CONTINUE

188 B(2,l) = 0.

189 В(4ДМ) = 0.

190 IF(ME .NE. 2)GOTO 31

191 С

192 С REDUCE TO TRIDIAGONAL FORM

193 с

194 DO 30 KM = 3,NXP

195 KMM = KM - 1

196 DUM = В(1ДМ)/В(2ДММ)

197 B(2,KM) = B(2,KM) - B(3,KMM)*DUM

198 B(3,KM) = B(3,KM) - B{4,KMM)*DUM

199 B(1,KM) = 0.

200 30 RRV{KM) = RRV(KM) - RRV(KMM)*DUM

201 B(l,2) = 0.

202 С

203 31 CALL BANFAC(B.NYPP,1)

204 CALL BANS0L(RRU,DDU,B,NYPP,1)

205 CALL BANS0L{RRV,DDV,B,NYPP,1)

206 С

207 с INCREMENT U AND V

208 С

209 DO 32 К = 2,NYP

210 KM К - 1

211 U(K,J) = U(K,J) + DDU(KM)

212 32 V(K,J) = V(K,J) + DDV(KM)

213 33 CONTINUE

214 ITER = ITER + I

215 IFdTER .LE. ITMAX)GOTO 21

216 С

217 С COMPARE WITH EXACT SOLUTION

218 С

219 34 SUMU = 0.

220 SUMV = 0.

.221 DO 36 J = 2,NXP

222 DO 35 К = 2,NYP

223 DU = U(K,J) - UE(K,J)

224 DV = V(K,J) - VE(K,J)

225 SUMU = SUMU + DU*DU

226 35 SUMV = SUMV + DV*DV

227 36 CONTINUE

228 RMSUE = SQRT(SUMU/AN)

229 RHSVE = SQRT(SUMV/AN)

230 WRITE(6, 37)ITER, RMSU, RMSV,RMSUE,RMSVE

231 37 FORMATC AFTERM3, ITERATIONS. RMS-RHS=\2E10.3,

232 1 RMS-ERR=\2E10.3)

233 DO 38 К = 1,NY

234 38 WRITE{6.39}(U(K,J),J=1,NX)

235 39 FORMATC U=\7F10.4)

236 DO 40 К = 1,KY

237 40 WRITE(6,41)(V(K,J),J=1,NX)

238 41 FORMATC V=\7F10.4)

239 STOP

240 END

10А.З. TWBURG: численное решение

Алгоритм, соответствующий схеме расщепления (10.70), (10.71), включен в качестве составной части в программу TWBURG (рис. 10.13), предназначенную для получения стационарных решений уравнения (10.52) в вычислительной области - 1 л: 1, О < / /тах, ГДС тах = Я/6Х. ГраНИЧНЫС

условия Дирихле получены из точного решения (10.63), (10.64). Начальные условия получены путем линейной интерполяции граничных условий в направлении оси х.

Основные параметры, используемые в программе TWBURG описываются в табл. 10.11. Имеется возможность реализовать решение уравнений (10.70), (10.71) в трех альтернативных вариантах. В случае ME = 1 для Lx, Lxx и т. Д. вводятся трехточечные дискретизации с центральными конечными разностями и, кроме того, полагается Ьх = Ьу = 0. При ME = 2 четырехточечная дискретизация со сдвигом вверх по потоку, применявшаяся ранее в п. 9.3.2, 9.4.3 и 9.5.2, используется как для Lx, так и для Ly, Для Lx используются четыре точки, начиная с (/ - 2, k) и кончая (/+1, k)\ для LJ/ - четыре точки, начиная с (/, k - 2) и кончая (/, k+ 1). Для типового точного решения, данные о котором приводятся на рис. 10.15, вблизи линии X = 1 величина и отрицательна. Строго говоря, для вычисления Lx в этой области должны использоваться четыре точки: (/-1, /г), (/ + 2, k). Для упрощения программирования

в программе TWBURG этого не сделано. При ME = 3 параметры массового оператора бд; и 8у выбираются эмпирическим путем, чтобы получить более точное решение.

Правая часть уравнения (10.70) вычисляется с помощью подпрограммы RHSBU (рис. 10.14), а точное решение (10.63), (10.64)-с помощью подпрограммы EXBUR (рис. 6.11) при количестве элементов в массивах UE и VE, измененном на 21 X 21. В дополнение к этому в подпрограмме EXBUR (строка 10) задается значение Хо согласно (10.65). Типовая выдача результатов, полученных по программе TWBURG, приводится на рис. 10.15.

Решения, полученные тремя указанными методами при выборе параметров согласно (10.65), приводятся в табл. 10.12 и 10.13 применительно к сетке 6X6. Соответствующие ошибки решения в зависимости от х приводятся на рис. 10.16 и 10.17. Даже на столь грубой сетке все методы дают сравнимую степень точности, причем наибольшие ошибки решения возникают в окрестности большого внутреннего градиента (рис. 10.12). Исходя из результатов § 9.4 и п. 10.1.4, можно ожидать, что больший разброс данных от разных методов