简介北京理工大学光电学院《光电成像器件计算机辅助设计CAD》课程实验报告核心内容是基于电场高斯定律、环路定律与电势差定理配合连续超张弛迭代算法与拉格朗日插值算法用C语言实现像管内电势计算的CAD程序设计。报告完整记录了从网格划分、有限差分方程建立到超张弛迭代收敛、边界封闭处理的数值求解流程并给出了圆柱形像管内任意点电位计算、对称轴电位拉格朗日插值的实现方法。资源共1个PDF文档压缩包743KB适合光电工程、微光夜视技术方向的高年级本科生与研究生阅读也可供从事电子光学系统仿真及数值计算工作的工程技术人员参考。已有171人学习下载。通过这份实验报告读者可掌握将电磁场理论模型转化为可运行程序的方法理解CAD软件在像管电场分布计算与等势线绘制中的具体应用为像管设计、电子轨迹分析与像质评估打下工程实践基础。1. 像管电场计算的一种落地解法用C语言把拉普拉斯方程变成可迭代的CAD工具拿到一份像管设计图阴极、阳极、荧光屏电极形状不规则手写拉普拉斯方程去求解析解几乎是不可能完成的任务。解析法只对少数对称理想模型有效实验测量法又受制于探头尺寸和场扰动精度很难上去。这时候用C语言写一个数值求解器在网格上把偏微分方程离散成差分方程通过连续超松弛迭代在毫秒级逼近真实电场分布再用拉格朗日插值求任意位置的电势就是一条非常务实的工程路径。这个项目正是一个完整的最小实现代码量不算大却把有限差分、连续超松弛迭代、等势线扫描绘制、叠加定理验证串成了一条完整的计算机辅助设计工具链。它不仅算出了像管内部的电势分布还把迭代因子、残差、等势线坐标全部落盘到文本文件非常适合想搞懂“电场数值计算到底怎么落地成软件”的开发者尤其是对电磁场数值解法刚上手的研究生或工程师。2. 理论基础与整体架构从三大静电场定理到五段式模块设计2.1 静电场三大定理为什么是数值算法的“免死金牌”在着手写代码之前必须先把物理边界条件想清楚。静电场有三大基本定理这不仅是理论推导的基石更是你调试程序时判断结果对错的标尺。第一个是唯一性定理。如果给定系统中所有电极的形状、排列以及每个电极的电位那么由这些电极产生的静电场将由拉普拉斯方程唯一确定。这句话的工程含义很直接只要你的网格足够密、边界条件给得正确不管你用迭代法还是直接法求解得到的电位分布应该收敛到同一个值。如果你在程序中尝试了多种初始猜测值最后迭代出来的结果不一样那一定是你的边界封闭处理或者差分格式写错了而不是算法发散的问题。第二个是相似性定理。若各电极的电位都增大为K倍空间各点的电位也增大为K倍等势面形状不变。这个定理在验证程序时非常有用你把输入文件里所有电极电压都乘2算出来的每个网格点电压也应该乘2等势线坐标则应该完全相同。如果做不到这一点说明你的程序中存在某些“硬编码”的常数影响了结果的线性比例关系。第三个是多电极系统的电位叠加定理。当各电极形状、相对位置确定后总电场分布等于各电极单独作用其余电极接地时电场的线性叠加。这个定理是程序中“叠加定理验证模块”的核心依据。实现上就是跑n遍单电极的迭代把结果加起来再跟一次跑完所有电极的完整结果做逐点对比。这既是对算法的验证也是对边界封闭处理正确性的交叉检查。2.2 不等距五点差分格式从拉普拉斯方程到线性方程组在无空间电荷的真空中静电场电位满足拉普拉斯方程。解析求解时这是偏微分方程但数值求解时需要把它离散成差分方程。程序采用旋转对称场下的“十”字形不等距五点差分格式。假设网格点坐标为(z, r)电位为φ0相邻上下左右四个点的间距分别为h1、h2、h3、h4取绝对值则差分公式可以写成φ0 (c1*φ1 c2*φ2 c3*φ3 c4*φ4) / c0其中系数c1到c4由相邻点的间距决定c0是它们的加权和。对于轴上点r0由于径向方向的间距存在0点需要用洛必达法则处理c3、c4的计算公式。这部分是程序里最容易写错的地方因为一旦网格划分不均匀系数的计算就不能用简单的“等间距差分模板”套必须根据实际网格坐标动态计算各点的间距和权重。实现时我的习惯是单独写一个系数计算函数输入是当前网格点的四个邻居坐标输出是c0到c4。这样等间距网格和不等距网格可以共用同一套差分迭代逻辑后续如果要换网格划分方式就不需要动主迭代循环了。2.3 程序模块划分输入、初始化、迭代、绘图、验证五段式这个程序的整体架构是典型的“串行流水线”风格模块边界非常清晰非常像早期C语言工程的标准组织方式。全程序一共分成五个功能模块输入模块负责读入像管参数、初始化部分变量并原样输出参数到日志文件。初始化模块负责划分网格步长初始化电位矩阵用线性插值法封闭像管边界。超松弛迭代模块这是整个CAD程序的心脏负责对像管内部电势分布进行数值计算。绘图与输出模块进行行列扫描得到等位线坐标用EasyX插件绘制电极和等位线图形再进行拉格朗日插值绘制轴上电位变化曲线。叠加定理验证模块将每个电极单独置为有效电位其余电极为0调用迭代模块计算单位电位分布函数最后叠加并与完整计算结果对比。这种五段式架构的优势在于“可替换性”。比如你不想用EasyX绘图只需要改绘图模块其他模块完全不用动你想把迭代算法换成共轭梯度法也只需要替换超松弛迭代模块本身输入输出格式和行为完全兼容。3. 核心代码拆解SOR迭代器与等势线扫描的实现细节3.1 差分方程系数计算与边界封闭的C语言实现这里给出系数计算的简化版本。需要注意的是实际工程代码通常会预先计算好所有网格点上的系数并存入数组中避免每次迭代都重复计算因为这是典型的“内存换时间”优化策略。// 计算不等距五点差分格式的系数 // 输入当前网格点坐标 (z, r) 以及上下左右四个邻居的坐标 // 输出差分系数 cc[0] 到 cc[4]其中 cc[0] 是中心点权重分母 void calc_coeff(double z, double r, double h1, double h2, double h3, double h4, double *cc) { double c1, c2, c3, c4; // 轴向差分系数直接按等距公式计算 c1 2.0 / (h1 * (h1 h2)); c2 2.0 / (h2 * (h1 h2)); // 径向差分系数需要处理 r0 的奇异性 if (r 1e-9) { // 轴上点使用洛必达法则后的极限公式 c3 4.0 / (h3 * (h3 h4)); c4 4.0 / (h4 * (h3 h4)); } else { // 非轴上点常规公式 c3 (2.0 * r h3) / (h3 * r * (h3 h4)); c4 (2.0 * r - h4) / (h4 * r * (h3 h4)); } cc[0] c1 c2 c3 c4; cc[1] c1; cc[2] c2; cc[3] c3; cc[4] c4; }代码里的关键点在于径向系数的来源。在柱坐标下拉普拉斯算子的径向部分带有1/r项直接做差分会在r0处产生奇点。程序的处理方式是对轴上点单独推导极限公式这能保证对称轴上的电位计算精度不会因为坐标奇异性而明显劣化。实际调试时如果发现轴上电位曲线出现毛刺或不光滑先检查这一段的系数计算是否正确。边界封闭的处理同样重要。像管电极并不是完全封闭的电极之间存在空白区域这些区域的边界电位是未知的。如果不做处理迭代计算就缺少边界值无法收敛。程序采用最简单的线性插值法在两个相邻电极的边界之间按照距离比例给空白区域内的边界点赋予电位值。这样做带来的误差在实际像管设计中是可以接受的因为远离电极的边界区域对轴上电位的影响本身就很小。3.2 连续超松弛迭代核心循环从初始猜测到动态调整连续超松弛迭代的核心思想是每次遍历网格时用当前点的四个邻居更新值来计算新的电位再引入一个松弛因子ω来加速收敛。当1 ω 2时收敛过程会明显加快但ω过大会导致震荡甚至发散。程序中的迭代策略分为三个阶段第一阶段用ω1进行纯高斯-赛德尔迭代作为基准启动第二阶段用固定的ω1.375迭代12次之后每12次迭代作为一个“轮次”在轮次结束时重新计算最优松弛因子并更新全局的ω值。// 连续超松弛迭代核心循环 // phi: 电位矩阵, R: 行数, C: 列数 // eps: 迭代控制精度, max_round: 最大轮次 void sor_solve(double **phi, int R, int C, double eps, int max_round) { double omega 1.0; double max_res 0.0, avg_res 0.0; for (int rnd 1; rnd max_round; rnd) { // 每轮迭代次数前两轮特殊处理后续每轮12次 int iters_per_round (rnd 2) ? (rnd 1 ? 1 : 12) : 12; for (int k 0; k iters_per_round; k) { max_res 0.0; avg_res 0.0; int cnt 0; // 遍历所有内部网格点 for (int i 1; i R - 1; i) { for (int j 1; j C - 1; j) { double old phi[i][j]; double new_val 0.0; // 这里需要根据网格点类型选取对应的差分系数 double c0, c1, c2, c3, c4; get_coeff_at(i, j, c0, c1, c2, c3, c4); new_val (c1 * phi[i-1][j] c2 * phi[i1][j] c3 * phi[i][j-1] c4 * phi[i][j1]) / c0; // SOR 更新公式 phi[i][j] (1.0 - omega) * old omega * new_val; double res fabs(phi[i][j] - old); max_res fmax(max_res, res); avg_res res; cnt; } } avg_res / cnt; // 输出本轮迭代的残差 printf(第 %d 轮第 %d 次迭代迭代因子: %f 平均残差%f , 最大残差%f\n, rnd, k 1, omega, avg_res, max_res); // 满足精度要求则提前退出 if (max_res eps) { printf(迭代收敛\n); return; } } // 每轮迭代结束后更新松弛因子 // 实际代码中omega的计算会用到谱半径估计这里简化为递增策略 if (rnd 1) omega 1.375; else if (rnd 2) omega 1.75; else omega fmin(1.95, omega 0.05); } }这段代码最需要注意的地方有两个。第一new_val的更新必须使用“最新值”也就是当前轮次已经计算过的节点值这是高斯-赛德尔迭代的基本要求。如果你用上一轮次的旧值做更新那就变成了雅可比迭代收敛速度会大幅下降。第二omega的更新策略需要根据实际残差变化做调整。在原始实验的迭代日志中可以看到omega从1.0逐步爬到1.85左右整个过程非常平滑。如果你在调试自己代码时发现omega调整后残差反而变大说明当前问题域需要用更保守的更新策略。3.3 行列扫描绘制等势线把电位矩阵变成坐标点集迭代完成之后你手里有一整张网格的电位值。要画等势线需要把“哪几个点电位相等”这个问题转化为算法可执行的操作。程序采用的方法非常经典行扫描加列扫描然后对相邻网格点做线性插值。// 扫描并记录所有等势线点 // phi: 电位矩阵, V_level: 当前要搜索的电位值 // points: 输出点集, max_points: 点集容量 int scan_equipotential(double **phi, int R, int C, double V_level, Point *points, int max_points) { int count 0; // 行扫描沿着每一行检查水平方向穿越 for (int i 1; i R - 1; i) { for (int j 1; j C - 1; j) { if (phi[i][j] V_level phi[i][j1] V_level) { // 在线性插值求精确z坐标 double t (V_level - phi[i][j]) / (phi[i][j1] - phi[i][j]); double z j t * 1.0; // 假设网格间距为1 double r i * 1.0; if (count max_points) { points[count].z z; points[count].r r; count; } } } } // 列扫描沿着每一列检查垂直方向穿越 for (int j 1; j C - 1; j) { for (int i 1; i R - 1; i) { if (phi[i][j] V_level phi[i1][j] V_level) { // 在线性插值求精确r坐标 double t (V_level - phi[i][j]) / (phi[i1][j] - phi[i][j]); double r i t * 1.0; double z j * 1.0; if (count max_points) { points[count].z z; points[count].r r; count; } } } } // 对点集按z坐标排序 qsort(points, count, sizeof(Point), compare_z); return count; }行扫描和列扫描都会产生等势线点但单独用其中一种会漏掉那些与网格线平行的等势线段。行列结合能最大程度保证等势线的连续性。扫描结束后得到的点集并不是有序的绘图路径所以程序还需要对点集按照z坐标进行排序然后按顺序连接。这个排序步骤在原始报告中描述得很清楚先按照z坐标增大顺序排序连接一次再按照z坐标减小顺序排序连接一次主要目的是处理等势线在像管底部发生回折的情况。4. 把程序跑起来dat输入文件与EasyX绘图环境配置4.1 输入文件参数详解从电极个数到迭代精度输入文件是后缀为.dat的文本文件每一行对应一类参数。文件的组织方式几乎是“无格式”的纯序列文本参数的顺序完全固定任何一行读错了后面的数据全部错位。参数含义示例值n电极总数不包括电位为0V的阴极9δ电极厚度mm占1个步长0.4Δzi相邻电极之间的距离mm一共n-1个5, 7.1, 7.1...Ni相邻电极之间要划分的步长数10, 12, 12...Vi电极电位V21, 38, 42...Δr1电极内孔径半径mm16M1Δr1范围内等步长划分的网格数25Δr2从电极内孔边沿到封闭边界处的径向距离mm10M2Δr2范围内等步长划分的网格数16ε迭代控制精度0.00002NST输出打印空间电位时网格点间隔数1或2INS轴上电位等距插值时的步长数80ΔV扫描搜索等势线电位间隔8m特定等势线扫描电位个数任意从输入格式可以看到程序对网格划分的约束是M1 M2 50所有电极厚度相同所有电极内孔径半径相同。这些约束在原始代码里被写成了硬编码的数组大小如果直接跑实验数据没问题但想改成像管结构就需要动全局常量。4.2 开发环境配置与编译运行步骤程序使用C语言编写绘图依赖EasyX插件所以运行的推荐环境是Windows系统加Visual Studio。具体配置步骤为第一步下载并安装Visual Studio社区版第二步从官网下载EasyX库运行安装程序它会自动识别当前环境中的Visual Studio版本并完成集成第三步新建一个空项目把源码文件加入项目在项目属性中确认“字符集”设置为“使用多字节字符集”因为早期代码习惯用char*处理文件名如果是Unicode字符集可能会报类型不匹配的编译错误第四步把输入文件放到与main.c相同的目录下直接编译运行。运行时程序会在控制台窗口打印迭代过程包括每一轮次、每一次迭代的迭代因子、平均残差和最大残差就像附录一展示的那样。全部迭代完成后会弹出EasyX绘图窗口显示等势线和轴上电位曲线同时会在输入文件同目录下生成一个.txt文件里面包含了完整的迭代结果和等势线坐标。由于绘图窗口是图形界面程序退出后窗口会消失如果看不到绘图结果需要在main函数末尾加上system(pause);或者getchar();来保留窗口。4.3 避坑指南迭代不收敛、边界封闭失效与绘图黑屏现象一迭代到第3轮时残差非但不下降反而突然增大。原因松弛因子ω更新过快跳过了最优点导致迭代过程进入震荡区间。原始日志中ω从1.375跳到1.755时最大残差从0.103反弹到0.280这是典型的过松弛震荡。解决把ω的更新步长调小每次只增加0.02到0.03或者引入自适应机制当检测到最大残差连续两次上升时自动回退ω值。现象二绘制的等势线在电极之间的缺口处断开或者明显扭曲。原因边界封闭处理不正确。如果两个相邻电极的电位差很大线性插值得到的边界点电位可能不符合实际的物理衰减规律导致靠近边界的网格点电位被错误牵引。解决检查边界封闭的赋值顺序确保左右边界、上下边界各自独立处理不要让插值公式跨边界引用。另外可以尝试在电极之间增加少量过渡网格点让插值更平缓。现象三EasyX绘图窗口一闪而过根本来不及看结果。原因程序运行完main函数后进程退出图形窗口被系统自动回收。解决在main函数的绘图代码后增加getchar();或system(pause);让程序挂起等待键盘输入。如果是在Visual Studio里按F5运行也可以改用CtrlF5开始执行但不调试这样控制台窗口结束后不会自动关闭。现象四修改网格数量M1或M2后程序直接崩溃或抛出越界错误。原因原始代码的数组大小是按50x113硬编码的M1M2超过50时就会越界写入内存。解决建议把这些数组定义改成动态分配使用malloc和free或者直接定义足够大的常量比如200并加上输入参数校验在读取M1M2时如果超过上限就报错退出而不是等到运行到一半才崩溃。5. 数据处理与误差分析残差日志、叠加定理验证与实际工程意义5.1 迭代残差日志解读平均残差与最大残差的玄学迭代收敛的过程不是一条光滑下降的曲线而是带着“毛刺”的震荡衰减。查看附录一里的残差日志你会发现第3轮第1次迭代的平均残差是0.056243比第2轮最后一次迭代的0.021593还要大。这正是调整松弛因子后的“回弹”现象。松弛因子的调整会让算法暂时偏离原来的路径但随后会以更快的速度收敛。解读残差日志时重点关注两点。第一是“最大残差”它代表当前网格中最不收敛的那个点的误差大小。如果最大残差一直不降说明某个区域可能被错误的边界条件污染了比如边界封闭线上有某个点的电位值没被正确初始化会持续向周围网格“发射”错误信号。第二是“平均残差”与“最大残差”的比值。正常情况下这个比值应该在0.2到0.5之间。如果平均残差很小但最大残差很大说明整个场分布已经基本稳定但有几个孤立点存在异常这种时候通常要回溯检查那几个点的坐标位置看看是否落在了电极边缘或边界线上。5.2 叠加定理验证模块把n次SOR结果叠起来叠加定理的验证逻辑非常朴素但工程价值很高。具体实现是先将所有电极电位清零只保留第1个电极的电位运行一次完整的连续超松弛迭代模块把结果存入一个独立的叠加数组然后清零所有电极只保留第2个电极的电位再运行一次迭代把结果累加到刚才的数组里重复这个过程直到所有电极都单独处理完成。// 叠加定理验证模块的核心伪代码 // base_phi: 正常计算得到的总电场电位 // super_phi: 叠加定理计算得到的电位 // single_phi: 单电极电场电位数组 void verify_superposition(double **base_phi, double **single_phi, double **super_phi, int n_electrode, int R, int C) { // 初始化叠加数组为零 for (int i 0; i R; i) for (int j 0; j C; j) super_phi[i][j] 0.0; // 循环处理每个电极 for (int k 1; k n_electrode; k) { // 设置电极电位只保留第k个电极的电压其余电极置零 set_electrode_potential(k); // 调用SOR迭代模块 sor_solve(single_phi, R, C, EPS, MAX_ROUND); // 累加到叠加数组 for (int i 0; i R; i) for (int j 0; j C; j) super_phi[i][j] single_phi[i][j] * V[k]; } // 对比叠加结果与直接计算结果 double max_diff 0.0, avg_diff 0.0, rms_diff 0.0; for (int i 0; i R; i) { for (int j 0; j C; j) { double diff fabs(base_phi[i][j] - super_phi[i][j]); max_diff fmax(max_diff, diff); avg_diff diff; rms_diff diff * diff; } } avg_diff / (R * C); rms_diff sqrt(rms_diff / (R * C)); printf(最大绝对残差: %f, 平均残差: %f, 均方根残差: %f\n, max_diff, avg_diff, rms_diff); }这段代码里有一个隐含的“坑”set_electrode_potential(k)函数需要重新初始化电位矩阵把非当前电极的电压全部清零。如果只是把电极电位值改掉而没有清零边界条件之前轮次的残余电位会污染下一次迭代的初值导致叠加结果完全错误。我一般会在set_electrode_potential里显式执行一次全矩阵初始化把边界条件重新赋值一遍确保每次迭代都在同一套边界设置下进行。叠加定理验证最大的意义在于它能以极低成本检查整个程序的边界处理和迭代逻辑是否自洽。如果叠加验证的最大绝对残差超过预设的迭代精度一个数量级那几乎可以肯定边界封闭或者迭代更新顺序里有隐藏的bug这种bug在单一组的迭代结果里很难被肉眼发现。5.3 误差来源拆解截断误差、精度误差和网格离散误差任何一种数值计算方法的误差都不能被完全消除只能被控制在可接受范围内。这个程序里的误差来源主要有三类。第一类是计算机本身的截断误差由于CPU位数有限浮点数在存储时只能保留有限的有效位数每一次加减乘除都会丢失最末尾几个bit的精度。这类误差会随着迭代次数累积但累积速度很慢在精度要求为0.00002V的数量级下完全可以忽略。第二类是给定精度带来的截断误差迭代结束条件是前后两次残差小于预设的ε这意味着所有小于ε的电位波动都被视为“已收敛”而忽略即使它可能是一个真实存在的微弱物理信号。第三类是网格离散误差这是数值方法的核心误差来源。网格永远是离散的电场在网格点之间的连续变化被线性近似替代网格越粗误差越大。程序用行列扫描加线性插值来画等势线本身就是对离散误差的一种补救但插值阶数只有一阶在曲率大的地方比如电极附近仍然会有明显的位置偏移。6. 二次开发思路从圆柱形像管到任意旋转对称结构的改造拿到这个程序后大部分人第一反应就是“怎么改代码能自定义电极形状”。原始程序假设了圆柱形像管所有电极都是完整的环形电极内孔径半径相同厚度相同。如果要做更复杂的结构和参数可以从以下三个方向入手。第一个方向是网格分布优化。把轴向网格改成非均匀分布在电极边缘和阴极附近加密网格在其他区域适当放宽网格间距。这需要修改初始化模块中的步长数组计算逻辑并且让差分系数计算函数支持任意疏密分布的网格。你可以对比一下在相同最大网格数的条件下均匀网格和非均匀网格算出来的轴上电位曲线谁更平滑。第二个方向是算法升级。把连续超松弛迭代换成多重网格法。原始报告参考文献里也提到了多重网格法的应用它在三维像管电场计算中能明显加速收敛。实现思路是先在最细的网格上迭代几次然后限制到粗网格上继续迭代得到修正项后再延拓回细网格。这个过程循环往复能大幅度消除低频误差分量。相比单纯调ω多重网格的收敛速度提升是指数级的。第三个方向是插值算法替换。程序目前用拉格朗日三项二次插值来计算对称轴上任意点的电位优点是代码简洁、容易实现缺点是精度受限。你可以把它替换成三次样条插值在已知网格点电位和电极边界值的前提下样条插值能保证二阶导数连续轴上电位曲线会更光滑计算电子轨迹时得到的偏导数也更可信。从投入产出比来看我建议你优先做第一项改造即非均匀网格因为它的工程量最小但对计算精度和速度的提升立竿见影。网格分布优化的代码改造点集中在初始化模块里只要把等间距的循环改为按比例函数分布即可。从那以后我每次把物理模型转成代码时都会强制走一遍“网格独立性验证”——用稀疏网格和加密网格各算一遍如果结果差异在可接受范围内才敢把程序交给别人用。希望帮到你。本文还有配套的精品资源点击获取