使用Matlab绘制LFP锂离子扩散动态示意图

今天分享如何运用Matlab编程实现LFP锂离子扩散动态示意图绘制。

一、前期准备工作（LFP简介+锂离子扩散机制+Matlab动画）

1.1 LFP简介

磷酸锂铁（分子式LiFePO4，Lithium Iron Phosphate ，按IUPAC命名法应为磷酸亚铁锂，简称LFP），是一种锂离子电池的正极材料，特色是地球资源丰富，价格低，不含钴等贵重元素。其工作电压适中（3.2V）、平台稳定、比容量高（170mAh/g）、高放电功率、可快速充电且循环寿命长，热稳定性高。

根据LFP晶体CIF文件，通过上一期分享教程可以绘制如下晶体结构示意图。

1.2 锂离子扩散机制

大量实验研究和理论计算表明，LFP晶体中锂离子主要通过一维通道实现固相扩散。第一性原理计算结果显示锂离子沿着[010]方向扩散活化能约0.27eV，远小于其他的锂离子扩散路径。如下图所示，我们可以理解到锂离子从原位置扩散至相邻锂空位，这个扩散过程按经典理论考虑，就是锂离子的扩散路径必然是沿着起点和终点之间最低的势能面上。因此，我们得到了模拟锂离子扩散过程的一个重要假设条件，即锂离子的扩散路径必然是沿着最低的势能面。

左图：Crystal structure of LiFePO4 and possible lithium pathways.

右图：Anisotropic harmonic lithium vibration in LiFePO4 shown as green thermal ellipsoids and the expected diffusion path.

参考文献：Nishimura, Si., Kobayashi, G., Ohoyama, K. et al. Experimental visualization of lithium diffusion in LixFePO4. Nature Mater 7, 707–711 (2008).

1.3 Matlab动画绘制

Matlab的动画制作主要有三种方式：

质点动画：最简单的动画产生方式，产生一个顺着曲线轨迹运动的质点来操作，使用comet、comet3函数，前者是二维，后者是三维；

电影动画：首先保存一系列的图形，然后按照一定的顺序像电影一样的播放，通过getframe函数将当前的图片抓取为电影的画面，再由movie函数将动画显示出来。

程序动画：在图形窗口中按照一定的算法连续擦除和重绘图形对象，可以采用重绘图形对象的方法来创建程序动画。改变对象的方法可以触发MATLAB对该对象进行重绘。

下面，我们重点研究如何使用Matlab建模，简单模拟锂离子在LFP晶体中沿着[010]方向扩散的动态示意图。

二、晶胞绘制数学模型搭建

2.1 模型假设

为了简化模型，我们对LFP中锂离子的扩散过程做如下合理的假设：

1）锂离子沿着[010]方向扩散，不考虑其他扩散路径；

2）由于Li与其最邻近的6个氧原子构成LiO6六面体，其他例如Fe，P较远，考虑简化计算复杂度，因此在计算势能时候只考虑最邻近O原子的作用；

3）Li与氧原子之间作用包括引力和斥力，其原子间相互作用势近似采用Lennard-Jones势能函数表示；

$\varphi =4\varepsilon ([\sigma /r]^{12}-[\sigma /r]^{6})$

其中，ε表示势阱深度，σ表示势能耗散。

由于未查到关于Li-O相互作用势能函数中对应参数，本实例中，通过多次调参，最终确定ε=1，σ=0.8。

4）不考虑扩散过程速度变化，假定锂离子沿[010]方向速度分量一致。

2.2 模型搭建

为了绘制锂离子在LFP晶体中扩散动画，首先必须获取锂离子在LFP晶体中扩散路径，通过扩散路径确定动画中每一帧对应的锂离子轨迹，从而实现动态过程的呈现。

如上图所示，我们需要计算1号锂原子到3号锂原子之间扩散轨迹（左图），因此问题可以简化为在锂原子1与锂原子3之间插入一系列垂直于[010]方向的平面，并计算锂原子在此平面上最低势能点处对应位置坐标（右图）。

三，Matlab函数代码

第一步：确定中心原子即1,2,3号锂原子的坐标Li1,Li2,Li3;

%% 绘制锂原子扩散示意图
%%导入晶体基本结构信息
clc;clear;
[cellpara,position,site,spacegroup]=cifread('lfp.cif');
siteposition = sitepositionfun(cellpara,position,site,spacegroup);
%%确定123锂原子坐标
Li2=[0.5 0.5 0.5].*cellpara(1:3)';
Li1=[0 0.5 0.5].*cellpara(1:3)';
Li3=[1 0.5 0.5].*cellpara(1:3)';

第二步，确定与中心Li原子相邻6个氧原子坐标，通过计算所有O原子与中心锂原子的距离，筛选其中距离最小的六个O原子，通过计算可得6个氧原子对应坐标；这里采用了一个偷懒技巧，应为Li-O成键键长不超过2.2Å，所以这里直接通过长度判断，找出相邻氧原子。

%%寻找与Li2最紧邻的6个氧原子
[len,wid]=size(siteposition);
sitelen=cell(len,wid);O2site=zeros(6,3);
kplus=0;
for k=1:lensitelen{k,1}=strcat('Li-',siteposition{k,1});[m,n]=size(siteposition{k,2}(:,:));for i=1:msitelen{k,2}(i,1)=norm(siteposition{k,2}(i,:)-0.5*cellpara(1:3)');if sitelen{k,2}(i,1)<=2.2&&contains(siteposition{k,1},'O')kplus=kplus+1;O2site(kplus,:)=siteposition{k,2}(i,:);endend
end

第三步，构建超晶胞，采用第二步类似方法，寻找分别与Li1和Li3最近邻的6个氧原子的位置;同时合并1，2，3号锂原子存在相互作用的所有氧原子，并剔除重复原子；

%%计算Li1 Li3最近邻氧原子位置
%%构建[3 1 1]超晶胞，并寻找与Li1 Li3最紧邻氧原子
%%绘制[3 1 1]超晶胞
figure(1)
down=[-1 0 0];
up=[2 1 1];
supercell=[3 1 1];
pos=plotsupercell([3 1 1],siteposition,cellpara,[-1 0 0],[2 1 1]);[len,wid]=size(pos);
sitelen1=cell(len,wid);O1site=zeros(6,3);
sitelen3=cell(len,wid);O3site=zeros(6,3);
kplus1=0;kplus3=0;
for k=1:lensitelen1{k,1}=strcat('Li1-',pos{k,1});sitelen3{k,1}=strcat('Li3-',pos{k,1});[m,n]=size(pos{k,2}(:,:));for i=1:msitelen1{k,2}(i,1)=norm(pos{k,2}(i,:)-Li1);sitelen3{k,2}(i,1)=norm(pos{k,2}(i,:)-Li3);if sitelen1{k,2}(i,1)<=2.2&&contains(pos{k,1},'O')kplus1=kplus1+1;O1site(kplus1,:)=pos{k,2}(i,:);elseif sitelen3{k,2}(i,1)<=2.2&&contains(pos{k,1},'O')kplus3=kplus3+1;O3site(kplus3,:)=pos{k,2}(i,:);endend
end
%%合并所有氧原子位，剔除重复项
Osite=[O1site;O2site;O3site];
Osite=uniquetol(Osite,'Byrows',true);

第四步，根据锂离子扩散轨迹分解方法，计算锂原子处于任意中间位置(xk yk zk)时，对应Li-O原子相互作用势能总和；

$\Phi = \sum_{k=1}^{N}\varphi (k)$

其中，ψ(k)表示Li与第k个氧原子之间相会作用势；N表示总计氧原子数。

求解Φ函数最小值，并得到相应Li位置坐标(xk yk zk)；在求解Φ函数最小值时使用了matlab内部优化算法GlobalSearch；

%%求解中间轨迹锂原子位置
%%设置中间轨迹数，构建最低势能函数方程
N=10;trace=zeros(N,3);for i=1:Ntrace(i,1)=0*(N-i)/(N+1)+cellpara(1)*(i)/(N+1);%%势能函数epsilon=1;sigma=0.85;phifun=@(x) sum(4*epsilon*((sigma/sqrt((trace(i,1)-Osite(:,1)).^2+(x(1)-Osite(:,2)).^2+(x(2)-Osite(:,3)).^2)).^12-...(sigma/sqrt((trace(i,1)-Osite(:,1)).^2+(x(1)-Osite(:,2)).^2+(x(2)-Osite(:,3)).^2)).^6));x0=[0.5*cellpara(2) 0.5*cellpara(3)];lb=[min(Osite(:,2)),min(Osite(:,3))];ub=[max(Osite(:,2)),max(Osite(:,3))];rng default % For reproducibilitygs = GlobalSearch;options=[];problem = createOptimProblem('fmincon','x0',x0,...'objective',phifun,'lb',lb,'ub',ub,'options',options);[xval,~,exitflag,~] = run(gs,problem);trace(i,2:3)=xval;
end

第五步，在晶胞中绘制中心锂原子对应的轨迹，检验计算结果是否合理；从结果来看，锂原子在1-3号之间扩散确实时曲线方式进行，但是由于原子相会作用势函数方程中参数不准，导致计算运动轨迹并非真实或者理论，但作为简单模拟应该是ok的。

%%在晶胞中绘制锂原子轨迹示意图
figure(2)
plotunitcell(siteposition,cellpara);
for i=1:length(trace(:,1))plotsphere(trace(i,:),0.2,'y');
end
saveas(figure(2),'trace.jpg');%%绘制原子函数
function plotsphere(pos,r,color)
[u,v,w]=sphere(100);
x=pos(1)+r.*u;
y=pos(2)+r.*v;
z=pos(3)+r.*w;
surf(x,y,z,'FaceColor',color,'EdgeColor','none','FaceAlpha',0.25);
colormap hot;
grid off;
end

第六步，根据每一帧下锂原子轨迹坐标计算锂原子坐标变化，通过平移原理，重复计算每一帧下对应所有原子轨迹坐标；最后根据获取每一帧对应所有原子轨迹坐标，绘制晶体示意图，通过getframe函数录制每一帧下图像，并通过imwrite函数保存为gif动画。

%%计算锂原子delta轨迹
dtrace=trace-Li1;
%%计算晶胞中所有原子轨迹
postrace=cell(len,wid,3*(N+1)-1);
postrace(:,:,1)=pos;
k=1;
for i=1:3for j=1:Nk=k+1;postrace(:,:,k)=pos;temp=pos{1,2}(:,1:3)+(i-1)*[cellpara(1) 0 0]+dtrace(j,:);temp(temp(:,1)>2*cellpara(1),:)=[];postrace{1,2,k}=temp;end
endfor i=1:30figure(3)plotsupercell(supercell,postrace(:,:,i),cellpara,down,up);axis([-cellpara(1),2*cellpara(1),-0.05*cellpara(2),1.05*cellpara(2),-0.05*cellpara(3),1.05*cellpara(3)]);light;CF=getframe(gcf);imind=frame2im(CF);[imind,cm] = rgb2ind(imind,256);if i==1imwrite(imind,cm,'test.gif','gif', 'Loopcount',inf,'DelayTime',1e-4);elseimwrite(imind,cm,'test.gif','gif','WriteMode','append','DelayTime',1e-4);endhold off
end

四，代码测试及演示

回顾晶体结构示意图绘制结合动画制作，可以实现动态视角观察晶体结构示意图。

下面我们先看一下锂离子扩散轨迹示意图；从图可以看出锂离子确实以近似sina函数曲线方式沿着[010]方向扩散，基本与文献报道相似。不过这里值得注意的是锂离子

那么，展示一下锂离子扩散动态示意图。

最后顺带提一下，希望晶体结构示意图绘制技巧能帮助大家在做PPT演示的时候用上，让大家收获老板的青睐，升职加薪走上人生巅峰

使用Matlab绘制LFP锂离子扩散动态示意图相关推荐

matlab 绘制椭圆锥波束指向示意图
参考了matlab 绘制任意方向和位置的空间圆锥体中绘制圆锥的思路,在此基础上加了我需要的应用.就是绘制波束扫描示意图. 1 原理根据波束宽度得到在归一化威力值上的弧长. 然后根据方位和俯仰向弧长产 ...
matlab 绘制三阶魔方-动态变化
三阶魔方绘制-动态变化魔方绘制魔方绘制可以参考链接: [手把手制作三阶魔方模拟器]用MATLAB绘制一个三阶魔方和 [手把手制作三阶魔方模拟器]用MATLAB让你的魔方动起来. 可以达到的效果为: ...
MATLAB绘制主函数动态图,matlab绘制动态图
mathematica绘制动态图,"绘图之王"争霸赛--Excel才是绘图王道,matlab绘制动态图,动态三维图绘制 matlab动态图画法_数学_自然科学_专业资料.Matla ...
使用matlab绘制示意图的示例
在本示例中,将给出如何用matlab绘制示意图,添加图例,添加标题等等,以及在matlab中调用LaTeX字体的命令. 代码如下: clear;clc;close all% Initialize sy ...
Matlab 绘制动态图
在写论文的过程中,我们经常需要用MATLAB绘制图形.论文中的图形都是图片格式,但是在展示和汇报时,如果将图形做成GIF动图,变量的变化过程就会非常直观,展示效果也会非常好.下面将本人利用MATLAB ...
matlab系统的根轨迹,实验五利用MATLAB绘制系统根轨迹
<实验五利用MATLAB绘制系统根轨迹>由会员分享,可在线阅读,更多相关<实验五利用MATLAB绘制系统根轨迹(6页珍藏版)>请在人人文库网上搜索. 1.实验五利用MAT ...
【转】MATLAB如何制作动画（动态图形演示movie）
https://jingyan.baidu.com/article/49711c6199002dfa451b7c52.html MATLAB如何制作动画(动态图形演示movie) 听语音 | 浏览:7 ...
利用matlab绘制简单IFS图形（Sierpinski三角形和BarnsleyFern巴恩斯利蕨）
利用matlab绘制简单IFS图形(Sierpinski谢尔宾斯基三角形和BarnsleyFern巴恩斯利蕨) 一.SierpinskiTriangle谢尔宾斯基三角形谢尔宾斯基三角形(英语:Sie ...
matlab中s_cplot,matlab系统模型建立和动态特性研究分析实验.docx
实验二MATLAB系统模型建立和动态特性分析实验一.实验目地 1掌握如何使用 MALAB进行系统模型地建立: 2 ?学习利用MALAB命令得阶跃响应曲线,分析系统动态特性; 3.利用MALAB求阶跃 ...
matlab绘制双缝、等厚劈尖干涉、牛顿环、迈克尔逊等倾干涉图样
matlab绘制双缝.等厚劈尖干涉.牛顿环.迈克尔逊等倾干涉图样 1.matlab绘制双缝干涉图样 2.matlab绘制等厚劈尖干涉图样 3.matlab绘制牛顿环干涉图样 4.matlab绘制迈克尔 ...

使用Matlab绘制LFP锂离子扩散动态示意图

使用Matlab绘制LFP锂离子扩散动态示意图相关推荐

最新文章

热门文章