华科微纳系统仿真作业

1、1、用有限差分法和有线元方法把以下问题变成数值方程，并说明两种方法的异同：边界条件：解： （1）有限差分法：先将拉普拉斯方程的完整形式列出，再使用差分形成代替微分形式，得到有限差分方程，之后代入边界条件进行求解。由导数的定义，二阶导数f,(x)可表示为d2fdx2=1xdfdx|x+x-dfdx|x 1hfx+h-fxh-fx-fx-hh =fx+h-2fx+f(x-h)h2对偏导数，则有：2ux2=ux+h,y-2ux,y+u(x-h,y)h2 2uy2=ux,y+h-2ux,y+u(x,y+h)h2 代入拉普拉斯方程，则有 即：（1.3）(2) 有限元法使用有限元法的计算流程为：Ø

2、; 求解区域离散化；Ø 构造插值函数形成分段光滑的坐标函数系；Ø 用 Ritz 方法求解微分方程对构造函数首先将整个区域离散为三角形的子区域如下图， 三角形微小子区域中的值由三角形节点值的插值结果表示，即其中为三角形的节点插值函数， 为函数在节点 处的函数值。可知每一个三角形子区域中泛函可由节点插值函数和节点函数值表示， 那么对于整个求解区域的泛函表达式为由变分原理可知计算这 n 个方程即可得到整个求解区域的值(3) 有限元法与差分法的区别相同的地方：两种方法都使用了离散化的思想，对求解域进行离散。不同的地方：a) 有限元法有更好的边界适应性，有限元法的单元不限于均匀规则单

3、元，单元形状可以有一定的任意性，大小也可以不同，且单元边界可以是曲线或曲面，不同形状单元可进行组合，因此可以处理任意复杂边界的结构；b) 有限元统一对待区域内的节点与边界上的节点，因此每个节点的计算精度总体上协调；而有限差分必须分别处理微分方程与定解条件；c) 有限元法是各种单元可以混合使用，所以写不出方程也可以求解；有限差分法划分的网格是规则的，对方程进行离散化，就是用很多个差分代替微分。2. 请看Microrobot的说明文件（这个文件可以在COMSOL3.5中打开Model Library-à Thermal Structure Interaction-àMicror

4、obot 3D），详细看懂说明书该问题的分析，用自己熟悉的3D工具进行建模，用COMSOL更高版本进行仿真。解：微型机器人腿构造如下：原理：热源加热使得局部温度升高，由于所使用有机物的热膨胀系数较高，故加热有有机物体积膨胀大，V型槽上部体积膨胀高于下部，故硅腿会向下弯曲。应用如下边界条件：（1） 对于热的部分，一个高导电性的边界条件替代铝的导电层： （2.1）（2）硬弹簧条件描述电阻层： （2.2）（3）结构部分包括热膨胀的壳单元，对两个薄层进行建模，COMSOL软件自动耦合壳与固体力学单元。2.2 COMSOL建模2.2.1几何模型1. 在画图菜单选择工作平面，单击y-z按钮；2. 右键几何

5、选择导入，单击浏览选择文件microrobot_3d.dxf，单击导入；图2.2 微型机械手COMSOL几何模型2.2.2材料微型机械手包含几种不同性质的材料，故对每一种材料定义子区域设定，使用子区域表达式变量更容易进入模型的子区域性质。定义下列材料参数：热导率，k比热容，c密度，rho杨氏模量，E泊松比，nu热膨胀系数，alpha同时，定义一个新变量Temperature，可以在后处理分析时显示所有区域的温度。2.2.3固体力学a.右键固体力学，选择线弹性材料，在线弹性材料设定中，杨氏模量、泊松比、密度均选择用户定义，其值分别输入E、nu、rho；b.右键固体力学，选择固定约束，在固定约束设

6、定中，手动选择边界74,75,149,1502.2.4设置固体传热2.2.5设置壳2.2.6设置固体传热22.2.7设置电流右键电流，选择电势，命名为电势1，边界选择子菜单手动选择28,50，电势子菜单电势V0值为flc1hs(t1/s-1e-3,5e-4)-flc1hs(t1/s-1e-2,5e-4)*3e-2。2.3计算结果与分析2.3.1温度场分布云图图2.4 t=20ms时微型机械手温度场分布 图2.4显示了微型机械手温度场分布，同时也显示了其变形情况，可以看出，最大温度为305K，正如所期望的那样，微型机械手沿温度升高的方向弯曲，当热源关闭时，机械手会返回到初始位置。2.3.2总位移

7、场云图 图2.5显示了微型机械手总位移场分布，：从位移变形图可以看出，腿向下弯曲，这与实际情况一致；从根部至末端，位移逐渐增大；末端最大为 即。图2.5 微型机械手总位移场分布云图2.3.3热量随时间的变化图2.6显示了热源随时间的变化情况，热量由热电阻产生，在仿真计算20ms中的前10ms内，一个大小为30mV的电势施加在每一个热电阻上，结果所产生的热源为2×1013W/m3，相应为100mW。图2.6 热源随时间的变化情况2.3.4总位移随时间的变化图2.7 机械手尖端总位移对时间的变化图2.7显示了机械手尖端总位移对时间的变化情况，位移先增大，在10ms11ms处增大到最大值约

8、为3.4×10-5m，最后减小。3. 请学习COMSOL4.4的COMSOL Reference Manual中Deformed Geometry and Moving Mesh的使用方法，用COMSOL4.4 仿真流体与结构的相互作用（问题描述请见文件:fluid_structure_interaction.pdf）。3.1模型的建立将SolidWorks建立好的模型导入comsol中，如下图所示3.1模型的建立和导入3.2添加fluid_structure_interaction模块，并设置好材料和已知条件3.3划分网格选择自由四面体网格。3.4结果分析在建立已知的模型后，输入流

9、体流速为1.5m/s图3.2 Fluid structure 位移分布云图图3.3 Fluid structure 应力分布云图分析：从位移和应力分布云图可以看出，流体与柱状体表面作用使得柱状体弯曲，与实际情况相符；顶端点为位移最大，应力最小，位移为 即，应力为；根部位移最小，应力最大，位移为，应力为 即。4、微型层流混合器 （1）请看COMSOL4.4 中Model LibrariesMEMS Module Micromixers-Lamella mixer PDF说明书， 用SOLIDWOKS工具建模，用COMSOL仿真； （2）自己设想一种可能的更有效的微型流体混合结构； （3）试试用其

10、他的软件进行仿真并比较说明 解：4.1模型描述 模型描述了一个微机电混合器中层流流体的混合。模型中采用分层流来提升混合效果，其中有薄片状的微通道。两股待混合的流体在两层中交替出现。压力促使流体从通道的后端到达前端。然后流体进入一个大的空间-混合腔。仿真过程中，既对在混合腔体中的溶性物质的对流和扩散进行了分析，还分析了流体流动的稳态条件。微型层流混合器如图3.1所示。包括了混合器的高度为30m的一个小的垂直截面。模型开始于向左弯曲的流道的中间平面，结束于向右弯曲的流道的中间平面。在混合器中，每一个微流道都有一个一边为20m的二次曲面截面。由于所选择的几何结构，在模型中微流道的高度只有10m。为了

11、能够造出弯曲的流道形状，模型中采用了同心圆，内圆半径分别为240m，270m，300m。除了曲率的方向，流道层的几何结构和尺寸都是相同的。流体从微流道流出进入一个混合腔体。腔体长200m，宽80m。图4.1 微型层流混合器几何模型模型中分析了流体的定常流动和混合器中溶解物质的对流和扩散过程。可以采用不可压缩的Navier-Stokes方程对流道内和腔体内的流体的流动进行求解。 （1）这里表示流体密度，u = (u, v, w)是流速场，p是流体压力，I是单位对角矩阵， 体动态粘滞度，F = (fx, fy, fz )是影响流体的体积力。在这个模型中，流体采用的是水，它的 = 

12、;1000 kg/m3； = 0.001 Pa·s； F = 0，因为此处没有体积力。系统在所有六个微流道的入口处施加10 Pa的压力以驱动流体流过腔体，腔体内部压力为0。在腔体出口处，流体流速只在边界垂直的方向上有分量。对流道和混合腔体的壁，采用无滑动边界条件。但是，在垂直方向上，由于几何结构限制，需要采用对称边界条件。下面的对流扩散方程描述了流体中溶解物质的浓度。 （2）这里，c是浓度，D为扩散系数， R是反应速率。在模型中，D =10-10 m2/s，R = 0，因为这里浓度不受任何反应的影响。在向左弯曲的流道入口边界

13、上，浓度为50 mol/m3，但向左弯曲的流道内浓度则为0。在混合腔体的出口边界处，物质流体通过对流流过边界。流道壁和腔体壁都因为溶解物质而隔开，而在上部和底部的边界处，采用对称边界条件。4.2 COMSOL建模4.2.1设置模型向导4.2.2几何1图4.2 微型层流混合器COMSOL几何模型4.2.3全局定义参数1. 在模型构建器窗口，右键全局定义，选择参数；2. 在参数设定窗口，定位到参数部分，在图表中设置如下：名称表达式描述p010PaDriving pressurec050mol/m3Input concentrationD_i1e-10m2/sIsotropic diffu

14、sion coefficient4.2.4材料1. 在模型构建器窗口模型1下方右键材料，选择材料浏览器；2. 在材料浏览器窗口定位到材料部分，选择基本材料>Water，liquid；3. 右键并选择增加材料至模型；4.2.5设置层流4.2.6设置稀物质传递4.2.7网格1图4.3 微型层流混合器网格划分图4.2.8求解1在模型构建器窗口，右键模型1>求解1，选择计算。4.3计算结果与分析4.3.1速度场（spf）1. 在模型构建器窗口，单击速度场（spf）>切面1；2. 在切面设定窗口，定位到平面数据部分，在面子菜单选择zx-平面；3. 在模型构建器窗口，右键速度场（spf）

15、，选择流线，命名为流线1；4. 在流线1设定窗口，定位到流线定位部分，在定位子菜单选择大小控制，最小距离设为0.01，最大距离设为0.02；5. 单击绘图按钮。图4.4 微型层流混合器切面速度场与流线图 如图3.4所示，详细描述了微型层流混合器流体的流动及其速度场，可以看出，沿着切面速度大小逐渐变化，这表明是层状抛物线流动；流线无法展示漩涡，且在流动方向只有很小的变化。图3.4也显示了微型层流混合器中最大流速发生在模型的对称边界处，在混合腔里面，只在沿x方向可见一个流速梯度；顶部与底部边界的对称条件使在z方向的速度均匀分布。在微通道中最大流速约为1.4mm/s，而在混合腔里面最大流速为0.5m

16、m/s。给出相应的长度（20m和80m），雷诺数Re分别是0.028和0.04，因此很明显是层流。4.3.2浓度图4.5 微型层流混合器切面速度场图4.6 微型层流混合器边界浓度分布图图4.6显示了微型层流混合器边界浓度分布情况，流入通道存在恒定的浓度值，分别为0或者50mol/m3,其值取决于通道。浓度在入口处有明显的分离，但是，朝向混合腔的末端浓度减小。在混合腔边缘流速较小，在边缘处的混合优于中心处。然而，这种混合不是完美的，因此，一个减小的流速，一个更长的混合腔或者一些其它提高混合质量的方式是更好的。4.3.3不同位浓度比较 图4.7显示了浓度在混合腔中心线上的分布曲线，可以看出，通道附

17、近浓度转变非常迅速，但是在接近混合腔末端的浓度分布曲线呈平坦的S形形状，在混合腔边缘中心，浓度分布曲线也呈相同的形状，但是其变化幅度约在17mol/m3到33mol/m3之间。图4.7 距离微通道不同位置（20m、80m、140m、200m）处沿混合腔中心线浓度在z轴方向的分布曲线4.4 更有效的微型流体混合的结构 2）设想结构并分析在原有模型基础上，对输入管道重新分配，左边三根，呈上二下一分布，右边三根，呈上一下二分布，并导入COMSOL对其进行分析，模型如图4.8所示。图4.8 重新设想的模型分析结果如图4.9所示。图4.9 分析结果（3）使用其他的软件进行仿真并比较说明这里用ANSYS软件CFX模块对其进行分析并比较。将使用UG建立的模型，保存为IGS格式导入ANSYS软件，采用20nodes95网格格式对模型进行网