| 查看: 7184 | 回复: 7 | ||
[求助]
lammps 计算粘度 lj单位 转化为 metal单位 出错 已有3人参与
|
|
最近做粘度的计算，修改lammps里面提供的例子，Muller-Plane的方法，但是计算的结果对不上，所以一步步排查，最后发现单位转换这里都没搞定。 原例子都是用lj单位写的，我想把它改成metal单位，就是按照manual里面lj单位的说明改，改了的地方都下划线标出了，想请大家帮忙看看是不是哪里错了。确实很繁琐，我把过程都写出来了，还请大家多多帮忙，也好帮助其他遇到这个问题的人。最后附上全文的word版，红色标出比较明显，还有lj单位和metal单位的in文件和log文件。  http://tieba.baidu.com/p/2483520076 按照上面这个帖子的方法，确定参数如下： sigma：3.405e-10 m epsilon： 1.67e-21 J 玻尔兹曼常数Kb：1.38065E-23 J/K 摩尔质量：39.948 gram/mol 单个原子质量：6.63352E-26 kg (用摩尔质量/阿伏伽德罗常数) 具体代码做了如下修改，lj部分为源代码，metal部分为我自己修改后的，源代码见lammps安装包：lammps/Examples/VISCOSITY/in.mp.2d lj: variable x equal 20 #x方向长度 variable y equal 20 #y方向长度 metal: 保持不变 lj: variable rho equal 0.6 #lj单位中的约化密度（相当于metal单位中的晶格常数） variable t equal 1.0 #温度 variable rc equal 2.5 #截断半径 metal: variable rho equal 6.217 #晶格常数 (由rho*=rho•sigma^dim，和rho=N/V； 二维模拟，dim=2，采用sq2结构所含原子数N=2，V=a^2，a为晶格常数； 得0.6=2•(sigma/a) ^2，所以a=6.217 埃) variable t equal 120.96 #温度 (由T*=T•Kb/epsilon ，1.0*1.67e-21/(1.38065E-23)=120.96 K) variable rc equal 8.5125 #截断半径 (由x*=x/sigma，2.5*3.405e-10 =8.5125埃) lj: units lj dimension 2 atom_style atomic neigh_modify delay 0 every 1 lattice sq2 ${rho} region simbox block 0 $x 0 $y -0.1 0.1 create_box 1 simbox create_atoms 1 box metal: 保持不变 lj: pair_style lj/cut ${rc} pair_coeff * * 1 1 metal: pair_style lj/cut ${rc} pair_coeff * * 1.042332128634570e-2 3.405 (势阱常数1.042332128634570e-2 eV，零势能距离 3.405 埃) lj: mass * 1.0 timestep 0.005 velocity all create $t 97287 metal: mass * 39.948 (质量改为相对摩尔质量) timestep 0.01 （lj单位默认时间步长0.005，由t* = t (epsilon / (m•sigma^2)^1/2，得lj下默认时间为10fs，即0.01ps) velocity all create $t 97287 lj: fix 1 all nve fix 2 all langevin $t $t 0.1 498094 fix 3 all enforce2d metal: fix 1 all nve fix 2 all langevin $t $t 0.2 498094 (由damp = damping parameter (time units)) fix 3 all enforce2d lj: # equilibration run thermo 1000 run 5000 unfix 2 # turn on Muller-Plathe driving force and equilibrate some more velocity all scale $t fix 4 all viscosity 100 x y 20 fix 5 all ave/spatial 20 50 1000 y center 0.05 vx & units reduced file profile.mp.2d # equilibration run variable dVx equal f_5[11][3]-f_5[1][3] thermo 1000 thermo_style custom step temp epair etotal press f_4 v_dVx run 20000 # data gathering run # reset fix viscosity to zero flux accumulation unfix 4 fix 4 all viscosity 100 x y 20 #100步交换一次动量，交换x方向的动量分量，另一个分量在y方向，共20层 metal: 保持不变 lj: variable visc equal -(f_4/(2*(step*0.005-125)*lx+1.0e-10))/(v_dVx/(ly/2)) #http://lammps.sandia.gov/threads/msg43465.html #http://www.52souji.net/lammps-command-fix-viscosity/ #f_4输出总的动量交换量，除以时间和截面积即为动量通量； #(step*0.005-125)为时间，0.005为lj单位默认的时间步长，125为到目前为止（20000+5000=25000）步的时间：25000*0.005=125，为防止25000步计算时分母为0,加上1.0e-10； #lx为截面积（二维，z方向为1）； #v_dvx为速度分量， 除以ly/2为速度梯度 variable visc equal -((f_4/(2*(step*0.005-125)*lx+1.0e-10))/(v_dVx/(ly/2)))*9.07812e-05 metal: variable visc equal -((f_4/(2*(step*0.01-250)*lx+1.0e-10))/(v_dVx/(ly/2)))*1.66e-5 lj: fix vave all ave/time 1000 1 1000 v_visc ave running start 26000 #ave running 每1000步输出一次到目前为止的v_visc的累积平均值，作为vave的值,从26000步开始 thermo_style custom step temp f_4 v_dVx v_visc f_vave run 50000 metal: 保持不变 |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com - 附件 1 : lj.docx
- 附件 2 : mp.rar
2015-04-11 19:20:29, 111.74 K
2015-04-11 23:41:38, 35.95 K
» 猜你喜欢
有多少人是今天查系统知道结果的？
已经有17人回复
-
为什么资助数各大高校都创新高，自己申请怎么就这么难
已经有8人回复
-
导师吐槽：我怎么摊上了这么个极品研究生！
已经有7人回复
-
青基如何破局
已经有7人回复
-
基金不中，共勉
已经有10人回复
-
国自然面上复盘~欢迎讨论
已经有11人回复
-
基金未中，这种答复是模板吗？
已经有6人回复
-
看板上这么多中的，有点像50人群里49个人都是骗子的那种感觉……
已经有6人回复
-
怎么看青基中了没有啊
已经有5人回复
-
梦想
已经有3人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- lammps中输出项Enthalpy的单位是什么
已经有8人回复
2楼2015-04-15 10:50:43
3楼2015-04-26 18:40:50
邹卓民
新虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 602.2
- 散金: 18
- 帖子: 209
- 在线: 124.9小时
- 虫号: 4881684
- 注册: 2016-07-30
- 性别: GG
- 专业: 交通工程
4楼2018-07-03 20:23:42
5楼2018-10-14 15:56:24
|
本帖内容被屏蔽 |
6楼2018-10-15 09:35:23
7楼2018-10-18 13:55:57
8楼2020-06-16 17:17:49
