汉化comsol算例：地球科学模块，微孔尺度流动

cora

微孔尺度流动

简介

这一多孔介质流动算例使用了非常规的方法，在多孔介质的空隙里应用了Navier-Stokes方程。本算例来源于Santa Barbara，University of California的Arturo Keller, Maria Auset, 和Sanya Sirivithayapakorn所进行的多孔介质流动实验。为了给该算例制作几何模型他们拍摄了电子显微镜扫描图像。这一非常规的COMSOL Multiphysics微孔尺度模型，加深了我们对大体积颗粒、胶体在几何形状不定的毛细孔微表面上的运动的认识。研究者中的一些人已经将他们的COMSOL Multiphysics模拟结果发表在了Water Resources Research杂志上(Ref. 1 and Ref. 2)。

Keller, Auset, and Sirivithayapakorn的实验室实验是基于岩石薄片的电子显微镜扫描（SEM）图像来设计的。他们用来将图像上的几何形状转化成实体的复杂过程类似于蚀刻硅片。(图 2-1)。 之后，图像被转成DXF文件并输入COMSOL Multiphysics。

图 2-1: 硅片上重复纹样的电子显微镜扫描图像。图片底部文字显示微孔的最大最小尺度大约在1 μm–100 μm数量级(Ref. 1).

cora

这一问题的经典解决方法是，将微表面上的流动当做连续性的求其平均，或求解主体流动的连续体性质，而不是求其具体形状，或追踪多孔介质中的每一个粒子运动方向。主体流动的性质可以采用一个方程插入，如达西定律，从而给出整个体积的平均流量。主流近似可以给出相当好的估算，对于较大区域内的流动是足够的，但是对于晶体颗粒之间的精微流动，主流近似则不能达到封闭Navier-Stokes方程组的解析描述精度。

本算例的对象是Keller, Auset,和Sirivithayapakorn二维电镜图像中的一张，应用笛卡尔坐标系下的Navier-Stokes方程求解微孔流动的速度和压强。算例以DXF文件输入几何形状，对不需求解的区域都没有划分网格。模拟结果图包括三维图，由二维图加上高度数据而得到。对边界做了积分，用于量化流量。

cora

模型定义
整个模型大小为640 μm* 320 μm。水主流方向是从右至左。在微孔内流动为层流，而且并不进入晶界。流动的进出口压力已知。假定流动对于上下边界是对称的。我们感兴趣的首要区域是左上角坐标为（0，0）μm到右下角坐标为(581.6, −265.0) μm之间的长方形区域。
假定微孔区的流体密度一定、温度不变，则流动可以以不可压缩的Navier-Stokes方程组来描述：

其中 η 表示动力学粘性系数 (kg/(m·s)), u 表示速度 (m/s), ρ 为密度 (kg/m3), p 表示压强 (Pa), I 为单位距阵。由于该问题的尺度非常小，所以模型中无需考虑重力引起的体积力。

cora

表 2-1: 边界设定 边界类型 边界设定 值
出口 出口: 压力，无粘滞应力 p=0
进口 进口, 压力，无粘滞应力 p=p0
晶界 壁面 -
对称边界 对称 -

p0为进出口压降
表 2-2 为模型相关数据。
表 2-2:  变量 值 描述
ρ0 1000 kg/m3 流体密度
η0 0.001 kg/(m·s) 动力学粘性系数
p0 0.715 Pa 压降

cora

结果
图 2-2 为 COMSOL Multiphysics应用Navier-Stokes 方程分析得到的微尺度多孔介质孔内流动的流体速度场。在进口处最狭窄的毛细孔内流速最高，在管道出现交叉伸展、流量增长的区域流速则有所降低。
图 2-3 为速度场图，高度代表压力水头。本图清楚地表现了大孔盲端的死区。
图中揭示了最高速度出现的位置，如所预期的，出现在具有高压降的狭窄孔径内。图2-4为出口附近某一区域的局部放大图，表现了高速也可能出现在压力梯度小而周围支管多的宽管内。
图 2-5: 出口边界上的法向速度。因为速度方向是与X轴反向的，所以速度是负的。

cora

使用图形化用户界面建模
模型导航视窗
1
打开 模型导航视窗, 在 空间维度 中选择 2D.
2
在应用模块菜单中选择 地球科学模块>流体流动>不可压缩 Navier-Stokes. 点击 确定.
几何建模
1
从 DXF 文件输入几何模型. 在 文件 菜单内, 选择 输入>从文件的CAD模型 并导航至 pore_scale.dxf.
2
保持几何图形输入的默认设置. 点击 确定.
3
点击绘图工具栏了上的 强迫成实体 按钮.
4
点击绘画工具栏上的 分离对象 按钮.
5
用鼠标选择对象CO2，CO2为连通的空隙流动区域。
6
按 Ctrl+C，复制该几何对象。
7
按 Ctrl+A 选择所有几何对象. 按Delete 键将它们全部删除.
8
按 Ctrl+V粘贴被复制的对象, 在弹出粘贴位置对话框中点击 确定 。

9
在主菜单工具栏中点击 缩放至视窗大小 按钮。
本例的几何尺度是微观的。可以利用COMSOL Multiphysics的单位转换功能，在几何建模中将之转换为SI长度单位米
10
点击绘图工具栏上的 比例刻度 按钮.
11
在比例系数对话框内，x和y框内写入1e-6。点击确定。.
12
再次点击缩放至视窗大小.

cora

选项与设定
1  在 物理量 菜单中, 选择 模型设定.
2  清除 简化表达式 选项, 点击确定.
3  在 选项 菜单中, 选择 常数.
4  写入下列常数名称、表达式和描述（可选的）; 点击 确定.
名称        表达式        描述
rho0        1000[kg/m^3]        流体密度
eta0        0.001[kg/(m*s)]        动力学粘性系数
p0        0.715[Pa]        压降
物理量设定
求解域设定
在 物理量 菜单中, 选择 求解域设定. 选择求解域1, 将下表内的物质属性写入对话框; 点击 确定.
变量        求解域1
ρ        rho0
η        eta0
边界条件
在 物理量 菜单中, 选择 边界设定. 按照下表将多种边界条件写入对话框; 点击 确定.
设定        边界 1–6        边界 7–25        边界 26, 27        其它
边界类型        出口        对称边界         进口        壁面
边界条件        压力, 无粘滞应力        -        压力, 无粘滞应力        无滑动
p0        0                 p0         
网格生成
1  在网格 菜单中, 打开 自由网格参数 对话框.
2  在 全局 页面上, 点击 自定义网格尺寸 按钮, 写入以下数据:
全局网格参数        表达式
最大单元尺寸        1e-5
单元增长率        1.6
网格曲率系数        0.3
3  在 边界 页面上, 选择边界 1–6.
4  在 最大单元尺寸 编辑框内写入 3.
5  在 高级 页面上, 将 几何松弛度 设为20.
6  点击 重划网格 按钮. 网格重新生成之后, 点击 确定

cora

计算求解
点击主菜单工具栏上的 求解 按钮.
后处理与绘图
生成 图 2-2 可根据以下步骤:
1  在 后处理 菜单中, 打开 绘图参数 对话框.
2  在 通用 页面上, 选择 表面 和 箭头 复选框.
3  点击 表面 分页. 在 表面数据 页面上的 内建物理量 下拉选项中, 选择 速度场. 在 单位 编辑框内, 输入 um/s 使速度的显示单位为 μm/s.
4  点击 箭头 分页.在 求解域数据 页面上的 内建物理量 下拉选项中, 选择 速度场.
5  在箭头位置 , x 点 和 y 点 均输入 40.
6  点击 应用.
生成 图 2-3 可根据以下步骤:
1  回到 绘图参数 对话框 并清除 箭头图 的选项.
2  在 表面 页面上, 点击 高度数据 分页. 选上 高度数据 选项, 在 表达式 编辑框内写入 p/(rho0*9.82[m/s^2]. 在 单位 下拉选项中, 选择 um 来显示动压水头。  
3  点击 确定.

生成图 2-4 可根据以下步骤:
首先定义感兴趣的区域. 在3D绘图工具栏上激活 移近/远 按钮, 或使用主菜单工具栏上的 选取视窗缩放 工具.
生成 图 2-5 可根据如下步骤:
1  在 后处理 菜单中, 打开 域图参数 对话框 (因为要显示边界上的物理量).
2  点击 线/拉伸 分页. 在 边界选择 框内选择边界 1–6.
3  在 y轴数据 框内的 内建物理量 下拉选单中, 选择 x速度.
4  在 x轴数据 选区内, 点击下方的选项，并点击 表达式 按钮
5  在 x轴数据 对话框内, 在 表达式 内填写y，并将 单位 改为 um . 点击 确定.
6  点击 确定 ， 域图参数 对话框关闭并生成图像.

基于Stokes法则跟踪粒子轨迹，可以根据以下步骤.
1  在 选项 菜单中, 打开 常数 对话框.
2  加入下列常数; 然后点击 确定.
名称 表达式 描述
p_mass 1[ug] 粒子质量
p_rad 1[um] 粒子半径
p_vol 4/3*pi*p_rad^3 粒子体积
g 0[m/s^2] 重力加速度
如前所述，本问题的微观尺度允许我们忽略重力. 将重力加速度常数g（为0）包含进来是，为了在重力作用不可忽略的条件下进行模拟的方便。
3  在 求解 菜单下, 选择 更新求解.
4  在 后处理 菜单中, 打开 绘图参数 对话框.
5  在 表面 页面, 清空 高度数据 分页上 高度数据 的选项.
6  在 通用 页面上, 在 绘图选择 区域内选择 箭头 和 粒子追踪 ( 表面 和 几何边 复选框也保持选择的状态).
7  点击 粒子追踪 分页. 在 质量 编辑框内, 输入 p_mass.
8  在 Fx 编辑框内, 输入 -6*pi*p_rad*eta0*(partu-u).
9  在 Fy 编辑框内, 输入 (p_mass-rho0*p_vol)*g-6*pi*p_rad*eta0*(partv-v).
10  按照下表设定 起始点:
起始点   
x 600e-6*ones(1,25)
y linspace(-80e-6,-275e-6,25)

11  在 初始值 页面上, 将 初始速度 分量设为 u, v, 和 0.
12  在 线颜色 页面上, 点击 颜色 按钮. 选择黄色,然后点击 确定.
13  点击 绘图参数 对话框上的 确定

cora

传上来不容易。。。
请大家多提意见多拍砖

eudemonia

楼主辛苦！学习学习！

cora

多谢楼上捧场
希望能得到反馈意见，有多少砖都接着哪。。。
如果帮忙提出意见的人多，还会陆续放上其它算例

TBE_Legend

你最好能整合到一起放一个pdf文档，这样看起来很不方便。

cora

是呀，这样发着更麻烦
我主要是希望直接发到版面可以方便网友指出错误
既然老大这样说了，就再放一个pdf，请帮忙看看

chbo76

好，楼主辛苦了。
支持一下。

cqulyb

大家有没有用LBM做过微流的？和comsol的结果对比一下？

kinggun

如果想让人改，不妨把HTML放上来。

cora

如果想让人改，不妨把HTML放上来。
kinggun 发表于 2009-4-10 21:57

那就直接上传html文件吧
希望有修改的反馈意见！

sdwife

支持楼主，现在在做这方面的东西，可是因英语水平有限，希望楼主多努力啊！

sdwife

最主要的是每个模型前面的理论性的东西

cora

最主要的是每个模型前面的理论性的东西
sdwife 发表于 2009-4-14 14:50

同感
就是理论部分希望有专业人士校对一下，后面的操作部分其实如果用英文版的还是直接看原文好