计算机在材料科学中的应用上机教程.doc

计算机在材料科学中的应用上 机 教 程（2011级高分子专业用）湖北工业大学材料科学与工程学院2013年2月计算机在材料科学中的应用上机教程 李纯清前 言当前，计算机在材料科学中的应用(compnter applicationln materials science)已经发展成为材料科学与工程的一个新的学科分支，涉及很多的学科和领域，需要把材料科学的专业知识与计算机技术、计算数学、信号与系统、物理、图形图像学、自动控制技术、机械等各方面的知识技术进行综合应用。计算机在材料科学中的应用是指以计算机为手段，通过理论和计算对材料的固有性质、结构与组分、使用性能及合成与加工进行综合研究的一门新学科方向，其目的在于使人们能主动地对材料进行结构、功能与工艺的优化与控制，以使按需要开发及制备新材料。把计算机应用于材料科学的思想产生丁20世纪50年代，其形成为一门独立的新兴学科则是20世纪80年代以后的事，在我国20世纪90年代逐步有学者进行这方面的研究尝试工作。近年来，现代科学(信息技术、量子力学、统计物理、固体物理、量子化学、计算科学、计算机图形学等)理论和方法技术的飞速发展，以及计算机能力的空前提高为材料计算机模拟提供了理论基础和有力手段。计算机在材料领域中的应用的发展将使材料科学从半经验的定性描述逐渐进入半定量到定量预测和控制的更为科学的阶段。材料的计算机模拟技术及应用已经成为现代材料科学研究中最为活跃的一个重要分支。目前计算机在材料领域中的应用在国际上还没有统一的通用术语，美国习惯称其为材料的“计算机分析与模型化”(computer-based analysis and modeling)，欧洲则称它为“计算材料科学”(computational material science)，日本则称其为“材料设计”(material design)，等等。尽管在用语上不完全一样，各国在研究领域和特色上也不尽相同，但基本含义是相同的。高技术新材料是现代知识经济的重要组成部分，其发展日新月异，这为材料的计算与设计提供了发展机遇和广阔空间。计算机在材料领域中的应用具有以下特点：(1)具有“前瞻性”；(2)具有“创新性”s；(3)可减少或替代实验室工作；(4)降低研究及生产成本。目前，随着材料科学研究工作的广泛开展和不断深入，有关“计算机在材料科学中的应用”的科技文献资料迅速增多。但这些文献大都出现在学术期刊或会议文集之中，还缺少既包括基本原理也包括最新研究成果的系统性、综合性书籍，教材(面向本科生、研究生)就更少。但21世纪新材料在国民经济建设中的地位越来越重要，虽然目前我国在这方面出版的教材及参考书还不多，但我们还是决定开设此门课程，目的是要学生了解材料科学的新进展，让更多人接触这一新的发展方向，缩小与国外的差距。本上机教程的教学目标如下：1、使材料科学与工程专业的学生了解计算机在材料科学中应用的现状和发展趋势，知道所要从事的材料科学研究领域中应当掌握哪些计算机技术和知识；2、通过本门上机教程程的学习引导，使他们初步掌握一些计算机地材料科学中应用的入门知识和技能，在面对这个新领域时能比较轻松；3、通过本门上机教程介绍计算机应用技术和软件在材料科学中的应用，希望提高材料科学工作者从事材料研究工作时的效率利和水平。由于时间匆忙，业务能力有限，书中难免有误，请读者批评指正，编者万分感谢!77目 录前 言i目 录iii上篇 数据处理及过程模拟1实验01：实验数据及模型参数拟合方法1单变量拟合1例1-1：1例1-22多变量拟合3例1-33例1-44例1-55例1-66例1-76例1-87实验02：非线性方程求解9例2-1 用对分法求f(x)=x3-7.7*x2+19.2*x-15.3在区间【1，2】之间的根。9例2-2求代数方程x3-2*x-5=0,在x0=2附近的根。10例2-310例2-4用牛顿迭代法求方程在x0=1附近的零点。11例2-5用割线法求方程的根（取x0=1.5,x1=4）。12例2-6：用牛顿迭代法解下列非线性方程组13实验03：线性方程组的迭代求解15简单迭代计算机算法15例3-1 用简单迭代格式解下列方程组15紧凑迭代计算机算法16例3-2 用简单迭代格式解下列方程组16松弛迭代计算机算法18例3-3 用简单迭代格式解下列方程组18实验04 基于matlab的微分方程数值解20例4-3 用庞格-库塔公式，求解下列初值问题。20例4-4：微分方程组数值解21实验作业：教材page 88 第5题21实验05 偏微分方程数值解22例5-1：用数值法求解下列偏微分方程：22例5_1:22例5_2: (第五章课后习题2)30实验06： 有限元法解微分方程34中篇 分析应用篇36实验07：泡沫材料泡孔（颗粒）大小及形状分析36实验08：简单流动polyflow模拟与分析（课堂讲授）41实验09：聚合物熔体与口模间传热的非等温流动polyflow模拟与分析42实验10：2.5d轴对称挤出polyflow模拟与分析43实验11：chemoffice2008在材料科学中的应用441、分子构象位垒442、红外光谱预测453、核磁共振谱预测514、计算过渡态能量55实验12：matrials studio在材料科学中的应用581. 建立全同立构pmma。582、聚合物与金属氧化物表面的相互作用60实验13：office、origin等在材料科技论文写作中的应用70实验14：材料科学中matlab图像处理基础71实验15：surface evolver在材料科学研究中的应用73计算机在材料科学中的应用上机教程 李纯清上篇 数据处理及过程模拟实验01：实验数据及模型参数拟合方法实验平台：matlab 2012b实验目的：lsqcurvefit的使用；实验原理：lsqcurvefit solves non-linear least squares problems.lsqcurvefit attempts to solve problems of the form:min sum (fun(x,xdata)-ydata).2 where x, xdata, ydata and the values returned by fun can be vectors or matrices. x = lsqcurvefit(fun,x0,xdata,ydata) starts at x0 and finds coefficients x to best fit the nonlinear functions in fun to the data ydata (in the least-squares sense). fun accepts inputs x and xdata and returns a vector (or matrix) of function values f, where f is the same size as ydata, evaluated at x and xdata. note: fun should return fun(x,xdata) and not the sum-of-squares sum(fun(x,xdata)-ydata).2). (fun(x,xdata)-ydata) is squared and summed implicitly in the algorithm.) 实验举例：单变量拟合例1-1：x=7:2:47 x = 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 y=9:3:69 y = 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 定义函数模型：f1=inline(a(1)+a(2)*x,a,x) f1 = inline function: f1(a,x) = a(1)+a(2)*x a1,resnormal1=lsqcurvefit(f1,1,1,x,y) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a1 = -1.5000 1.5000resnormal1 = 7.0682e-28 例1-2定义二次函数模型f2=inline(a(1)+a(2)*x+a(3)*x.2,a,x) f2 = inline function: f2(a,x) = a(1)+a(2)*x+a(3)*x.2 x2=-3:3 x2 = -3 -2 -1 0 1 2 3 y2=4 2 3 0 -1 -2 -5 y2 = 4 2 3 0 -1 -2 -5 a2,resnormal2=lsqcurvefit(f2,0.5, -1, 0,x2,y2) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a2 = 0.6667 -1.3929 -0.1310resnormal2 = 3.0952 多变量拟合例1-3x3=100 200 300 500 100 700 800;2 4 1 0.3 5 3 4 x3 = 100.0000 200.0000 300.0000 500.0000 100.0000 700.0000 800.0000 2.0000 4.0000 1.0000 0.3000 5.0000 3.0000 4.0000 y3=1.127 2.416 2.205 2.312 1.484 6.038 7.325 y3 = 1.1270 2.4160 2.2050 2.3120 1.4840 6.0380 7.3250 定义传热实验模型函数f3=inline(c(1)*x(1,:).c(2).*x(2,:).c(3),c,x) f3 = inline function: f3(c,x) = c(1)*x(1,:).c(2).*x(2,:).c(3) a3,resnormal3,residual3=lsqcurvefit(f3,0.02 0.5 0.2,x3,y3) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a3 = 0.0230 0.8000 0.3000resnormal3 = 2.6650e-07residual3 = 1.0e-03 *0.1358 0.1653 -0.1808 -0.0250 -0.2256 0.3005 -0.2151 例1-4利用解矛盾方程的方法，用二次多项式函数拟合以下数据.序号xy序号xy1-3451-12-2262-23-1373-5400a=1 -3 9;1 -2 4;1 -1 1;1 0 0;1 1 1;1 2 4 ;1 3 9 a = 1 -3 9 1 -2 4 1 -1 1 1 0 0 1 1 1 1 2 4 1 3 9 y=4 2 3 0 -1 -2 -5 y = 4 2 3 0 -1 -2-5 inv(a*a)*(a*y) ans = 0.666666666666667 -1.392857142857143 -0.130952380952381 例1-5x5=-3 -2 -1 2 4 x5 = -3 -2 -1 2 4 y5=14.3 8.3 4.7 8.3 22.7 y5 = 14.3000 8.3000 4.7000 8.3000 22.7000 f5=inline(a(1)+a(2)*x.3,a,x) f5 = inline function: f5(a,x) = a(1)+a(2)*x.3 a5 resnorm5=lsqcurvefit(f5,10,0,x5,y5) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a5 = 10.6751 0.1368resnorm5 = 112.0147 例1-6x6=1 2 3 4 5 6 7 8 x6 = 1 2 3 4 5 6 7 8 y6=15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6 y6 = 15.3000 20.5000 27.4000 36.6000 49.1000 65.6000 87.8000 117.6000 f6=inline(a(1)*exp(a(2)*x),a,x) f6 = inline function: f6(a,x) = a(1)*exp(a(2)*x) a6,resnormal6=lsqcurvefit(f6,10 ,0,x6,y6) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a6 = 11.4241 0.2914resnormal6 = 0.0119 例1-7x7=1 1 2 3;1 3 5 -1;1 -1 2 5 x7 = 1 1 2 3 1 3 5 -1 1 -1 2 5 y7=2 -1 1 -2 y7 = 2 -1 1 -2 f7=inline(a(1)*x(1,:)+a(2)*x(2,:)+a(3)*x(3,:),a,x) f7 = inline function: f7(a,x) = a(1)*x(1,:)+a(2)*x(2,:)+a(3)*x(3,:) a7,resnormal7=lsqcurvefit(f7,0,0,0,x7,y7) local minimum found.optimization completed because the size of the gradient is less thanthe default value of the function tolerance.a7 = -1.5917 0.5899 0.7572resnormal7 = 5.4712 例1-8x8=-23.7 -10 0 10 20 30 40+273.15 x8 = 249.4500 263.1500 273.1500 283.1500 293.1500 303.1500 313.1500 y8=0.101 0.174 0.254 0.359 0.495 0.662 0.880/0.101325*760 y8 = 1.0e+03 * 0.7576 1.3051 1.9052 2.6927 3.7128 4.9654 6.6005 fun8=inline(exp(a(1)-a(2)./(x+a(3),a,x) fun8 = inline function: fun8(a,x) = exp(a(1)-a(2)./(x+a(3) a8,resnormal8=lsqcurvefit(fun8,14,2000,-15,x8,y8) local minimum possible.lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance.a8 = 1.0e+03 * 0.0172 2.5921 -0.0035resnormal8 = 842.1467 y88=exp(a8(1)-a8(2)./(x8+a8(3) y88 = 1.0e+03 * 0.7544 1.3155 1.9048 2.6861 3.6991 4.9865 6.5934 (y88-y8)./y8 ans = -0.0042 0.0079 -0.0002 -0.0024 -0.0037 0.0042 -0.0011 实验02：非线性方程求解例2-1 用对分法求f(x)=x3-7.7*x2+19.2*x-15.3在区间【1，2】之间的根。方法1：x=solve(x3-7.7*x2+19.2*x-15.3,x) x = 1.7 3.0 3.0 方法2：新建一个script文件，文件名为half2root.m，输入如下语句。clc;clear all;a=1;b=2;fa=a3-7.7*a2+19.2*a-15.3;fb=b3-7.7*b2+19.2*b-15.3;c=(a+b)/2;fc=c3-7.7*c2+19.2*c-15.3;while abs(fc)0.5*10(-4) c=(a+b)/2; fc=c3-7.7*c2+19.2*c-15.3; if fb*fc0 b=c; fb=fc; else a=c; fa=fc; endendformat longfx=fc, x=c half2root fx = 2.062949538128578e-05x = 1.700012207031250 例2-2求代数方程x3-2*x-5=0,在x0=2附近的根。方法1：fsolve(x3-2*x-5,2) equation solved.fsolve completed because the vector of function values is near zeroas measured by the default value of the function tolerance, andthe problem appears regular as measured by the gradient.ans = 2.094551481698583 方法2：例2-3方法1：先建立文件名为nonlineequation.m的脚本文件。文件内容如下：function f=nonlineequation(x)f=0.02*x(1)2-x(1)+0.1*x(2)2+1;0.1*x(1)2-x(2)+0.01*x(2)3+2;fsolve(x)nonlineequation(x),0,0) equation solved.fsolve completed because the vector of function values is near zeroas measured by the default value of the function tolerance, andthe problem appears regular as measured by the gradient.ans = 1.6319 2.4055 方法2:例2-4用牛顿迭代法求方程在x0=1附近的零点。%newton_method_solve_equation.mclear allsyms xf=x3-7.7*x2+19.2*x-15.3x0=1err=1while(err0.00001) x_tmp=x0-subs(f,x,x0)/subs(diff(f,x,1),x,x0); err=norm(x_tmp-x0,inf) x0=x_tmp endx0errnewton_method_solve_equation f =x3 - (77*x2)/10 + (96*x)/5 - 153/10x0 = 1err = 1err = 0.4118x0 = 1.4118err = 0.2115x0 = 1.6232err = 0.0691x0 = 1.6923err = 0.0076x0 = 1.6999err = 8.9619e-05x0 = 1.7000err = 1.2357e-08x0 = 1.7000x0 = 1.7000err = 1.2357e-08 例2-5用割线法求方程的根（取x0=1.5,x1=4）。%chord_method_solve_equation.mclear all;syms x;f=x3-7.7*x2+19.2*x-15.3x0=1.5x1=3err=1;while(err0.00001) x_tmp=x1-subs(f,x,x1)*(x1-x0)/(subs(f,x,x1)-subs(f,x,x0); err=norm(x_tmp-x1,inf); x0=x1; x1=x_tmp; endx1err f =x3 - (77*x2)/10 + (96*x)/5 - 153/10x0 = 1.500000000000000x1 = 3x1 = 2.999999999999988err = 1.199040866595169e-14 例2-6：用牛顿迭代法解下列非线性方程组%newton_method_solve_equation_group.mclear all;syms x y;f1=4-x2+y2;f2=1-exp(x)-y;jac=diff(f1,x) diff(f1,y);diff(f2,x) diff(f2,y);xy0=1 -1.7;err=1;while(err0.00001) xy_tmp=-inv(subs(jac,x y,xy0)*subs(f1 f2,x y,xy0); err=norm(xy_tmp,inf); xy0=xy0+xy_tmp; endxy0errnewton_method_solve_equation_group xy0 = -2.1882 0.8879err = 2.4786e-06 实验03：线性方程组的迭代求解简单迭代计算机算法简单迭代公式：例3-1 用简单迭代格式解下列方程组clc clear all a=2 -1 -1;1 5 -1;1 1 10 a = 2 -1 -1 1 5 -1 1 1 10 b=-5 8 11 b = -5 8 11 d=diag(diag(a) d = 2 0 0 0 5 0 0 0 10 b=eye(3)-inv(d)*a b = 0 0.5000 0.5000 -0.2000 0 0.2000 -0.1000 -0.1000 0 f=inv(d)*b f = -2.5000 1.6000 1.1000 x0=1 1 1 x0 = 1 1 1 x0=b*x0+f (重复多次)x0 = -1.0000 2.0000 1.0000 紧凑迭代计算机算法紧凑迭代公式：例3-2 用简单迭代格式解下列方程组clc clear all a=2 -1 -1;1 5 -1;1 1 10 a = 2 -1 -1 1 5 -1 1 1 10 b=-5 8 11 b = -5 8 11 d=diag(diag(a) d = 2 0 0 0 5 0 0 0 10 l=-tril(a,-1) l = 0 0 0 -1 0 0 -1 -1 0 u=-triu(a,1) u = 0 1 1 0 0 1 0 0 0 g=(d-l)-1*u g = 0 0.5000 0.5000 0 -0.1000 0.1000 0 -0.0400 -0.0600 f=(d-l)-1*b f = -2.5000 2.1000 1.1400 x0=1 1 1 x0 = 1 1 1 x0=g*x0+f （重复多次）x0 = -1.0000 2.0000 1.0000 松弛迭代计算机算法紧凑迭代公式：例3-3 用简单迭代格式解下列方程组clc clear all a=2 -1 -1;1 5 -1;1 1 10 a = 2 -1 -1 1 5 -1 1 1 10 b=-5 8 11 b = -5 8 11 omiga=1.1 omiga = 1.1000 d=diag(diag(a) d = 2 0 0 0 5 0 0 0 10 l=-tril(a,-1) l = 0 0 0 -1 0 0 -1 -1 0 u=-triu(a,1) u = 0 1 1 0 0 1 0 0 0 g=(d-omiga*l)-1*(1-omiga)*d+omiga*u) g = -0.1000 0.5500 0.5500 0.0220 -0.2210 0.0990 0.0086 -0.0362 -0.1714 f=omiga*(d-omiga*l)-1*b f = -2.7500 2.3650 1.2524 x0=1 1 1 x0 = 1 1 1 x0=g*x0+f （重复多次）x0 = -1.0000 2.0000 1.0000 实验04 基于matlab的微分方程数值解例4-3 用庞格-库塔公式，求解下列初值问题。本实验用matlab中的相关命令来求解。1、求理论解 y=dsolve(dy-y2*cos(x),y(0)=1,x) y =-1/(sin(x) - 1) syms x subs(y,x,0.1:0.1:0.8) ans = 1.1109 1.2479 1.4195 1.6378 1.9210 2.2970 2.8107 3.5380 2、求数值解第一种方法：先定义m函数文件，如下：function li=lichun(x,y)li=y2*cos(x);再在windows command下执行如下语句：x1,y1=ode45(lichun,0, 0.8,1) x2,y2=ode45(lichun,0.0:0.1:0.8,1) 第二种方法：直接在wondows command下直接执行如下语句：clear all; %定义内嵌函数f=inline(y2*cos(x),x,y) f = inline function: f(x,y) = y2*cos(x) x1,y1=ode45(f,0.0, 0.8,1) x2,y2=ode45(f,0.0:0.1:0.8,1) 例4-4：微分方程组数值解先定义名为exam4_4.m的文件，其内容如下：function fun=exam4_4(t,y)fun=0.09*y(1)*(1-y(1)/20)-0.45*y(1)*y(2); 0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2)然后在matlab的command window下输入以下命令：clear all ;t,y=ode45(exam4_4,0:0.5:4,1.6,1.2) ; t ,y ans = 0 1.6000 1.2000 0.5000 1.2689 1.2327 1.0000 0.9995 1.2663 1.5000 0.7817 1.3010 2.0000 0.6068 1.3366 2.5000 0.4674 1.3732 3.0000 0.3571 1.4108 3.5000 0.2706 1.4494 4.0000 0.2033 1.4890 plot(t,y(:,1),t,y(:,2); 实验作业：教材page 88 第5题用ode45命令求解初值问题：实验05 偏微分方程数值解例5-1：用数值法求解下列偏微分方程：写出matlab程序。本微分方程的物理意义是表示：在不考虑流体本身热传导时的套管传热微分方程。例5_1: (第五章课后习题1)例5_2: (第五章课后习题2)具体求解matlab程序，可网上下载。网址：实验06： 有限元法解微分方程clear allsyms x real x0=1x1=1.25x2=1.5x3=1.75x4=2 q00=(x-x1)*(x-x2)/(x0-x1)/(x0-x2)q01=(x-x0)*(x-x2)/(x1-x0)/(x1-x2)q02=(x-x0)*(x-x1)/(x2-x0)/(x2-x1)q10=(x-x3)*(x-x4)/(x2-x3)/(x2-x4)q11=(x-x2)*(x-x4)/(x3-x2)/(x3-x4)q12=(x-x2)*(x-x3)/(x4-x2)/(x4-x3) v0=q00 q01 q02v0v0=v0*v0k0=int(x*v0v0,1,1.5)f0=int(-x*v0,1,1.5)v1=q10 q11 q12v1v1=v1*v1k1=int(x*v1v1,1.5,2)f1=int(-x*v1,1.5,2) dv0=diff(v0,x)dv1=diff(v1,x) dv0v0=dv0*dv0dv1v1=dv1*dv1 dk0=int(x2*dv0v0,1,1.5)dk1=int(x2*dv1v1,1.5,2) kk0=k0+dk0kk1=k1+dk1 kkk0=zeros(5)kkk0(1:3,1:3)=kk0kkk1=zeros(5)kkk1(3:5,3:5)=kk1 k=kkk0+kkk1fff0=zeros(5,1)fff0(1:3)=f0fff1=zeros(5,1)fff1(3:5)=f1 f=fff0+fff1 syms u1d u2d realk_sym=sym(k)f_sym=sym(f)k_sym(1,1)=k_sym(1,1)-(2-u1d)/2k_sym(5,5)=k_sym(5,5)+(4-4*u2d)/4f_sym(1)=f_sym(1)-2f_sym(5)=f_sym(5)+4u01234=inv(k_sym)*f_sym u1d u2d=solve(u01234(1)-2,u01234(5)-4,u1d,u2d) u01234=eval(subs(u01234)中篇 分析应用篇实验07：泡沫材料泡孔（颗粒）大小及形状分析实验目的：了解图像处理一般过程方法技巧及原理熟悉常用图像处理软件的用法及功能了解颗粒大小及形状分析实验原理：粒径分析中的图像处理方法及算法实验软件：/ij/download.htmlwindowsdownload imagej 1.47bundled with with 64-bit java(25mb; requires 64-bit windows),32-bit java(28mb) orwithout java(5mb). (instructions)ij147-jdk6-setup.exe实验素材：自行准备一张有泡孔或颗粒填料的电镜图片。或用绘图软件自行绘制一张。（jpg or bmp） 实验内容：统计测量泡孔或颗粒大小、数量及形状实验步骤：1 二值化adjustauto cesssmoothprocessbinarymake binaryanalyzeanalyze particlesanalyzesummarize13296625525525522725525525538255255255426255255255522552552556402552552557625525525582255255255964255255255101625525525511425525525512325525525513312552552551410255255255153255255255161255255255175255255255181325525525519242552552552032552552552162552552552242552552552382552552552432552552552511255255255262232552552552715255255255281425525525529392552552553011255255255311725525525532525525525533282552552553442552552553522552552553642552552553729255255255381255255255391255255255405525525525541225525525542225525525543172552552554432552552554532552552554639255255255478255255255484255255255499255255255504525525525551125525525552225525525553102552552555452552552555522552552555612225525525557232552552555812552552555912552552556042552552556118255255255621825525525563425