博文
Hamaker方程
||
若粒子之间的van der Waals吸引势与距离的6次方成反比, Hamaker方程给出了两个宏观球形粒子之间的吸引能. 当粒子之间的吸引势取其他形式时, 可利用下式进行计算$:$
E=−π2q2λz∫z+R2z−R2[R22−(z−R)2]dR∫R+R1R−R1R21−(R−r)2rn−1dr
令 α=R1,β=R2,κ=−π2q2λz
E=κ∫z+βz−β[β2−(z−R)2]dR∫R+αR−αα2−(R−r)2rn−1dr=κ∫z+βz−β[β2−(z−R)2]In(R)dR
In=∫R+αR−αα2−(R−r)2rn−1dr=R2−α2n−21rn−2−2Rn−31rn−3+1n−41rn−4(n≥4)
当 n≤4 时, 此积分需要单独考虑.
n=1I1E1=43α3=κ43α343β3=−Q1Q2λz
这种情形下, 两个球形粒子之间的势能与将其视为处于中心的质点是一样的. 引力, 静电相互作用就属于这种情况, 所有我们可以将星球与荷电小球视为质点, 利用万有引力定律或是库伦定律计算势能.
n=2I2E2/κ=2αR+(R2−α2)lnR−αR+α=∫z+βz−β2αR[β2−(z−R)2]dR+∫z+βz−β(R2−α2)[β2−(z−R)2]lnR−αR+αdR=83αβ3z+∫z+βz−β[−R4+2zR3+(α2+β2−z2)R2−2zα2R+α2(z2−β2)]lnR−αR+αdR
相关的积分公式如下
∫R4lnR−aR+adR∫R3lnR−aR+adR∫R2lnR−aR+adR∫R1lnR−aR+adR∫lnR−aR+adR=110(−2a5ln(R−a)−a(2a4ln(a+R)+2a2R2+R4)+2R5lnR−aR+a)=110[−aR4−2a3R2−2a5ln(R2−a2)+2R5lnR−aR+a]=112(−3a4ln(R−a)+3a4ln(a+R)−6a3R+3R4lnR−aR+a−2aR3)=112[−2aR3−6a3R+3(R4−a4)lnR−aR+a]=13(−a3ln(R−a)−a(a2ln(R+a)+R2)+R3lnR−aR+a=13[−aR2−a3ln(R2−a2)+R3lnR−aR+a]=12(−a2ln(R−a)+R2lnR−aR+a+a(aln(a+R)−2R))=12[−2aR+(R2−a2)lnR−aR+a]=−aln(R−a)+RlnR−aR+a−aln(a+R)=−aln(R2−a2)+RlnR−aR+a
=−−+∫(A4R4+A3R3+A2R2+A1R+A0)lnR−aR+adRA410aR4−A36aR3−(A45a2+A23)aR2−(A32a2+A1)aR(A45a4+A23a2+A0)aln(R2−a2)[A45R5+A23R3+A0R+A34(R4−a4)+A12(R2−a2)]lnR−aR+a
下面的eigenmath代码用于计算定积分
# Language: eigenmathA4 =-1A3 = 2zA2 = a^2+b^2-z^2A1 =-2 z a^2A0 = a^2(z^2-b^2)I4(R)=-A4 a R^4/10I3(R)=-A3 a R^3/6I2(R)=-(A4 a^2/5+ A2/3) a R^2I1(R)=-(A3 a^2/2+ A1) a RI00(R)=-(A4 a^4/5+ A2 a^2/3+ A0) a # ln(R^2-a^2)I01(R)=(A4 R^5/5+ A2 R^3/3+ A0 R + A3 (R^4-a^4)/4+ A1 (R^2-a^2)/2)# ln((R-a)/(R+a))I41(R)= I5+I4+I3+I2I41(z+b)-I41(z-b)simplifyI00(z+b)simplifyI00(z-b)simplifyI01(z+b)simplifyI01(z-b)simplify最终结果
E2/κ=83αβ3z+215αβz[z2+11α2−9β2]+215α3[α2+5(z2−β2)]ln(z−β)2−α2(z+β)2−α2+215β3[β2+5(z2−α2)]ln(z−α)2−β2(z+α)2−β2+z30[−z4+10(α2+β2)z2+15(α2−β2)2]lnz2−(α−β)2z2−(α+β)2=215αβz[z2+11(α2+β2)]+215α3[α2+5(z2−β2)]ln(z−β)2−α2(z+β)2−α2+215β3[β2+5(z2−α2)]ln(z−α)2−β2(z+α)2−β2+z30[−z4+10(α2+β2)z2+15(α2−β2)2]lnz2−(α−β)2z2−(α+β)2
n=3I3E2/κ=2RlnR+αR−α−4α=∫z+βz−β−4α[β2−(z−R)2]dR+∫z+βz−β2R[β2−(z−R)2]lnR+αR−αdR=−163αβ3+∫z+βz−β[2R3−4zR2+2(z2−β2)R]lnR−αR+adR=−23αβ[z2+3(α2+β2)]−43α3zln(z−β)2−α2(z+β)2−α2−43β3zln(z−α)2−β2(z+α)2−β2+16[z4−6(α2+β2)z2−3(α2−β2)2]lnz2−(α−β)2z2−(α+β)2
n=4I4E2/κ=2αRR2−α2+lnR−αR+α=∫z+βz−β2αRR2−α2[β2−(z−R)2]dR+∫z+βz−β[β2−(z−R)2]lnR+αR−αdR=∫z+βz−β2αβ2RR2−α2dR−∫z+βz−β2αR(z−R)2R2−α2dR+∫z+βz−β[β2−(z−R)2]lnR+αR−αdR=4αβz+α(β2−α2−z2)ln(z+β)2−α2(z−β)2−α2−2α2zlnz2−(α+β)2z2−(α−β)2−83αβz+α(z2+α23−β2)ln(z+β)2−α2(z−β)2−α2+z(−z23+β2−α2)lnz2−(α−β)2z2−(α+β)2+23β2ln(z−α)2−β2(z+α)2−β2=43αβz+23α3ln(z−β)2−α2(z+β)2−α2+23β3ln(z−α)2−β2(z+α)2−β2+z(13z2+α2+β2)lnz2−(α−β)2z2−(α+β)2
n=5E2/κ=23αβ+z2−α2−β23lnz2−(α−β)2z2−(α+β)2
n=6E2/κ=z6lnz2−(α−β)2z2−(α+β)2−2αβz3α2+β2−z2[α2−(z−β)2][α2−(z+β)2]=z6lnz2−(α−β)2z2−(α+β)2+αβz3[1z2−(α+β)2+1z2−(α−β)2]
下面的 matlab 代码可以帮助推导 n>4 时的结果
# Language: matlabclc;syms n z r R a b In I f E real% In=(R^2-a^2)*r^(2-n)/(n-2)-2*R*r^(3-n)/(n-3)+ r^(4-n)/(n-4);%% m=6;%% I=subs(In, n, m);% I=In;% I=subs(I, r, R+a)-subs(I, r, R-a);% I=simple(I);% I%pretty(I)I=2*(a*n+exp(2*a*n)*(a*n-1)+1)*exp(-n*(a+R))/n^3f=(b^2-(z-R)^2)*Ipretty(f)I=int(f,R);I=simple(I);% I=subs(I,'sqrt(-1)','j');% I=subs(I,'atan((R*j)/a)','j*log((a+R)/(a-R))/2');% I=collect(I,'j');% I=subs(I,'j^2','-1');Ipretty(I)E=subs(I, R, z+b)-subs(I, R, z-b);E=simple(E);% E=collect(E,'a^2+b^2-z^2')Epretty(E)可见, 即便对于均匀分布的球形粒子, 其作用力的一般形式也很复杂, 但综合看起来, 表达式中都含有对数形式, 这可能预示着在构造势函数时要考虑对数形式的函数, 但常见的势能函数中却很少使用这种形式.
推广一下, 当粒子之间的相互作用为指数势 e−nr 时, 得到的结果更简单,
IE2/κ=∫R+αR−α[α2−(R−r)2]re−nrdr=2e−n(α+R)n(α(n(α+R)+3)+R)−e2nα((n(α(nα−nR−3)+R)+3)+3)n4=Ae−nR(BR−C)=∫z+βz−β[β2−(z−R)2]IdR=e−n(z+β)[23βB−C+zB+(C+3βB−zB)e2nβn3−2βC−βB−zB+(C+βB−zB)e2nβn2−6Be2nβ−1n4]
其中
A=2n4e−nα,B=n[1+nα+(nα−1)e2nα],C=n2α2+3nα−e2nα(n2α2−3nα+3n+3)
图片表格公式代码完整版请参看Hamaker方程https://blog.sciencenet.cn/blog-548663-845023.html
上一篇：GROMACS中文手册翻译
下一篇：数据展示：请选择好的颜色映射方案
