数值密度SPH(Number Density SPH)的一些推导

本文提到的方法见于:Tartakovsky, Alexandre M., and Paul Meakin. “A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh–Taylor instability.” Journal of Computational Physics 207.2 (2005): 610-624.

普通SPH形式

密度: \rho_i=\sum_{j}m_jW_{ij}.

其中\sum_ji的核半径h内求和,W_{ij}=W(\bm{x}_i-\bm{x}_j,h)\bm{x}_i为粒子i的位置,

梯度公式(差分):

\nabla A_i=\sum_j\frac{m_j}{\rho_j}(A_j-A_i)\nabla_i W_{ij}.

梯度公式(对称):

\nabla A_i=\rho_i\sum_jm_j\left(\frac{A_i}{\rho_i^2}+\frac{A_j}{\rho_j^2}\right)\nabla_i W_{ij}.

拉普拉斯算子(一般形式):

\nabla^2 A_i=-\sum_j\frac{m_j}{\rho_j}(A_i-A_j)\frac{2\Vert\nabla_i W_{ij}\Vert}{\Vert \bm{r}_i-\bm{r}_j\Vert}

数值密度形式

数值密度定义:n_i=\rho_i/m_i.

采用数值密度的SPH算法并不区分每个粒子的质量,即认为m_i=m_j,故有

n_i=\sum_j W_{ij}.

则梯度公式(差分)为

\nabla A_i=\sum_j\frac{1}{n_j}(A_j-A_i)\nabla_i W_{ij}.

梯度公式(对称)为

\frac{m_i}{\rho_i}\nabla A_i=\sum_j m_im_j\left(\frac{A_i}{\rho_i^2}+\frac{A_j}{\rho_j^2}\right)\nabla_i W_{ij}.

注意到m_i=m_j,故有

\nabla A_i=n_i\sum_j\left(\frac{A_i}{n_i^2}+\frac{A_j}{n_j^2}\right)\nabla_i W_{ij}.

相应的拉普拉斯算子(一般形式):

\nabla^2 A_i=-\sum_j\frac{1}{n_j}(A_i-A_j)\frac{2\Vert\nabla_i W_{ij}\Vert}{\Vert \bm{r}_i-\bm{r}_j\Vert}

关于SPH通用核函数的一些推导

值的表示

对于一维情况,现假设我们有一个定义在非负实数域上的标量权值函数\phi_1(x),截断至1,即满足:\int{\phi_1(x)}\mathrm{d}x=1,且\phi_1(x)=0\forall x>1成立。

我们希望基于它设计一个关于\vec{r}的对称核函数W_1(\vec{r}),截断至h,那么这个核函数的形式应为

W_1(\vec{r})=\alpha_1 \phi_1\left(\frac{|\vec{r}|}{h}\right).

它的积分应为1,因此

\int W_1(x)\mathrm{d}x=h\alpha_1\int\phi_1\left(\frac{x}{h}\right)\mathrm{d}\frac{x}{h}=h\alpha_1=1.

\alpha_1=\frac{1}{h}

W_1(\vec{r})=\frac{1}{h} \phi_1\left(\frac{|\vec{r}|}{h}\right)

对于二维情况,我们需要先找到一个截断至1的\phi_2,满足\iint\phi_2(\vec{r})\mathrm{d}x\mathrm{d}y=1,其中\vec{r}=(x,y).

那么类似地,

\iint W_2(x)\mathrm{d}x=h^2\alpha_2\int\phi_2\left(\frac{|\vec{r}|}{h}\right)\mathrm{d}\frac{x}{h}\mathrm{d}\frac{y}{h}=h^2\alpha_2=1.

W_2(\vec{r})=\frac{1}{h^2} \phi_2\left(\frac{|\vec{r}|}{h}\right)

三维同样:

W_3(\vec{r})=\frac{1}{h^3} \phi_3\left(\frac{|\vec{r}|}{h}\right)

梯度的表示

我们知道对于任意维数均有

\nabla |\vec{r}|=\frac{\vec{r}}{|\vec{r}|}.

因此对于任意W(\vec{r})=\alpha_d\phi_d\left(\frac{|\vec{r}|}{h}\right),我们有

\nabla W=\frac{\alpha_d}{h}\frac{\vec{r}}{|\vec{r}|}\frac{\mathrm{d}\phi_d}{\mathrm{d}x}\left(\frac{|\vec{r}|}{h}\right).

Robert Saye:VIIM论文算例简介

3.10收敛性测试和数值验证

首先是每算例都做三组:

  • 固定[latex]\epsilon[/latex],改变网格尺寸h,检查[latex]h\to 0[/latex]的收敛情况。
  • 固定比例关系[latex]\epsilon=\alpha h[/latex],其中[latex]\alpha[/latex]为一常数,检查[latex]h\to 0[/latex]的收敛情况。
  • 交换极限,先算一个[latex]\epsilon\to 0^+[/latex]的内极限,然后检查[latex]h\to 0[/latex]的收敛情况。

所有测试都是在规则正方形/立方体网格上完成的。

3.10.1二维基础测试

一个圆环以固定速度向外扩展。

一个正方形以固定速度向外扩展。

一个正方形以固定速度收缩。

一个圆环在表面张力的作用下减速收缩。

这些都可以和解析解做对比。

3.10.2二维分歧点

在表面张力下,T-junction逐步演化成Y-junction:

这个没有解析解,进行收敛性分析。

3.10.3二维von Neumann-Mullins定理

使用Young-Laplace定理,沿着二维泡泡单个cell表面积一圈分,得到其面积变化率:

也就是说单个cell面积变化仅与其边数有关,而且面积随时间的变化函数总是分段线性的。

然后随机初始态算了一下,说这个符合von Neumann-Mullins定理,各不同边数cell的面积变化如下:

然后又算了一下各边数的收敛状况:

3.10.4三维分歧点

和二维的类似。

3.10.5三维von Neumann-Mullins定理

Reuleaux tetrahedron收缩。这个有解析解。

推导如下:

PARIS模拟器算例介绍

文章链接:https://doi.org/10.1016/j.cpc.2021.107849

算例部分在第四章。

4.1基础测试

4.1.1泊肃叶流动

设置一个8*8*2网格,模拟二维正方形水槽内的泊肃叶流动。参数[latex]||\nabla p||=\mu=\rho=1[/latex].在水槽左端右端设置压强边界条件。模拟至流场以1e-3精度稳定,约在1000步后,时刻1达成。使用显式viscous diffusion.上下边界处使用粘性边界条件u=0.

4.1.2二维斯托克斯流

和泊肃叶流动相同,但计算区域中心有一直径为0.5的圆盘。不进行advection以保证稳定。使用显式粘性力。如果使用隐式粘性,则收敛速度较快,但有一个量级为[latex]\tau[/latex]的误差项。除此之外还有两个测试。第一个是圆盘在重力作用下下落,第二个是设置流场入口/出口的速度边界条件而非用压差驱动。

4.1.3水滴对流

似乎是设置一个水滴让它在计算区域内移动,以证明VOF方法是正确的,能保证水滴的形状不变。

4.1.4圆柱对流

暂缺

4.1.5计算速度

大致每秒1.6e6个格子(单处理器)。

4.2表面张力波

二维计算域,两种密度不同的流体,在中间生成一个表面张力波。

4.3水滴和泡泡震荡

水滴是内部水外部空气,泡泡是外部水内部空气。一开始初始化为椭球,让它震荡变成球。用两个不同的分辨率,对比误差。

4.4雨滴下落

暂缺。特地选择了一个比较低的分辨率,因为低分辨率容易崩。

关于MAC网格不可压流体求解中投影矩阵的笔记

现假设我们有一流体格i,和与之相邻的格子j。我们研究投影矩阵A中,代表(i,j)格子对的元素(下简记为A[i,j])和i格子所对应的对角元(下简记为A[i,i])的数值情况。注意,i和j都是格子,本身坐标是一个向量,不要与矩阵坐标相混淆。

现列举j的几种情况。

情况一:j亦是流体格。此时,j格压强也是未知数,A[i,j]=-1,而A[i,i]中有一项1,如下图。(或在porous flow中,面积分数为a,那就是-a和a)

情况二:j亦是流体格,但i和j中间有一Neumann边界条件(虽然这很奇怪,我怀疑它会不会出现)。此时格子i的方程中,i-j相关联的项直接被取消,因此A[i,j]和A[i,i]的相应项都是0,如下图。

情况三:j是固体格。此时,i-j之间必有一Neumann边界条件,否则边界条件设置错误。此时,应强行将界面速度设为0,而且在矫正过程中,给j格施加一人工压强,值和i格的压强相等。从数学上,这样就使得矫正步不会改变界面处的速度。因此,i格对应的方程(也就是投影矩阵A的一行)中,i和j的压强项相互抵消,因此A[i,j]和A[i,i]中的相应项均为0,如下图。

情况四:j在计算域之外,而i和j之间没有边界条件。此时,意味着i-j之间是计算域的自由边界,一般认为j是空气,压强为常数(为方便起见,一般认为是0).因此,i格对应的方程中,有i的压强项而无j的压强项。这意味着A[i,j]=0,A[i,i]有一项1,如下图。

情况五:j是空气格。此时,求解的是一个带自由表面的不可压流体问题。如果没有LevelSet等辅助信息追踪表面,则可简单地归结为情况三,也就是A[i,j]=0,A[i,i]有一项1.而如果有LevelSet,则可根据i,j两格的phi函数,求出来一个修正系数theta,详见:http://43.135.142.98/2021/07/11/free-surface-theta/

所以仍然A[i,j]=0,A[i,i]的那一项是那个修正系数theta,如下图。

总的规律是,当且仅当i,j都是纯流体格且中间无边界条件时,A[i,j]=-1(在porous flow中是-a),其余情况下都是0.

至于A[i,i]的这一项,若有Neumann边界条件则为0。如果是水-气界面,则为theta,否则就是1 (在porous flow中是a) .

此外还需注意一个事实:如果i和j都是流体格,则i的方程中关于j的一项,和j的方程中关于i的一项是对称的,都是1(或a).因此我们可以在MAC网格的面节点上,把所有这样的A[i,j]存下来。这里注意,有的面可能两边都不是流体格,此时这个A[i,j]的值无意义,不出现在最后的泊松方程组中。在编程时可以随便赋个值,比如说0.

Robert Bridson书中关于自由表面求解的笔记

Fluid Simulation for Computer Graphics第二版,8.2节。

在求解自由表面流体时,简单地把计算域划分为液体格和气体格,会导致artifact.如果能更精确地追踪表面,可以修改一下在水-气界面附近计算压强差的方式。

根据Gibou et al. 2002的ghost fluid方法,可以把更新流体速度节点[latex]u_{i+1/2,j,k}[/latex]的方程写作

假设(i,j,k)在流体内部,而(i+1,j,k)在空气中,即[latex]\phi_{i,j,k}\leq 0<\phi_{i+1,j,k}[/latex].认为流体表面处压强p=0。我们可以算出来一个

(注意其分子分母都是负数),界面位置在[latex](i+\theta)\Delta x,j,k)[/latex]的地方。

然后,假设我们在(i+1,j,k)处有一个ghost pressure,可以求出来:

所以方程就是:

这个额外的系数可以叫做[latex]\theta[/latex].

注意[latex]\phi{i,j,k}[/latex]小的话,可能会把对角项搞得很大,但这没关系,这样更正定了,岂不美哉。

Ubuntu 20.04 Server软路由配置折腾记录

记录一下最近的一次折腾,使用一台安装了Ubuntu 20.04 Server的工控机搭建软路由,并成功配置家庭NAS的内网访问。计算机网络没好好学,全靠本能debug,感谢lcy同志的全程技术科普(

本文主要参考了https://blog.lcy.im/2017/09/14/%E4%BB%8Eubuntu-server-%E5%BC%80%E5%A7%8B%E9%85%8D%E7%BD%AE%E8%87%AA%E5%8A%A8%E5%88%86%E6%B5%81%E7%9A%84%E8%BD%AF%E8%B7%AF%E7%94%B1/

网络拓扑设计

瞎画的,是这意思就行(逃

根节点是联通的WO-67光猫路由一体机,管理子网192.168.1.x(24是子网掩码的长度,意思是前24位是这个子网固定的地址前缀),其中包括一台桥接模式运行的电力猫,提供WIFI覆盖。我们希望安装一部软路由,管理192.168.2.x子网。软路由后面桥接一台TP-LINK路由器,该路由器提供WIFI覆盖(因为软路由没有无线功能)。此外,还需要在软路由上接一台群晖NAS,并希望在整个网络中提供对NAS的访问。

软路由系统准备

我买的主机就长这样。参数是,x86-64系统,4张网卡,4G内存,32G硬盘。下文中用“软路由”称呼这台机器。

一般这种软路由都是装lede/openwrt之类的系统,不过我要装一个比较正常的ubuntu,原因是后面还要配置自动分流。和店家交流后,他们表示能装ubuntu 16.04,是正常的桌面系统。不过我在一开始折腾的时候把这个系统搞挂了,于是决定换成ubuntu 20.04 server,这样的另一个好处是没有图形系统比较省电,跑桌面系统会肉眼可见地烫手,但是ubuntu server就不会。

这个软路由是带USB+HDMI+VGA接口的,也就是说插上屏幕键盘就能当台电脑用,所以装系统就是很普通的流程:从ubuntu官网上下载amd64版本的.iso文件,用win32diskimager把它刻进U盘,然后把这个安装盘插到软路由的USB接口上,接屏幕键盘启动,启动时按F7进入启动选项,选择由U盘启动,然后跟着流程走即可。这样软路由就装上了ubuntu 20.04 server。

这时需要解决一个奇怪的问题:系统不插VGA监视器无法启动,插HDMI都不行(似乎)。不确定是我软路由硬件的原因还是ubuntu server的原因。这个必须得处理,至少VGA那个分辨率看着也太费劲了……方法是按照这个教程:https://askubuntu.com/questions/825687/what-could-prevent-an-ubuntu-server-from-booting-without-a-vga-connected-monitor,在/etc/default/grub把GRUB_TERMINAL=console一行取消注释,指示系统启动时无需有图形界面。然后用sudo update-grub更新引导选项。在折腾的时候我还在BIOS里面把启动选项设置成了legacy,关闭UEFI,估计应该没用。

至此软路由在系统层面已经准备完毕(按需用apt装包略去不表),接下来开始配置网络。

网络接口配置

在较早的linux系统中,网络接口配置是通过/etc/network/interfaces文件完成的,但在18.04以后改成了一个叫做netplan的东西。netplan的问题是,它没有“allow-hotplug”选项的简单对应,因此还是换回了原来的配置方法:https://askubuntu.com/questions/1031709/ubuntu-18-04-switch-back-to-etc-network-interfaces,主要应该就是安装ifupdown包,然后卸载掉netplan相关的东西:

# systemctl stop systemd-networkd.socket systemd-networkd networkd-dispatcher systemd-networkd-wait-online
# systemctl disable systemd-networkd.socket systemd-networkd networkd-dispatcher systemd-networkd-wait-online
# systemctl mask systemd-networkd.socket systemd-networkd networkd-dispatcher systemd-networkd-wait-online
# apt-get --assume-yes purge nplan netplan.io

然后设置文件/etc/network/interfaces:

source /etc/network/interfaces.d/*

auto lo
iface lo inet loopback

auto enp1s0
allow-hotplug enp1s0
iface enp1s0 inet static
    address 192.168.1.10/24
    gateway 192.168.1.1

iface enp2s0 inet manual
iface enp3s0 inet manual
iface enp4s0 inet manual

auto br0
iface br0 inet static
    bridge_ports enp2s0 enp3s0 enp4s0
        address 192.168.2.1/24

这里面lo就是localhost,不用管。enp1s0~enp4s0是软路由的四个网口,似乎有些机器上,它们的名字会叫作eth0,eth1之类,总之可以用命令ip a查询。我们用enp1s0作为WAN口,其他三个作为LAN口,把这三个桥接起来。

enp1s0需要有allow-hotplug,这个意思是,在软路由启动的时候即使WAN口没插线也能启动。如果不加,它自己倒是没什么问题,只是重启的时候会因为一个什么网络服务没有就绪等待两分钟。但麻烦的是,后面的dhcp服务器会无法正常启动,至少我是这样。可以看到,把WAN口在一级子网中的地址固定为了192.168.1.10,网关是光猫路由一体机。enp2s0~enp4s0的配置都是manual,因为它们的具体地址需要靠网桥br0完成。

auto br0这一句必不可少,它保障了在软路由启动的时候可以自动创建网桥(作为一个虚拟设备,好像是这么个意思),并且恰当地初始化三个LAN口。后面可以看到,bridge_ports这一句把三个LAN口桥接起来,并且规定地址(也就是软路由在二级子网中的地址)是192.168.2.1.

这样配置完/etc/network/interfaces之后,可以sudo service networking restart重启网络服务,然后用ifconfig查看网络信息,应该能看到lo,网桥和四个网口。可以重启软路由验证设置的正确性,重启之后不需要进行任何操作,networking服务会自动读取这个文件里面的配置,并创建相应的接口。

配置sysconf

在/etc/sysctl.conf里面添加如下内容,然后重启机器

net.ipv4.conf.default.rp_filter=0
net.ipv4.conf.all.rp_filter=0
net.ipv4.ip_forward=1
net.ipv6.conf.all.forwarding=1

配置dhcp服务器

首先装包sudo apt install isc-dhcp-server

然后配置/etc/default/isc-dhcp-server,把ipv4和ipv6的接口都设置成我们的虚拟网桥br0:

INTERFACESv4="br0"
INTERFACESv6="br0"

然后设置具体的dhcp规则:

host synology-nas {
    hardware ethernet 00:11:32:93:0d:23;
    fixed-address 192.168.2.9;
}

subnet 192.168.2.0 netmask 255.255.255.0 {
    range 192.168.2.20 192.168.2.100;
    option routers 192.168.2.1;
    option domain-name-servers 192.168.2.1;
}

我这里添加了两条规则。第一条,是群晖NAS的ip地址固定为192.168.2.9,当然后缀2.9也没啥深意,只要不是192.168.2.1(软路由的ip地址)都行。hardware ethernet后面接的是NAS的MAC地址,这个如何查看呢?答案是把NAS用网线连到软路由上之后,用arp -n查看软路由的连接信息,其中的HWaddress项就是MAC地址。

第二条则是子网内的dhcp服务,地址池为192.168.2.20~192.168.2.100,网关是192.168.2.1,而最后一句dns服务器,这里是因为我使用了freedns-go:https://github.com/tuna/freedns-go,所以可以把dns服务器设置成软路由自身。如果没有,那就设成114.114.114.114,1.1.1.1,8.8.8.8之类的。

这样dhcp服务器设置完成,sudo service isc-dhcp-server restart重启dhcp服务器。这时候用网线连到路由器上,用ifconfig就能看到自动分配的地址。

配置自动nat访问互联网

然后还需要nat以访问互联网。命令是:

iptables -t nat -A POSTROUTING -o enp1s0 -j MASQUERADE

当然不能每次重启都敲这么一句,需要把它设置成一个自动启动的服务。ubuntu 20用systemd完成这件事情:https://www.howtogeek.com/687970/how-to-run-a-linux-program-at-startup-with-systemd/

首先我们创建一个文件start_nat.sh,内容如下:

#!/bin/sh
iptables -t nat -A POSTROUTING -o enp1s0 -j MASQUERADE

并且用sudo chmod +x start_nat.sh赋予执行权限,然后把它放到/usr/local/bin目录下。然后再创建一个文件/etc/systemd/system/start_nat.service,内容为

[Unit]
Description=Start NAT

[Service]
ExecStart=/usr/local/bin/start_nat.sh
RemainAfterExit=true

[Install]
WantedBy=multi-user.target

Description随便填。如果需要在.sh里面启动tmux(并且有可能在tmux里面运行某个程序),这个RemainAfterExit就必须要有。ExecStart是这个服务在启动的时候执行什么命令,也就是前面那个脚本start_nat.sh。然后我们把它的权限改成640:

sudo chmod 640 /etc/systemd/system/start_nat.service

这样我们就创建了一个名叫start_nat的服务。然后把它注册为启动时运行:

sudo systemctl daemon-reload
sudo systemctl enable start_nat
sudo systemctl start start_nat

最后那个start是现在就启动的意思。你也可以重启一次软路由。如果一切配置恰当,重启之后,连接软路由的电脑应该能直接连接到一个LAN并且正常访问互联网,无需进行任何操作。

NAS端口转发配置

到现在为止我们拥有了一个正常功能的路由器。但是我还希望在整个家庭局域网内访问NAS的各项服务,所以还需要额外折腾一下。具体来说有四个:web控制页面5051端口,webdav远程硬盘5006端口,synology drive群晖本地云6690端口,smb共享文件(比如说在‘我的电脑’里敲个\\192.168.2.9访问共享文件夹)445端口。

因此我们需要在软路由上配置四个端口转发,前三个随便更改了一下端口名称:

  • 软路由的16000端口转发至NAS的5051端口
  • 软路由的16006端口转发至NAS的5006端口
  • 软路由的16090端口转发至NAS的6690端口
  • 软路由的445端口转发至NAS的445端口

注意之前把NAS的ip地址固定为192.168.2.9了。每一个转发可以通过两条iptables命令实现,因此总共是八条,把它们添加到之前/usr/local/bin/start_nat.sh的末尾,这样每次软路由一启动就会自动完成转发设置。总的下来,start_nat.sh如下:

#!/bin/sh
iptables -t nat -A POSTROUTING -o enp1s0 -j MASQUERADE
iptables -t nat -A PREROUTING -p tcp --dport 16000 -j DNAT --to-destination 192.168.2.9:5051
iptables -t nat -A POSTROUTING -p tcp -d 192.168.2.9 --dport 5051 -j MASQUERADE
iptables -t nat -A PREROUTING -p tcp --dport 16006 -j DNAT --to-destination 192.168.2.9:5006
iptables -t nat -A POSTROUTING -p tcp -d 192.168.2.9 --dport 5006 -j MASQUERADE
iptables -t nat -A PREROUTING -p tcp --dport 16090 -j DNAT --to-destination 192.168.2.9:6690
iptables -t nat -A POSTROUTING -p tcp -d 192.168.2.9 --dport 6690 -j MASQUERADE
iptables -t nat -A PREROUTING -p tcp --dport 445 -j DNAT --to-destination 192.168.2.9:445
iptables -t nat -A POSTROUTING -p tcp -d 192.168.2.9 --dport 445 -j MASQUERADE

这样,比如说,在家庭局域网任意位置访问192.168.1.10:16000即可访问NAS的5051端口。

路由器桥接,收尾工作

在折腾之前家里的路由器不是桥接模式而是路由模式。需要更改一下。总的来说,从路由模式切换成桥接模式需要做三件事:

  • 把WAN口网线拔掉插到LAN口
  • 登录路由器管理界面,关闭dhcp
  • 在路由器管理界面手动设置一个LAN口ip地址,不要冲突

其中我遇到了一件奇怪的事情,就是设置TP-LINK路由器的时候,一开始忘了关dhcp,后来才关掉,但此时连接它WIFI的电脑会始终显示网关是192.168.2.2,也就是这台路由器的地址,而非正确的软路由地址192.168.2.1。然后通过在电脑上手动设置ip,dns服务器和网关(设为软路由地址),然后再把ip改回自动模式,解决了问题,这个方法有点奇怪但总之work了,原理不明。

最后结果如图:

Real-time exploration of regular volume data by adaptive reconstruction of isosurfaces读书笔记

本文讨论使用Octree(八叉树)重建表面的问题。

3 Continuous isosurfaces

讨论连续性问题。

如果在Octree上储存数据,那么在不同大小的格子交界的地方,数据可能会不连续。

本文采用的方法:在大小格子交界处,让大格子在小格子上sample数值:

图1

也就是用方块节点插值出来三角节点的值,把原来的扔掉。这样一来,小格子上的数值就是线性的,那么marching cubes就会在大小格子上得到同样的intersection point,也就是说数值是连续的。

但是,虽然数值连续,但算出来的表面不一定连续。如图3:

图3

本文的方法,首先算法保证相邻的叶子格子不会相差超过一层(one generation,应该是这个意思)。也就是都是图3这种情况:一个大格子和四个小格子相邻。然后如图6所示:在marching cubes用大格子算出来的三角形重心处加一个点,把这个三角形切成一堆“扇叶”,把小格子算出来的那些交点加到这堆“扇叶”的外圈,切得更细就行了。

图6

关于马拉戈尼效应(Marangoni Effect)

本文基于维基百科内容写作:

https://en.wikipedia.org/wiki/Marangoni_effect

https://en.wikipedia.org/wiki/Marangoni_number

示例视频(需科学上网):

https://upload.wikimedia.org/wikipedia/commons/3/38/Marangoni_effect_experimental_demonstration.ogv

马拉戈尼效应由表面张力的梯度产生。在简化情况下,流体被马拉戈尼效应驱动而流动的特征速度[latex]u\sim\Delta\gamma/\mu[/latex],其中[latex]\delta\gamma[/latex]是表面张力的特征差距(difference in surface tension),[latex]\mu[/latex]就是流体的粘度。至于为什么,见后面马拉戈尼数的部分。

这里注意,Pa=N/m^2,所以[latex]Pa\cdot s[/latex]的单位是[latex]N\cdot s/m^2[/latex].所以[latex]\Delta\gamma/\mu[/latex]的单位就是速度单位N/m。

在室温下,水的表面张力高达约0.07N/m,而粘度则约为[latex]10^{-3}Pa\cdot s[/latex].由于表面张力和粘度之间的巨大差距,只要水的表面张力能稍微变化一点(例如,从0.070到0.069),这个特征速度就能达到1m/s这个量级。

然后我们看马拉戈尼数(Marangoni number).它的记号是Ma,注意不要跟马赫数搞混了。

马拉戈尼数是一种Peclet数,定义如下:

[latex]Ma=\frac{马拉戈尼效应产生的对流输运速率}{表面张力源产生的扩散输运速率}[/latex]

这个扩散输运,我的理解,比如说如果表面张力之差是由表面活性剂(surfactant,如肥皂)的浓度梯度产生,那么它指的就是肥皂分子的扩散速率。

举这样一个例子:假设有一水平面,表面张力在L的距离上变化了[latex]\delta\gamma[/latex]。假设水足够深,深度至少为L,那么L就是问题中的唯一特征尺度。我们可以把这个问题叫做马拉戈尼流(Marangoni flow).

这样,就可以用Stokes flow计算这个问题。也就是:令流体柯西张量的梯度等于粘性扩散。表面张力产生的张量,尺度为[latex]\Delta\gamma/L[/latex].而粘性张量的尺度为[latex]\mu u/L[/latex],u就是特征速度。让这二者相等,得到流速[latex]u=\Delta\gamma/u[/latex].

这样的话,马拉戈尼数里面,对流输运的特征速率就是特征速度乘以特征尺度,即uL.而扩散输运的特征速率那就是菲克定律(Fick’s law)里面的扩散系数D。所以总的马拉戈尼数就是:

[latex]Ma=\frac{uL}{D}=\frac{\Delta\gamma L}{\mu D}[/latex].

我们可以用量纲分析法小小地验算一下这个式子。

表面张力[latex]\gamma[/latex]的单位是N/m,其中N的量纲是[latex]MLT^{-2}[/latex],那么表面张力之差[latex]\Delta\gamma[/latex]的量纲,和表面张力一样,就是[latex]MT^{-2}[/latex],马拉戈尼数分子的量纲就还是[latex]MLT^{-2}[/latex].

另一边,扩散系数(或称质量扩散率)D的单位是[latex]m^2/s[/latex],量纲是[latex]L^2T^{-1}[/latex].而前面说过,粘度的单位是[latex]N\cdot s/m^2[/latex],量纲是[latex]ML^{-1}T^{-1}[/latex].用扩散系数的量纲乘以粘度的量纲,就得到分母的量纲[latex]MLT^{-2}[/latex],分子分母恰好抵消。

OpenGL屏幕像素坐标问题

传给glutMouseFunc()的(x,y)是左上至右下的坐标:x表示离屏幕左边缘的像素数,y表示离屏幕上边缘的像素数,如图:

而glReadPixels()采用的坐标系统则与此不同,它是从左下到右上的。也就是说,x表示离屏幕左边缘的像素数,y表示离屏幕下边缘的像素数,如图:

所以,如果想把glutMouseFunc()取得的坐标传给glReadPixels(),就需要做一个变换:

y=Height-y-1