數(shù)值傳熱學(xué)-CMAC算法-源程序代碼

1、using system; using system.collections.generic; using system.componentmodel; using system.data; using system.drawing; using system.text; using system.windows.forms; using system.io; using system.threading; using system.collections; using zedgraph; using zedgraph.web; using system.drawing.imaging; us

2、ing excel = microsoft.office.interop.excel; namespace cmac publicpartialclassform_camc : form public form_camc() initializecomponent(); privatevoid btn_js_click(object sender, eventargs e) int m = 41; /x 方向網(wǎng)格數(shù)int n = 41;/y 方向網(wǎng)格數(shù)/ 定義數(shù)組double , p = new double m + 1, n + 1;/ 壓力 p值double , p1 = new doub

3、le m + 1, n + 1;/ 計(jì)算過(guò)程中的 p值，用于過(guò)渡double , u = new double m + 1, n + 1;/ 速度 u,v 的值（這里存儲(chǔ)值的方式是以定位網(wǎng)格線(xiàn)來(lái)存儲(chǔ)，與我們通常用的行和列剛好相反，請(qǐng)?zhí)貏e注意！）double , v = new double m + 1, n + 1; double , u1 = new double m + 1, n + 1;/ 計(jì)算過(guò)程中的 u,v 值，用于過(guò)渡double , v1 = new double m + 1, n + 1; double , u11 = new double m + 1, n + 1;/ 記錄尖

4、括號(hào)的 u值double , v11 = new double m + 1, n + 1;/ 記錄尖括號(hào)的 v值double xx = new double m; / 網(wǎng)格 x坐標(biāo)double yy = new double n; / 網(wǎng)格 y坐標(biāo)double xp = new double m + 1;/ 壓力 p控制點(diǎn) x坐標(biāo)double yp = new double n + 1;/ 壓力 p控制點(diǎn) y坐標(biāo)double xu = new double m + 1;/ 速度 u控制點(diǎn) x坐標(biāo)double yu = new double n + 1;/ 速度 u控制點(diǎn) y坐標(biāo)double x

5、v = new double m + 1;/ 速度 v控制點(diǎn) x坐標(biāo)double yv = new double n + 1;/ 速度 v控制點(diǎn) y坐標(biāo)double , dux1 = new double m + 1, n + 1;/u 的動(dòng)量方程中對(duì)流相的第一項(xiàng)，x方向double , dux2 = new double m + 1, n + 1;/u 的動(dòng)量方程中對(duì)流相的第二項(xiàng)，y方向double , ksx1 = new double m + 1, n + 1;/u 動(dòng)量方程擴(kuò)散相的全部（兩相一起）double , ksx2 = new double m + 1, n + 1; doub

6、le , dvy1 = new double m + 1, n + 1;/v 的動(dòng)量方程中對(duì)流相的第一項(xiàng)，x方向double , dvy2 = new double m + 1, n + 1;/v 的動(dòng)量方程中對(duì)流相的第二項(xiàng)，y方向double , ksy1 = new double m + 1, n + 1;/v 動(dòng)量方程擴(kuò)散項(xiàng)的第一項(xiàng)double , ksy2 = new double m + 1, n + 1;/v 動(dòng)量方程擴(kuò)散項(xiàng)的第二項(xiàng)double ap, ae, aw, as, an;/ 壓力泊松方程系數(shù)double , b = new double m + 1, n + 1;/ 壓

7、力方程源項(xiàng)double t, dt, dx, dy;/ 定義時(shí)間、時(shí)間步長(zhǎng)、x方向步長(zhǎng)、 y方向步長(zhǎng)int i, j, int_p; int da; / 外部循環(huán)次數(shù)int k; / 求解壓力方程的循環(huán)次數(shù)double re = 100.0;/ 雷諾數(shù)double x = 1.0;/ 方腔長(zhǎng)度double y = 1.0;/ 方腔高度double errp = 0.0;/ 壓力 p 殘差double erru = 0.0; / 速度 u 殘差double errv = 0.0;/ 速度 v 殘差 t = 0;/ 時(shí)間 dt = 0.0001;/ 時(shí)間步長(zhǎng) dx = x / (m - 1);/x

8、 方向步長(zhǎng) dy = y / (n - 1);/y 方向步長(zhǎng) #region/ 初始化變量for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) pi, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) p1i, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) ui, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) vi, j = 0.0; fo

9、r (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) u1i, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) v1i, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) u11i, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) v11i, j = 0.0; for (i = 0; i m + 1; i+) for (j = 0; j n + 1

10、; j+) bi, j = 0.0; #endregion #region/ 坐標(biāo)初始化/ 網(wǎng)格坐標(biāo)初始化for (i = 0; i m; i+) xxi = dx * i;/ 網(wǎng)格 x坐標(biāo) for (j = 0; j n; j+) yyj = dy * j;/ 網(wǎng)格 y坐標(biāo) / 壓力 p控制點(diǎn)坐標(biāo)初始化 xp0 = 0.0;/ 壓力邊界點(diǎn) x坐標(biāo) xpm = 1.0; for (i = 1; i m; i+) xpi = (xxi - 1 + xxi) / 2.0;/ 壓力內(nèi)點(diǎn) x坐標(biāo) yp0 = 0.0;/ 壓力邊界點(diǎn) y坐標(biāo) ypn = 1.0; for (j = 1; j n; j+)

11、 ypj = (yyj - 1 + yyj) / 2.0;/ 壓力內(nèi)點(diǎn) y坐標(biāo) / 速度 u控制點(diǎn)坐標(biāo)初始化for (i = 0; i m + 1; i+) xui = xpi;/ 速度 u控制點(diǎn) x坐標(biāo) for (j = 0; j n + 1; j+) yuj = ypj;/ 速度 u控制點(diǎn) y坐標(biāo) / 速度 v控制點(diǎn)坐標(biāo)初始化for (i = 0; i m + 1; i+) xvi = xpi;/ 速度 v控制點(diǎn) x坐標(biāo) for (j = 0; j n + 1; j+) yvj = ypj;/ 速度 v控制點(diǎn) y坐標(biāo) #endregion #region/ 邊界條件for (j = 0;

12、j n + 1; j+) u0, j = 0.0;/ 為速度 u的w 、e 邊界賦初值。 um, j = 0.0; for (i = 0; i m + 1; i+) ui, 0 = 0.0;/ 為速度 u的n、s 邊界賦初值。 ui, n = 1.0; for (j = 0; j n + 1; j+) v0, j = 0.0;/ 為速度 v的w 、e 邊界賦初值。 vm, j = 0.0; for (i = 0; i m + 1; i+) vi, 0 = 0.0;/ 為速度 v的n、s 邊界賦初值。 vi, n = 0.0; #endregion #region/ 求解過(guò)程for (da =

13、1; da = 20000; da+)/ 外層循環(huán) for (i = 0; i m + 1; i+)/ 將u的值賦值給 u1 for (j = 0; j n + 1; j+) u1i, j = ui, j; for (i = 0; i m + 1; i+)/ 將v的值賦值給 v1 for (j = 0; j n + 1; j+) v1i, j = vi, j; #region/ 求解尖括號(hào) u / 求解 u的動(dòng)量方程中對(duì)流相的第一項(xiàng)，x方向， dux1 （非鄰邊界的內(nèi)節(jié)點(diǎn)）for (i = 2; i m - 1; i+) for (j = 1; j n; j+) dux1i, j = ui,

14、j * (ui + 1, j - ui - 1, j) / (2.0 * dx); / 求解邊界鄰點(diǎn)for (j = 1; j n; j+) / 左邊界鄰點(diǎn) dux11, j = u1, j * (-4.0 * u0, j + 3.0 * u1, j + u2, j) / (3.0 * dx); / 右邊界鄰點(diǎn) dux1m - 1, j = um - 1, j * (4.0 * um, j - 3.0 * um - 1, j - um - 2, j) / (3.0 * dx); / 求解 u的動(dòng)量方程中對(duì)流相的第二項(xiàng)，y方向， dux2 （非鄰邊界的內(nèi)節(jié)點(diǎn)）for (j = 2; j n -

15、1; j+) for (i = 1; i m; i+) dux2i, j = vi, j * (ui, j + 1 - ui, j - 1) / (2.0 * dy); / 求解邊界鄰點(diǎn)for (i = 1; i m; i+) / 下邊界節(jié)點(diǎn) dux2i, 1 = vi, 1 * (-4.0 * ui, 0 + 3.0 * ui, 1 + ui, 2) / (3.0 * dy); / 上邊界節(jié)點(diǎn) dux2i, n - 1 = vi, n - 1 * (4.0 * ui, n - 3.0 * ui, n - 1 - ui, n - 2) / (3.0 * dy); / 求解 u動(dòng)量方程，擴(kuò)散相第

16、一項(xiàng)，x方向， ksx1 （非邊界鄰點(diǎn)）for (i = 2; i m - 1; i+) for (j = 1; j n; j+) ksx1i, j = (ui + 1, j - 2.0 * ui, j + ui - 1, j) / (dx * dx); / 求解邊界鄰點(diǎn)for (j = 1; j n; j+) / 左邊界鄰點(diǎn) ksx11, j = (2.0 * u0, j + u2, j - 3.0 * u1, j) / (3.0 / 4.0 * (dx * dx); / 右邊界鄰點(diǎn) ksx1m - 1, j = (2.0 * um, j + um - 2, j - 3.0 * um - 1

17、, j) / (3.0 / 4.0 * (dx * dx); / 求解 u動(dòng)量方程，擴(kuò)散相第二項(xiàng)，y方向， ksx2 （非邊界鄰點(diǎn)）for (j = 2; j n - 1; j+) for (i = 1; i m; i+) ksx2i, j = (ui, j + 1 - 2.0 * ui, j + ui, j - 1) / (dy * dy); / 求解邊界鄰點(diǎn)for (i = 1; i m; i+) / 下邊界鄰點(diǎn) ksx2i, 1 = (2.0 * ui, 0 + ui, 2 - 3.0 * ui, 1) / (3.0 / 4.0 * (dy * dy); / 上邊界鄰點(diǎn) ksx2i, n

18、 - 1 = (2.0 * ui, n + ui, n - 2 - 3.0 * ui, n - 1) / (3.0 / 4.0 * (dy * dy); / 求尖括號(hào) u的值（內(nèi)點(diǎn)）for (i = 1; i m; i+) for (j = 1; j n; j+) u11i, j = ui, j - dt * (dux1i, j + dux2i, j) + (dt / re) * (ksx1i, j + ksx2i, j);/尖括號(hào) u的值 / 求尖括號(hào) u邊界值for (j = 0; j n + 1; j+) u110, j = u0, j;/ 左邊界值 u11m, j = um, j;/

19、右邊界值 for (i = 0; i m + 1; i+) u11i, 0 = ui, 0;/ 下邊界值 u11i, n = ui, n;/ 上邊界值 #endregion #region/ 求解尖括號(hào) v / 求解 v的動(dòng)量方程中對(duì)流相的第一項(xiàng)，x方向， dvy1 （非鄰邊界的內(nèi)節(jié)點(diǎn)）for (i = 2; i m - 1; i+) for (j = 1; j n; j+) dvy1i, j = ui, j * (vi + 1, j - vi - 1, j) / (2.0 * dx); / 求解邊界鄰點(diǎn)for (j = 1; j n; j+) / 左邊界鄰點(diǎn) dvy11, j = u1, j

20、 * (-4.0 * v0, j + 3.0 * v1, j + v2, j) / (3.0 * dx); / 右邊界鄰點(diǎn) dvy1m - 1, j = um - 1, j * (4.0 * vm, j - 3.0 * vm - 1, j - vm - 2, j) / (3.0 * dx); / 求解 v的動(dòng)量方程中對(duì)流相的第二項(xiàng)，y方向， dvy2 （非鄰邊界的內(nèi)節(jié)點(diǎn)）for (j = 2; j n - 1; j+) for (i = 1; i m; i+) dvy2i, j = vi, j * (vi, j + 1 - vi, j - 1) / (2.0 * dy); / 求解邊界鄰點(diǎn)fo

21、r (i = 1; i m; i+) / 下邊界節(jié)點(diǎn) dvy2i, 1 = vi, 1 * (-4.0 * vi, 0 + 3.0 * vi, 1 + vi, 2) / (3.0 * dy); / 上邊界節(jié)點(diǎn) dvy2i, n - 1 = vi, n - 1 * (4.0 * vi, n - 3.0 * vi, n - 1 - vi, n - 2) / (3.0 * dy); / 求解 v動(dòng)量方程，擴(kuò)散相第一項(xiàng)，x方向， ksy1 （非邊界鄰點(diǎn)）for (i = 2; i m - 1; i+) for (j = 1; j n; j+) ksy1i, j = (vi + 1, j - 2.0 *

22、 vi, j + vi - 1, j) / (dx * dx); / 求解邊界鄰點(diǎn)for (j = 1; j n; j+) / 左邊界鄰點(diǎn) ksy11, j = (2.0 * v0, j + v2, j - 3.0 * v1, j) / (3.0 / 4.0 * (dx * dx); / 右邊界鄰點(diǎn) ksy1m - 1, j = (2.0 * vm, j + vm - 2, j - 3.0 * vm - 1, j) / (3.0 / 4.0 * (dx * dx); / 求解 v動(dòng)量方程，擴(kuò)散相第二項(xiàng)，y方向， ksy2 （非邊界鄰點(diǎn)）for (j = 2; j n - 1; j+) for

23、(i = 1; i m; i+) ksy2i, j = (vi, j + 1 - 2 * vi, j + vi, j - 1) / (dy * dy); / 求解邊界鄰點(diǎn)for (i = 1; i m; i+) / 下邊界鄰點(diǎn) ksy2i, 1 = (2.0 * vi, 0 + vi, 2 - 3.0 * vi, 1) / (3.0 / 4.0 * (dy * dy); / 上邊界鄰點(diǎn) ksy2i, n - 1 = (2.0 * vi, n + vi, n - 2 - 3.0 * vi, n - 1) / (3.0 / 4.0 * (dy * dy); / 求尖括號(hào) v的值（內(nèi)點(diǎn)）for (i

24、 = 1; i m; i+) for (j = 1; j n; j+) v11i, j = vi, j - dt * (dvy1i, j + dvy2i, j) + (dt / re) * (ksy1i, j + ksy2i, j);/尖括號(hào) u的值 / 求尖括號(hào) v邊界值for (j = 0; j n + 1; j+) v110, j = v0, j;/ 左邊界值 v11m, j = vm, j;/ 右邊界值 for (i = 0; i m + 1; i+) v11i, 0 = vi, 0;/ 下邊界值 v11i, n = vi, n;/ 上邊界值 #endregion #region/ 求

25、解壓力泊松方程/ 壓力方程求解for (k = 1; k =10000; k+) for (i = 0; i m + 1; i+)/ 將p的值賦值給 p1 for (j = 0; j n + 1; j+) p1i, j = pi, j; for (i = 1; i m; i+) for (j = 1; j 1 & j 1 & j 1 & i 1 & i m - 1) bi, j = (v11i, j - 1 + v11i, j) / (2.0 * dy) + (u11i - 1, j - u11i + 1, j) / (2.0 * dx); if (i = 1)

26、 aw = 0; bi, j = (v11i, j - 1 + v11i, j) / (2.0 * dy) - (u11i, j + u11i + 1, j) / (2.0 * dx); if (i = m - 1) ae = 0; bi, j = (v11i, j - 1 + v11i, j) / (2.0 * dy) + (u11i - 1, j + u11i, j) / (2.0 * dx); ap = ae + aw + an + as; / 求解下一時(shí)層 p值 pi, j = (ae * pi + 1, j + aw * pi - 1, j + as * pi, j - 1 + an

27、 * pi, j + 1 + bi, j) / ap; / 設(shè)置收斂條件double max = 0; for (i = 1; i m; i+) for (j = 1; j max) max = math.abs(p1i, j); errp = max;/ 前后兩次迭代壓力的最大誤差if (errp 0.00001) break; else continue ; if (k = 10000) messagebox .show( 壓力方程不收斂，請(qǐng)檢查！, 提示 , messageboxbuttons .ok, messageboxicon.exclamation); return ; #end

28、region #region/ 求解下一時(shí)層速度并判斷收斂條件/ 求解下一時(shí)層速度u for (i = 1; i m; i+) for (j = 1; j n; j+) ui, j = u11i, j + dt * (pi - 1, j - pi + 1, j) / (2.0 * dx); / 求解下一時(shí)層速度v for (i = 1; i m; i+) for (j = 1; j n; j+) vi, j = v11i, j + dt * (pi, j - 1 - pi, j + 1) / (2.0 * dy); / 設(shè)置收斂條件double max_u = 0; double max_v

29、= 0; for (i = 1; i m; i+) for (j = 1; j max_u) max_u = math.abs(u1i, j); erru = max_u;/ 前后兩次迭代速度u的最大誤差for (i = 1; i m; i+) for (j = 1; j max_v) max_v = math.abs(v1i, j); errv = max_v;/ 前后兩次迭代速度v的最大誤差if ( math.max(erru, errv) 0.00001) messagebox .show( 計(jì)算結(jié)束請(qǐng)檢查！ , 提示 , messageboxbuttons .ok, messagebo

30、xicon.exclamation); break; else continue ; if (da = 20000) messagebox .show( 速度不收斂，請(qǐng)檢查！ , 提示 , messageboxbuttons .ok, messageboxicon.exclamation); return ; #endregion / 求解邊界壓力 , 本程序采用絕熱邊界條件for (j = 0; j n + 1; j+) p0, j = p1, j; pm, j = pm - 1, j; for (i = 0; i m + 1; i+) pi, 0 = pi, 1; pi, n = pi,

31、n - 1; #endregion / 壓力場(chǎng)datatable dt = new datatable (); for (i = 0; i m + 1; i+) datacolumn newcolumn = new datacolumn(xpi.tostring(), p0, 0.gettype(); dt.columns.add(newcolumn); for (j = 0; j n + 1; j+) datarow newrow = dt.newrow(); for (i = 0; i m + 1; i+) newrowxpi.tostring() = pi, j; dt.rows.add

32、(newrow); bindingsource bindingsource = new bindingsource (); bindingsource.datasource = dt; datagridview1.datasource = bindingsource; / 速度 u datatable dt2 = new datatable (); for (i = 0; i m + 1; i+) datacolumn newcolumn = new datacolumn(xui.tostring(), u0, 0.gettype(); dt2.columns.add(newcolumn);

33、for (j = 0; j n + 1; j+) datarow newrow = dt2.newrow(); for (i = 0; i m + 1; i+) newrowxui.tostring() = ui, j; dt2.rows.add(newrow); bindingsource bindingsource2 = new bindingsource(); bindingsource2.datasource = dt2; datagridview2.datasource = bindingsource2; / 速度 v datatable dt3 = new datatable ()

34、; for (i = 0; i m + 1; i+) datacolumn newcolumn = new datacolumn(xvi.tostring(), v0, 0.gettype(); dt3.columns.add(newcolumn); for (j = 0; j n + 1; j+) datarow newrow = dt3.newrow(); for (i = 0; i m + 1; i+) newrowxvi.tostring() = vi, j; dt3.rows.add(newrow); bindingsource bindingsource3 = new bindin

35、gsource(); bindingsource3.datasource = dt3; datagridview3.datasource = bindingsource3; #region/ 將數(shù)據(jù)導(dǎo)出為文件stringbuilder sb1 = new stringbuilder();/ 儲(chǔ)存要保存的字符串 sb1.appendformat(title=p ); sb1.appendformat(rn); sb1.appendformat(variable=x,y,t ); sb1.appendformat(rn); sb1.appendformat(zone=t,i=11,j=11,c=b

36、lack ); sb1.appendformat(rn); for (i = 0; i m + 1; i+) for (j = 0; j n + 1; j+) string myprint = ; myprint = xpi + + ypj + + pi, j + rn; / 分別是 x坐標(biāo)、 y坐標(biāo)、數(shù)據(jù)點(diǎn)。 sb1.appendformat(myprint); / 數(shù)據(jù)導(dǎo)出filestream fs1 = new filestream ( d:/p.dat, filemode .create);/ 文件流對(duì)象streamwriter sw1 = new streamwriter(fs1, encoding .getencoding( gb2312);/ 寫(xiě)入器 sw1.write(sb1.tostring();/ 操作/ 關(guān)閉 sw1.close(); fs1.close(); stringbu