华科微纳系统仿真大作业

1、1、用有限差分法和有线元方法把以下问题变成数值方程，并说明两种方法的异 同：Ñ2j(x, y) = 0边界条件：j(x,0) = j(x,1) = 0j(0, y) = j(1, y) = 0解：（1）有限差分法 有限差分法的基本思想是将问题求解域划分为均匀的差分网格，用有限个网格节点代替连续的求解域。基于 Taylor 级数展开等方法，把描述问题的微分方程中的 微分用网格节点上的函数值的差分来代替，从而将微分方程转化为以网格节点上函 数值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数 值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法，Ñ

2、2j(x, y) = 0该方程为拉普拉斯方程 Ñ2u = u+ uyyxx将运算符离散化，则 f ¢(x) 的近似值公式为f ¢(x) = f ( x + h) - 2 f ( x) + f ( x - h) + O(h2 )h2（1-1）由于定义域为 0 £ x £ 1,0 £ y £ 1, 将定义域划分为均匀网域，x，y 方向的步长都记为 h，则Ñ2u = u( x + h, y) + u( x - h, y) + u( x, y + h) + u( x, y - h) - 4u( x, y) + O(h2 )

3、h2将 O(h2 ) 这个高阶无穷小舍去，则（1-2）Ñ2u = u( x + h, y) + u( x - h, y) + u( x, y + h) + u( x, y - h) - 4u( x, y) = 0h2（1-3）其中在所 有的内部网格点 (x, y) = (xi , y j ), i = 2,L, n -1; j = 2,L, m -1, 精度为O(h2 ) ，用 ui, j近似表示 u(xi , y j) ，则式（1-3）可以表示为u+ u+ u+ u- 4uÑ2ui, j» i +1, ji -1, ji , j +1i , j -1i , j

4、= 0h2（1-4）这就是拉普拉斯方程的五点差分公式，将 h2 消去可得ui+1, j + ui-1, j + ui, j +1 + ui, j +1 + ui, j -1 - 4ui, j = 0（1-5）根据边界条件j(x,0) = j(x,1) = j(0, y) = j(1, y) = 0,将定义域划分为 5×5 网格，如下图所示：最终得到由 9 个方程组成的 AP=B 线性方程组，表示为：- 4 p1 + p2p1 - 4 p2 + p3p2 - 4 p3+ p4+ p5= 0= 0+ p6= 0p1- 4 p4 + p5p2+ p4 - 4 p5 + p6+ p7= 0+

5、 p8= 0p3+ p5 - 4 p6p4- 4 p7 + p8+ p9= 0= 0p5+ p7 - 4 p8 + p9= 0p6+ p8 - 4 p9 = 0通过高斯消去法可以得到解向量P = p1 , p2 ,L, p9 （2）有限元方法 使用有限元发的计算流程为：A、求解区域离散化；B、 构造插值函数形成分段光滑的坐标函数系；2C、 用 Ritz 方法求解微分方程。 对j(x, y) 构造函数1 éæ ¶ jö2æ ¶ jö ùPj(x, y) = òòêçç

6、÷÷2 êëè ¶x ø+ çç÷÷è ¶y øúdxdyúû首先将整个区域离散为三角形的子区域如图 1，三角形微小子区域中的值由三 角形节点值的插值结果表示，即图 1 求解域离散为 n 三角形j(x, y) = Niji + N jjj + Nkjk其中 Ni , N j , Nk 为三角形的节点插值函数，ji ,jj ,jk 为函数j(x, y) 在节点 i, j, k处的函数值。可知每一个三角形子区域中泛函可由节点插值函

7、数和节点函数值表示，那么对 于整个求解区域的泛函表达式为,Pj(x, y)= P(j1,j2,×××,jn )由变分原理可知，¶P = ¶P= ××× = ¶P = 0¶j1¶j2¶jn计算这 n 个方程即可得到整个求解区域的值。（3）有限差分法与有限元法的异同相同的地方：两种方法都使用了离散化的思想，对求解域进行离散。并且都是 将求解域划分成有限个网格进行近似求解。不同的地方：有限差分法是利用级数的概念将连续函数离散化，正如高等数学 上所学的连续函数用泰勒级数表达一样，网格

8、上的结点就是级数中的一个取值点， 这样以级数和的形式求得最终的解，这个解是近似解，其余项就是误差。有限元法 是利用插值原理对求域进行近似求解，将求解域划分网格，每个网格看作一个单元 进行求解，这样可以得到若干有限个单元的解，这些解的集和构成整体函数的解。 就是说每个单元一个解，这些解分布在整个求解域上，构成不同区域解的变化，如 力的变化，温度的变化，这样就可以宏观上看到在不同点上不同的值了。其中具体的表现在以下几个方面：1、有限元单元的形状没有限制，可以做到单元形状与边界重合，处理无规边 界的问题；2、有限元统一对待区域内的节点与边界上的节点，因此节点的计算精度总体 上协调；而有限差分必须分别

9、处理微分方程与定解条件；3、有限元法是各种单元可以混合使用，所以写不出方程也可以求解；差分法 划分的网格是规则的，对方程进行离散化，就是用很多个差分代替微分。2、运 行 COMSOL 软 件 ， 打 开 ModelLibrary- à Thermal Structure Interaction-àMicrorobot 3D 的示例，详细看懂说明书该问题的分析，用自己熟 悉的 3D 工具进行建模，用 COMSOL 进行仿真，谈谈 COMSOL 多场耦合的体会。 解：（1）仿真过程：微型机器人的腿部使用的是由聚酰亚胺构成的 V 型凹槽结构。由于聚酰亚胺 具有很高的热膨胀系数，当

10、其受到加热时，腿部的变形程度也会相对较低。Microrobot leg 模型用到了热转移模块和结构力学模块，需要对它进行瞬态热 转移分析和准静态热变形分析。首先在 SolidWorks 中进行建模。图 1 SolidWorks 建模效果图 上图为腿部组装完成的效果图。在实际建模过程中，没有考虑由 Al 和 Si2O 构成的组合层。下图 2 为细节部分。图 2 腿部细节然后将 SolidWorks 模型导入到 Comsol 中，得到如下的模型：图 3 导入到 Comsol 后的效果 实验过程中，边界的冷却效应对仿真结果有着重大的影响。根据经验，很大一部分的热量是从腿部本身，也就是硅结构上散走；而

11、在腿部的末端也会产生少量的 散热现象。热量直接在腿部进行传导。由于各种材料的热膨胀系数的差异，导致腿 部的弯曲变形。下图 4 为仿真过程中，在 20ms 时的弯曲程度和边界上的温度分布。从中可以得知，仿真过程中温度变化量的最大值为 304.7K、最小值为 274.7K。图 4 在 20ms 时的温度分布示意图下图 5 为在腿部末端，时间与其变形量的关系曲线。其中，最大的位移变形量 接近 3.5X10-5m，时间点为 0.01s。可以得知，当加热温度上升时，变形量也会增加； 而当加热停止一段时间后，又开始向初始状态恢复。（2）体会：图 5 腿部末端处的变形量曲线多场耦合能最大程度地体现分析对象的

12、真实工况，在真实环境下，研究对象往 往同时受多个场的作用。体现多场耦合也有很多途径和方法，大体来说可分为多场 直接耦合和多个场顺序耦合。一般情况下，这两种耦合都采用有限元方法来计算。 COMSOL 的多场耦合是基于加权残值法的，且是顺序耦合，因为 comsol 的各个模 块的单元都是独立的，在计算各个场变量的时候是由顺序的。而直接耦合，由于每个单元同时包括各个场的自由度，所以它的即时和同时的，一次计算的各场变量是 没有顺序的，这也是直接耦合和顺序耦合的区别。comsol 的多场耦合，虽然也是有顺序性，但它是通过耦合变量来求解各个场 的反演方程和各个积分方程，耦合变量会在设定的模型中与各个独立场

13、变量同时求 解，comsol 的多场耦合在求解一个场的时候,耦合变量会调用下一场去反算,通过引 入耦合变量耦合变量来体现双向耦合。从原理上看，直接耦合最接近实际情况，因 为采用耦合单元，它的每步计算同时地，无中间条件地体现了多场作用和影响。个 人认为 comsol 的多场耦合结果的可靠和合理程度，应该居于直接耦合和顺序耦合 之间。comsol 采用的耦合方法是有很大通用性和开放性的，只要有场的控制方程， 就能计算场的耦合。因而其使用范围较广。3、微型层流混合器（1）在 COMSOL MODEL 中 MEMS-Microfluidics-lamella mixer;（2）用熟悉的 3D 工具建模

14、，用 COMSOL 仿真；（3）设想一种更有效的微型流体混合的结构；（4）试试用其他的软件进行仿真并比较说明 解：（2）仿真过程： 在宏观领域，流体的混合一般是通过机械激励或者湍流模式实现的；而在微观领域，这些方法不再适用。本模型展示的是微观领域的薄片层流混合器，采用的是 分层组件来提高混合效率。此模型分析的是流体的稳态和对流、扩散情况。流体通过微型弯曲通道进入混 合腔内。首先在 SolidWorks 中进行建模。图 6 SolidWorks 建模效果图然后将 SolidWorks 模型导入到 Comsol 中，得到如下的模型：图 7 导入到 Comsol 后的效果下图 8 为混合器内流体流动

15、状况仿真效果图，可以看出内部颜色的渐变效果， 这表明着层流的抛物线流动。其中流线并没有发生旋转等状况，而且在流动方向上 只是发生了轻微的改变。从下图中还可以得出，微型通道内部的最大流体速度发生 在边界处，值大致为 1.32mm/s；而混合腔内，可以很清楚的看到流速在水平方向 上的变化，最大的流速值大致为 0.5 mm/s。图 8 流体流速分布图下图 9 展示的是模型边线上的流体浓度分布效果。在微型通道的进口端，流体 浓度基本可以维持在固定值，0 或者 50mol/m3。当流体进入容腔后，开始混合。在 容腔的入口处，浓度差异最大，随着流体深入到内部，其差异性也慢慢降低。由于 流速的影响，在容腔内

16、部，边界处的混合效果要比中心处的好。但是混合的效果并 不是很理想，可以尝试降低流速、增长容腔尺寸、或者改变结构来改善混合效果。图 9 流体浓度分布图下图 10 展示的是在混合容腔中心部位的 Z 轴方向上的浓度曲线。从中可以得 知，在微型管道附近，其变化很大；在容腔的尾端，呈现出扁平状的 S 型曲线变化。 在容腔的上下两边呈现出对称趋势，在中心处幅值大约在 25 mol/m3。图 10 混合容腔中心部位的 Z 轴方向上的浓度曲线（3）设想的微型流体混合器的结构： 从仿真的结果可以知道，混合的效果并不是很理想，可以采用增大容腔长度的方式来改善混合效果。即：现行的方案为：0.2mmX0.08mm，可

17、以将长度改为 0.3mm，使流体混合的距离变长。 同时，现在仿真的方案，其结构还是相对比较简单的，可以对其结构进行改善。 方案 1：可以增加微型管道的弯曲度，进而可以降低流体的流速，以提高在容腔内的混合效果。方案 2：可以将微型管道的形状改造成蜿蜒、曲折式的。流体在这种微型管道 中流动时，能够使得处于层流状态的流体产生“无序”的流动方式。而这种无序的 流动方式，使得两流体在接触时，接触面能够被扭曲和拉长，从而更大程度上提高 混合的效果。方案 3：可以在微型管道的底部加工一些不对称、相互错开的人字形结构。这 种倾斜的、不对称的结构会对低雷诺数的粘性流体产生各项异性的流阻，使得流畅 产生横向的速度

18、分量。这样整个流场就会产生旋转形的流线，当流体混合时，可以 增强混合效果。（4）与其他仿真软件比较：ANSYS 流体单元能进行流体动力学分析，分析类型可以为瞬态或稳态。分析 结果可以是每个节点的压力和通过每个单元的流率。并且可以利用后处理功能产生 压力、流率和温度分布的图形显示。另外，还可以使用三维表面效应单元和热流 管单元模拟结构的流体绕流并包括对流换热效应。本文尝试采用，ansys 进行分析。首先将模型导入到 ansys 中，并对其进行设置，然后网格化后进行求解，其分 析结果如下：图 11 ansys 仿真结果在单纯的流体分析领域，ansys 也能产生一定的结果，但是其操作过程明显要 比 comsol 来的复杂。而且 ansys 的分析结果没有 comsol 来的好。总体感觉 comsol擅长于做流体分析，而 ansys 擅长于做结构。Comsol 更加容易上手，不过如果边界 条件设定的有问题，可能求解会发生较大的偏差。课程感受经过一学期的微纳系统仿真的课程，首先对微纳系统方面有了更深入的了解。 对仿真的过程比如前处理，求解，后处理的三个过程有了一定的认识，拓宽了自己 的视野，明白在微观条