液滴高度求解模型
求解模型液滴高度液滴高度求解模型
摘要
对液滴外形轮廓进行研究从而得到液滴的表面张力、接触角，体积以及液滴高度等物理参数，一直是流体研究的重要内容。本文针对水平固体表面上液滴高度随体积的变化规律建立数学模型。通过拉普拉斯方程对液滴极限高度进行了求解，并计算了液滴接触角变化时饱和高度与极限高度的比、液滴直径及饱和体积。通过建立饱和液滴轮廓图的微分方程模型，画出了它的正视轮廓图。
对于问题一，分析液滴的高度变化过程，建立拉普拉斯方程模型。在液滴体积逐渐变大时，液滴高度增加到饱和高度此时的拉普拉斯方程中的曲率半径相等，利用流体的静态平衡Laplace 公式求解出了饱和高度h max =
2γls (1−cos θ)
。
ρg
液滴的体积继续增大，高度下降。达到极限高度时拉普拉斯方程的曲率半径R 2无穷大，在利用流体的静态平衡公式求解出了极限高度h min =
γls (1−cos θ)
。
ρg
对于第二问，在第一问的模型基础上我们可以求出对应的饱和高度与极限高度的比、液滴直径及饱和体积的表达式，再由Matlab 编程从0度到180度每隔10度计算得到结果见表1。并对求出的结果绘制了随接触角θ变化的关系图。分析其变化趋势，补充说明了模型的计算结果。
对于第三问，建立饱和液滴轮廓图的微分方程模型，依据附加压强的
Laplace 公式对液滴的正视图轮廓线上任意一点p 进行了受力分析，得到各变量
之间的关系式。然后，为了便于确定微分方程组的初值条件，以弧长s （s 为顶点到p 点的曲线弧长）为自变量，建立轮廓线上的任意一点p 的纵坐标z ，横坐标x 以及与x 轴的夹角ϕ关于s 的微分方程组模型。运用龙格-库塔法编程求解，绘制出了饱和液滴正视轮廓图，并附上图形各参数。结果见图8、表2。
关键词：液滴高度，接触角，饱和高度，极限高度，Lapace 方程，龙格-库塔法
一、问题重述
在物理实验中发现一个有趣的现象如下：
测量放在一种固体材料的水平平面上具有不同体积的液滴在静态时的高度时，发现当该液滴与固体的接触角θ>0不变的情况下，随着液滴体积的递增，液滴的高度递增，直到液滴体积达到某个(在此称之为)饱和体积时，液滴高度达到最大值(在此称之为饱和高度).当液滴体积从饱和体积开始递增时，液滴的高度递减,而且随着体积的增大高度递减量越来越小，液滴高度似乎趋于一(在此称之为)极限高度：
1.请对于一般的接触角θ∈(0,π](弧度)建立数学模型解释以上现象，并给出极限高度的表达式.
2.对于您所建立的数学模型，请对于接触角θ从10到180每间隔10度°计算出对应的饱和高度与极限高度的比、液滴直径(单位:毛细长度)及饱和体积(单位:毛细体积)，最终计算的结果均舍入到四位有效数字(用计算机的科学表示法，如123.4写成1.234e2)要求与精确值的绝对误差在最后一位的半个单位之内.3.画接触角为180度、液滴体积等于饱和体积时液滴的正视轮廓图.
二、问题分析
对液滴外形轮廓进行研究从而得到液滴表面张力，接触角，体积等参数一直是流体力学研究的重要内容。对于水平固体平面上的液滴，其形态主要受表面张力与重力作用，表面张力使液滴有聚合趋势，重力作用使其向周围扩散，在两种力的综合作用下，便形成一个类似于凸平透镜的液垫。由文献[1]知，液滴处在固体表面上时有两种平衡形态。对于水平固体平面上的液滴，其形态主要受表面张力与重力作用，表面张力使液滴有聚合趋势，重力作用使其向周围扩散，因此，当表面重力占主导地位时，会使液面形成一个类似于图1的液垫，而当表面张力占主导地位时，由于液滴在张力作用下总是力图使自己保持最小的表面积，结合几何知识，可将液滴模型表面简化为如图2的球冠状。
针对问题一，根据题意，我们将高度的变化分为两个阶段，第一个阶段为液滴达到饱和之前，考虑到液滴在固体表面上不同平衡形态，我们又将第一阶段分为两种情况分别进行求解，对于第一种情况（如图1) ，将其视为球缺，利用求解表面张力的lapace 方程，建立模型进行求解。对于第二种情况(如图2) ，根据势能关系、铺展系数等建立模型进行求解。在第二阶段液体达到饱和体积后，我们
根据拉普拉斯定律，热力学和流体力学相关知识建立微分方程模型来求液滴
的极限高度。
图1液滴的球冠模型图2液滴的液垫模型
针对问题二，要求给出接触角θ从10到180每间隔10度对应的饱和高度与极限高度的比、液滴直径及饱和体积，分析可知，要求液滴直径，故我们可将液滴等效为问题一中的球冠模型进行求解。
针对问题三，对液滴的截面图线上任意一点P 做受力分析，根据附加压强的拉普拉斯公式，以及适当选取正交截面，建立平面直角坐标系如图3所示，得到曲面上每一点的纵坐标，横坐标，去路半径与z 轴的夹角关于自变量s （s 为接触点到P 点的曲线弧长）的微分方程组。根据得到的微分方程组，可以运用龙格-库塔法，得出其数值解，是一些离散的点。这时，一组（x ，z ），便可唯一确定一个液滴轮廓线，其中，x 表示液滴轮廓线的横坐标，z 表示液滴轮廓线的纵坐标。
三、问题假设
1、固体材料表面为理想表面。2、固体材料表面水平。
3、固体材料水平面是各向均匀的。
4、液体的性质（比如密度分布）是均匀的。5、温度始终保持室温20摄氏度。6、大气压分布是均匀的。
7、液体体积增加的过程为准静态过程，液体的高度取平衡状态下的高度。8、液体不会被固体材料吸收或发生化学反应。
四、符号说明
符号
符号含义
液滴接触面的曲率半径液体内的压力、液体外的压力液滴表面张力系数
液滴与接触面的接触角液滴饱和高度，极限高度液体密度重力加速度液滴直径
饱和液滴体积毛细长度
液滴轮廓线在p 点处的切线与x 轴夹角
R 1R 2R p 1p 2γθh max ，h min
ρg
D V
h ∗ϕ
五、模型的建立
5.1问题一模型的建立
对于液体表面张力的研究，其满足拉普拉斯方程[2]。
拉普拉斯方程表示液面曲率与液体压力之间的关系的公式。一个弯曲的表面称为曲面，通常用相应的两个曲率半径来描述曲面，即在曲面上某点作垂直于表面的直线，再通过此线作一平面，此平面与曲面的截线为曲线，在该点与曲线相重合的圆半径称为该曲线的曲率半径R 1。通过表面垂线并垂直于第一个平面再作第二个平面并与曲面相交，可得到第二条截线和它的曲率半径R 2，用R 1与R 2可表示出液体表面的弯曲情况。若液面是弯曲的，液体内部的压力p 1与液体外的压力p 2就会不同，在液面两边就会产生压力差∆p =p 1−p 2, 其数值与液面曲率
大小有关，可表示为：∆p =γ(拉斯方程。
11
+) 式中γ是液体表面张力。该公式成为拉普R 1R 2
我们知道当该液滴与固体的接触角θ>0不变的情况下，随着液滴体积的递增，液滴的高度递增，直到液滴体积达到某个(在此称之为)饱和体积时，液滴高度达到最大值(在此称之为饱和高度).当液滴体积从饱和体积开始递增时，液滴的高度递减,而且随着体积的增大高度递减量越来越小，液滴高度似乎趋于一(在
此称之为)极限高度，此过程可用图1描述。
液滴体积增大，高度增加到饱和高度h max
液滴体积继续增大，高度减小到，直到极限高度h min
图3液滴高度变化过程图
由普遍形式的拉普拉斯定律，液体受表面张力作用的压力p 可表示为：
∆p =γ(
11
+) R 1R 2
（1）
当液滴由小变大时，高度h 增加，达到饱和高度h max ，此二力的平衡导致了对于球状液滴的拉普拉斯定律：
p =2γls /R
由图像中的几何图1关系可得：
（2）
R =
h max 1−cos θ
（3
）
图4液滴曲率半径几何图
联立（2）（3）式得到
p =
由流体静力学压强公式得到：
2γls (1−cos θ)
h max
（4）
p =ρgh max
整理（4）和（5）式可得饱和高度：
（5）
h max =
2ls (1−cos )
ρg
（7）
液滴的体积继续增大，高度由饱和高度减小，直到达到极限高度h min ，此时液层面积充分大，R 2将充分大，拉普拉斯方程可表示为：
p =
由几何学知识得：
γls R h min 1−cos θ
（8）
R =
将（9）代人（8）式可得：
（9）
p =
由流体静力学压强公式得到：
γls (1−cos θ)
h min
（10）
p =ρgh min
将上式整理得出极限高度表达式：
（11）
h min =
由球缺的面积积分可知：
γls (1−cos θ)
ρg
（12）
dV =πx 2dy =π(R 2−y 2) dy
V 缺=∫dV =∫
R −h max
R R
R −h max
πx 2dy =∫
13
(R 2−y 2）dy =πRh 2−πh max max
R −h max
3
R
（13）
5.2问题二模型的计算
首先对于不同的液体，他们的表面张力不同，所以求解结果不同。在这里我们去自然界最常见的液体水，在20℃状态下的求解。
根据资料可知水在20℃时的表面张力为γls =72.8×10-3，密度
ρ=0. 998×103kg /m 3，取重力加速度g =9. 80665m /s 2。
根据模型1的结果由（7）、（12）式可知
h min =
γls (1−cos θ)
ρg
h max =2γls (1−cos θ) /ρg
故其高度比可以表示为
h max h min
由（9）式可知液滴半径为R =
h min
，那么直径
1−cos θ
2h min
D =2R =
1−cos θ
由（13）式可知饱和体积为
13
V =πRh 2−πh max max
3
根据以上公式对接触角θ从10°到180°每间隔10度计算,编程计算（编程见附件1）得到下列结果表1：
表1：计算结果
θ
10°20°30°40°50°60°70°80°90°
高度比1.4141.4141.4141.4141.4141.4141.4141.4141.414
直径6.26e-23.14e-22.1e-21.60e-21.30e-21.10e-20.96e-20.84e-20.78e-2
饱和体积2.2106e-84.3379e-86.3035e-88.0393e-89.4911e-810.6223e-811.4153e-811.8731e-812.0177e-8
θ
100°110°120°130°140°150°160°170°180°
高度
直径
比
1.4140.72e-21.4140.66e-21.4140.62e-21.4140.60e-21.4140.58e-21.4140.56e-21.4140.56e-21.4140.54e-21.4140.54e-2
饱和体积11.8890e-811.5412e-811.0390e-810.4527e-89.8535e-89.3079e-88.8734e-88.5941e-88.4978e-8
对以上结果我们做一定分析。
由Matlab
编程绘制出液滴接触角变化时，饱和高度，极限高度的变化关系图
图5接触角变化时饱和高度极限高度变化图
从图3我们可以看出接触角由0度到180度增大时，液滴的饱和高度、极限高度都在增大，并且这种变化趋势在减小。变化关系是一致的，这也符合高度比为1.414
。
图6液滴半径和液滴饱和体积变化图
从图4我们可以看出对于接触角的增大，液滴的半径是在减小的，液滴的体积是先增大，后减小。并且在接触角为90度时液滴的饱和体积最大。
5.3问题三5.3.1模型建立
由文献[3]可知对于给定的液滴与固体表面，当液滴稳定存在于固体表面之上的时候，接触角θ即可确定。根据液面的拉普拉斯方程可知，在不考虑重力时，弯曲液面上任意一点p 在液面两边产生的压力差即附加压强差∆p 为：
⎛11⎞⎜⎟∆p =γls +
⎜R p 1R p 2⎟⎝⎠
（14）
其中∆p 为接触面的压力差，根据流体力学知识∆p 为大于0的常数，γls 为液面与气体接触面的表面张力，R p 1，R p 2为p 点的第一和第二曲率半径。
对于弯曲液面上任意一点p 和x 0，他们的附加压强差等于内部重力长生的压强差，即：
⎛11
γls ⎜+⎜R x 1R x 2
0⎝0
⎞⎛⎞
⎟=γls ⎜1+1⎟−ρgz 1
⎜R ⎟⎟
⎝p 1R p 2⎠⎠
（15）
其中R x 01，R x 02分别为参考点（本问题中为液面的顶点）的第一和第二曲率半径，z 1为p 点距离参考平面的垂直距离。
又根据题目对模型的要求，以毛细长度为单位进行计量，对公式（15）进行
无量纲化处理，单位长度为毛细单位：
h *=
对（15）进行变换得：
γls ρg
（16）
⎞⎛⎞γls ⎛⎜1+1⎟=h *2⎜1+1⎟−z 1
⎜R ⎟⎟ρg ⎜⎝p 1R p 2⎠⎝R x 01R x 02⎠
（17）
⎛11⎞⎜⎟，ρ，g 均为常量，故上式左端为常数，而右端h *2为单又由于γls +
⎜R x 1R x 2⎟
0⎝0⎠位长度的平方。为便于计算，我们用k 代替上式左端，化简得
k =
由微分几何知识得：
11
+−z 1R p 1R p 2
（18）
1−x ' ' (z ) =
R 11+x ' 2(z )(（19）
11=
R 2x (z )1+x ' 2(z )
3
（20）
对于液滴的体积，利用微分法得：
dV =πx 2
dz
（21）
图7液滴轮廓分析图
基于以上的推导，我们参照图7所示建立的坐标系，建立如下的微分方程模型，首先，为了便于计算求解，我们将x ，z 分别表示成关于s 的函数，即
x 1=x 1(s 1) z 1=z 1(s 1)
其中s 1为从原点到该点的弧线长度。由图5可知
（22）
dx 1=ds 1cos ϕdz 1=ds 1sin ϕ
对（21）式变换得：
（23）（24）
dV =πx 2sin ϕds
由（18）式得：
11=k −+z R p 1R p 2
由式（20）得
11
=
R p 2x (z )1+x ' 2(z )
（25）
（26）
=1
⎛dx ds ⎞x [1+⎜⎟]ds dz ⎝⎠
1
2=
=
将1代入式（26）得：R p 2x [1+(cosϕ) 2]sin ϕx 1sin ϕ=k −+z R p 1x
又对弧长s 有：（27）
ds =R p 1d ϕ即1d ϕ=R p 1ds
（28）故1sin ϕd ϕ=k −+z =R p 1x ds
由（9），（10），（11），（14）可得如下饱和液滴轮廓图的微分方程模型
⎧dx ⎪ds =cos ϕ
⎪dz ⎪=sin ϕ⎪ds ⎨dV ⎪=πx 2sin ϕ⎪ds sin ϕ⎪d ϕ=k +z −⎪x ⎩ds
5.3.2模型轮廓图的绘制
根据得到的微分方程模型，可以运用龙格-库塔法编程求解（程序见附录3），得出其数值解，并根据这些解绘制液滴的正视轮廓图如图5。
图8接触角为180度时饱和液滴的正视轮廓图表2饱和液滴正视图的各项参数
常数K 值
0.0159液滴饱和高度2.1066液滴与固体接触面直径10.3654液滴体积1.46e02
六、模型的评价
针对本文所建立的模型，我们在模型建立中没有将液滴高度与体积变化过程完整的展现出来，只是根据临界条件计算出液滴的饱和高度与极限高度，进而得出液滴的体积。这样就不能直接解释题目中所说明的现象，只能根据所得的结果进行间接的推导。而且在用数值解法画液滴轮廓图中，用等步长搜索法求解合适的常数K，使得用MATLAB 编程进行求解时，计算效率极低，方法还需要改进。
七、参考文献
[1]陈逸云. 液滴高度的变化. 科技创新导报。2013年
[2]任素贞. 王旭珍. 物理化学（第三版）. 上海：上海科学技术出版社.2012年12月
[3]宁乔. 图像法求液滴表面张力和接触角. 空间科学. 2008年
[4]RotenbergY，BoruvkaL ，NeumannA W ．Determination of surface tension and contact angle from the shapes of axisymmetric fluid interface ．J ．Colloid
InteSci ．1983，93(1)：169~182
八、附录
附录1：
问题二的数值计算程序
clear clc
format long
for i=1:1:18
theta=i*pi/18;
hmin(i)=sqrt((72.8e-3*(1-cos(theta)))/(0.998e3*9.80665));
hmax(i)=sqrt(2*(72.8e-3*(1-cos(theta)))/(0.998e3*9.80665));
R(i)=hmax(i)/(1-cos(theta));
V(i)=pi*R(i)*hmax(i)^2-(1/3)*pi*hmax(i)^3;
end
hmax,hmin,bizhi=hmax./hmin,R,V
附录2：
图5，图6的绘制程序
clear clc
format long
for i=1:1:18
theta=i*pi/18;
hmin(i)=sqrt((72.8e-3*(1-cos(theta)))/(0.998e3*9.80665));
hmax(i)=sqrt(2*(72.8e-3*(1-cos(theta)))/(0.998e3*9.80665));
R(i)=hmax(i)/(1-cos(theta));
tiji(i)=(pi*R(i)*hmax(i)^2-(1/3)*pi*hmax(i)^3);
end
hmax,hmin,bizhi=hmax./hmin,R,tiji
theta2=10:10:180
figure(1)
plot(theta2,hmax,'ro',theta2,hmin,'b*',theta2,hmax,'r',theta2,hmin,'b')legend('液滴饱和高度',' 液滴极限高度')
xlabel('接触角度')
ylabel('液滴高度')
title('液滴高度随接触角变化')
figure(2)
subplot(2,1,1)
plot(theta2,R,'r',theta2,R,'ro')
xlabel('接触角度')
ylabel('液滴半径')
title('液滴半径随接触角变化')
subplot(2,1,2)
plot(theta2,tiji,'b',theta2,tiji,'b*')
xlabel('接触角度')
ylabel('液滴体积')
title('液滴体积随接触角变化')
附录3
问题三中轮廓图绘制程序
1. 等步长计算K 值
sl=100;
theta=pi;
a=0.002;
b=0.2;
k=[a:(b-a)/sl:b];
for i=1:(sl+1)
[h(i),V(i),D(i),s(i)]=solveeq(k(i),theta);
end
c=find(h==max(h));
hm=h(c)
K=k(c)
S=s(c)
Hb=hm/(2*sin(theta/2))
D=D(c)
v=V(c)
2. 等步长计算K 值附属函数
function [h,v,D,s]=solveeq(k,theta)
Options=odeset('Abstol',1e-8,'Events',@eventsn);
m=[theta,k];
[TOUT,YOUT,TE,YE,IE]=ode45(@odefun12,[06],[0000],Options,m);upb=6.1;
while numel(YE)==0
[TOUT,YOUT,TE,YE,IE]=ode45(@odefun12,[0upb],[0000],Options,m);upb=upb+0.1;
end
h=YE(2);
v=YE(4);
s=TE;
D=2*max(YOUT(:,1));
End
3. 控制ode43函数终止条件附属函数
function [value,isterminal,direction]=eventsn(s,y,m)
value=y(3)-m(1);
isterminal=1;
direction=1;
4. 根据计算出的K 值画轮廓图
options=odeset('Abstol',1e-8,'Events',@eventsn);
m=[pi0.0159];
[TOUTYOUT TE YE IE]=ode45(@odefun12,[010],[0000],options,m);x=YOUT(:,1);
z=YOUT(:,2);
z=max(z)-z;
plot([-66],[max(z)max(z)],'k--')
hold on
plot(x,z,'r')
hold on
plot([-1010],[00],'k')
hold on
plot(0,max(z),'o')
x2=-x;
plot(x2,z,'r')
axis([-1010-26])
text(0,max(z)+0.5,num2str(max(z)));
title('接触角为180度时液滴轮廓图')
xlabel('横坐标X')
ylabel('纵坐标Z')
hold off