数值传热学-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 方向网格数int n = 41;/y 方向网格数/ 定义数组double , p = new double m + 1, n + 1;/ 压力 p值double , p1 = new doub

3、le m + 1, n + 1;/ 计算过程中的 p值，用于过渡double , u = new double m + 1, n + 1;/ 速度 u,v 的值（这里存储值的方式是以定位网格线来存储，与我们通常用的行和列刚好相反，请特别注意！）double , v = new double m + 1, n + 1; double , u1 = new double m + 1, n + 1;/ 计算过程中的 u,v 值，用于过渡double , v1 = new double m + 1, n + 1; double , u11 = new double m + 1, n + 1;/ 记录尖

4、括号的 u值double , v11 = new double m + 1, n + 1;/ 记录尖括号的 v值double xx = new double m; / 网格 x坐标double yy = new double n; / 网格 y坐标double xp = new double m + 1;/ 压力 p控制点 x坐标double yp = new double n + 1;/ 压力 p控制点 y坐标double xu = new double m + 1;/ 速度 u控制点 x坐标double yu = new double n + 1;/ 速度 u控制点 y坐标double x

5、v = new double m + 1;/ 速度 v控制点 x坐标double yv = new double n + 1;/ 速度 v控制点 y坐标double , dux1 = new double m + 1, n + 1;/u 的动量方程中对流相的第一项，x方向double , dux2 = new double m + 1, n + 1;/u 的动量方程中对流相的第二项，y方向double , ksx1 = new double m + 1, n + 1;/u 动量方程扩散相的全部（两相一起）double , ksx2 = new double m + 1, n + 1; doub

6、le , dvy1 = new double m + 1, n + 1;/v 的动量方程中对流相的第一项，x方向double , dvy2 = new double m + 1, n + 1;/v 的动量方程中对流相的第二项，y方向double , ksy1 = new double m + 1, n + 1;/v 动量方程扩散项的第一项double , ksy2 = new double m + 1, n + 1;/v 动量方程扩散项的第二项double ap, ae, aw, as, an;/ 压力泊松方程系数double , b = new double m + 1, n + 1;/ 压

7、力方程源项double t, dt, dx, dy;/ 定义时间、时间步长、x方向步长、 y方向步长int i, j, int_p; int da; / 外部循环次数int k; / 求解压力方程的循环次数double re = 100.0;/ 雷诺数double x = 1.0;/ 方腔长度double y = 1.0;/ 方腔高度double errp = 0.0;/ 压力 p 残差double erru = 0.0; / 速度 u 残差double errv = 0.0;/ 速度 v 残差 t = 0;/ 时间 dt = 0.0001;/ 时间步长 dx = x / (m - 1);/x

8、 方向步长 dy = y / (n - 1);/y 方向步长 #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/ 坐标初始化/ 网格坐标初始化for (i = 0; i m; i+) xxi = dx * i;/ 网格 x坐标 for (j = 0; j n; j+) yyj = dy * j;/ 网格 y坐标 / 压力 p控制点坐标初始化 xp0 = 0.0;/ 压力边界点 x坐标 xpm = 1.0; for (i = 1; i m; i+) xpi = (xxi - 1 + xxi) / 2.0;/ 压力内点 x坐标 yp0 = 0.0;/ 压力边界点 y坐标 ypn = 1.0; for (j = 1; j n; j+)

11、 ypj = (yyj - 1 + yyj) / 2.0;/ 压力内点 y坐标 / 速度 u控制点坐标初始化for (i = 0; i m + 1; i+) xui = xpi;/ 速度 u控制点 x坐标 for (j = 0; j n + 1; j+) yuj = ypj;/ 速度 u控制点 y坐标 / 速度 v控制点坐标初始化for (i = 0; i m + 1; i+) xvi = xpi;/ 速度 v控制点 x坐标 for (j = 0; j n + 1; j+) yvj = ypj;/ 速度 v控制点 y坐标 #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/ 求解过程for (da =

13、1; da = 20000; da+)/ 外层循环 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/ 求解尖括号 u / 求解 u的动量方程中对流相的第一项，x方向， dux1 （非邻边界的内节点）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); / 求解边界邻点for (j = 1; j n; j+) / 左边界邻点 dux11, j = u1, j * (-4.0 * u0, j + 3.0 * u1, j + u2, j) / (3.0 * dx); / 右边界邻点 dux1m - 1, j = um - 1, j * (4.0 * um, j - 3.0 * um - 1, j - um - 2, j) / (3.0 * dx); / 求解 u的动量方程中对流相的第二项，y方向， dux2 （非邻边界的内节点）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); / 求解边界邻点for (i = 1; i m; i+) / 下边界节点 dux2i, 1 = vi, 1 * (-4.0 * ui, 0 + 3.0 * ui, 1 + ui, 2) / (3.0 * dy); / 上边界节点 dux2i, n - 1 = vi, n - 1 * (4.0 * ui, n - 3.0 * ui, n - 1 - ui, n - 2) / (3.0 * dy); / 求解 u动量方程，扩散相第

16、一项，x方向， ksx1 （非边界邻点）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); / 求解边界邻点for (j = 1; j n; j+) / 左边界邻点 ksx11, j = (2.0 * u0, j + u2, j - 3.0 * u1, j) / (3.0 / 4.0 * (dx * dx); / 右边界邻点 ksx1m - 1, j = (2.0 * um, j + um - 2, j - 3.0 * um - 1

17、, j) / (3.0 / 4.0 * (dx * dx); / 求解 u动量方程，扩散相第二项，y方向， ksx2 （非边界邻点）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); / 求解边界邻点for (i = 1; i m; i+) / 下边界邻点 ksx2i, 1 = (2.0 * ui, 0 + ui, 2 - 3.0 * ui, 1) / (3.0 / 4.0 * (dy * dy); / 上边界邻点 ksx2i, n

18、 - 1 = (2.0 * ui, n + ui, n - 2 - 3.0 * ui, n - 1) / (3.0 / 4.0 * (dy * dy); / 求尖括号 u的值（内点）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);/尖括号 u的值 / 求尖括号 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/ 求解尖括号 v / 求解 v的动量方程中对流相的第一项，x方向， dvy1 （非邻边界的内节点）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); / 求解边界邻点for (j = 1; j n; j+) / 左边界邻点 dvy11, j = u1, j

20、 * (-4.0 * v0, j + 3.0 * v1, j + v2, j) / (3.0 * dx); / 右边界邻点 dvy1m - 1, j = um - 1, j * (4.0 * vm, j - 3.0 * vm - 1, j - vm - 2, j) / (3.0 * dx); / 求解 v的动量方程中对流相的第二项，y方向， dvy2 （非邻边界的内节点）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); / 求解边界邻点fo

21、r (i = 1; i m; i+) / 下边界节点 dvy2i, 1 = vi, 1 * (-4.0 * vi, 0 + 3.0 * vi, 1 + vi, 2) / (3.0 * dy); / 上边界节点 dvy2i, n - 1 = vi, n - 1 * (4.0 * vi, n - 3.0 * vi, n - 1 - vi, n - 2) / (3.0 * dy); / 求解 v动量方程，扩散相第一项，x方向， ksy1 （非边界邻点）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); / 求解边界邻点for (j = 1; j n; j+) / 左边界邻点 ksy11, j = (2.0 * v0, j + v2, j - 3.0 * v1, j) / (3.0 / 4.0 * (dx * dx); / 右边界邻点 ksy1m - 1, j = (2.0 * vm, j + vm - 2, j - 3.0 * vm - 1, j) / (3.0 / 4.0 * (dx * dx); / 求解 v动量方程，扩散相第二项，y方向， ksy2 （非边界邻点）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); / 求解边界邻点for (i = 1; i m; i+) / 下边界邻点 ksy2i, 1 = (2.0 * vi, 0 + vi, 2 - 3.0 * vi, 1) / (3.0 / 4.0 * (dy * dy); / 上边界邻点 ksy2i, n - 1 = (2.0 * vi, n + vi, n - 2 - 3.0 * vi, n - 1) / (3.0 / 4.0 * (dy * dy); / 求尖括号 v的值（内点）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);/尖括号 u的值 / 求尖括号 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; / 求解下一时层 p值 pi, j = (ae * pi + 1, j + aw * pi - 1, j + as * pi, j - 1 + an

27、 * pi, j + 1 + bi, j) / ap; / 设置收敛条件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( 压力方程不收敛，请检查！, 提示 , messageboxbuttons .ok, messageboxicon.exclamation); return ; #end

28、region #region/ 求解下一时层速度并判断收敛条件/ 求解下一时层速度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); / 求解下一时层速度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); / 设置收敛条件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( 计算结束请检查！ , 提示 , messageboxbuttons .ok, messagebo

30、xicon.exclamation); break; else continue ; if (da = 20000) messagebox .show( 速度不收敛，请检查！ , 提示 , 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 / 压力场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/ 将数据导出为文件stringbuilder sb1 = new stringbuilder();/ 储存要保存的字符串 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坐标、 y坐标、数据点。 sb1.appendformat(myprint); / 数据导出filestream fs1 = new filestream ( d:/p.dat, filemode .create);/ 文件流对象streamwriter sw1 = new streamwriter(fs1, encoding .getencoding( gb2312);/ 写入器 sw1.write(sb1.tostring();/ 操作/ 关闭 sw1.close(); fs1.close(); stringbu