大家好，又见面了，我是你们的朋友全栈君。
帮别人实现rubber friction的运算公式。公式不是我专业的内容，完全不懂，只是照着公式用python计算并画图。做出的图，与论文中的图进行对比，完全一致。
之所以用python，是因为matlab实在太大了。。。
python科学计算需要用到numpy和scipy，画图需要用到matplotlib。这三个模块的安装，windows下很容易，mac下稍微麻烦一些。
matplotlib需要注意的是，因为公式涉及的三个图，分别是对数-对数，对数-常数，常数-对数坐标，所以不能仅仅使用plot。
步骤一：
公式：![Python画图[通俗易懂]](http://images0.cnblogs.com/blog/120731/201302/28131821-089306e0b3194537904f5e386d0cb615.png)
![Python画图[通俗易懂]](https://javaforall.net/wp-content/uploads/2020/11/2020110817443450.jpg)
代码：
![Python画图[通俗易懂]](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
![Python画图[通俗易懂]](https://images.cnblogs.com/OutliningIndicators/ExpandedBlockStart.gif)
View Code
1 def EOmega(E1,omega,tau,a1): 2 a=complex(1,-omega*tau) 3 b=complex(a1+1,-omega*tau) 4 return E1*a/b 5 6 def DrawEOmega(E1,tau,a1,xmin,xmax): 7 re=[] 8 im=[] 9 10 omega = np.logspace(xmin, xmax, 1000) 11 for o in omega: 12 re.append(EOmega(E1,o,tau,a1).real) 13 im.append(abs(EOmega(E1,o,tau,a1).imag)) 14 15 plt.figure() 16 plt.loglog(omega,re,label="$Re(E)$",color="red",linewidth=2) 17 plt.loglog(omega,im,"b--",label="$Im(E)$") 18 plt.ylim(10**0,10**10) 19 plt.xlabel("frequency omega[rad/s]") 20 plt.ylabel("E-modulus[Pa]") 21 plt.title("Complex E-modulus for rheological model") 22 plt.legend(loc="center right") 23 plt.show() 24 25 26 E1=10.0**9 27 a1=1000.0 28 tau=1.0/10**3 29 xmin=-2 30 xmax=8 31 DrawEOmega(E1,tau,a1,xmin,xmax)
图：
![Python画图[通俗易懂]](http://images0.cnblogs.com/blog/120731/201302/28132235-fb62bfea62454707983a6f56521bccd6.png)
步骤二：
公式：![Python画图[通俗易懂]](http://images0.cnblogs.com/blog/120731/201302/28132504-230a167836444d3da32c7bedbf978cf5.png)
（注：公式有些错误，G(q)的第二个积分没有Im）
代码：
View Code
1 def EOmegaExpand(E1,tau,a1,q0,zeta,v,phi): 2 omega=q0*zeta*v*np.cos(phi) 3 return EOmega(E1,omega,tau,a1) 4 5 def GQPartTwo(phi,E1,tau,a1,q0,zeta,v,gamma,sigma0): 6 a=EOmegaExpand(E1,tau,a1,q0,zeta,v,phi) 7 b=(1.0-gamma**2)*sigma0 8 c=a/b 9 re=c.real 10 im=c.imag 11 d=re**2+im**2 12 return d 13 14 def GQPartOne(zeta,E1,tau,a1,q0,v,gamma,sigma0,H): 15 a=quad(GQPartTwo, 0, 2*np.pi,args=(E1,tau,a1,q0,zeta,v,gamma,sigma0))[0] 16 b=(zeta**(1.0-2*H))*a 17 return b 18 19 def GQ(E1,tau,a1,q0,v,sigma0,gamma,H,h0,q): 20 a=quad(GQPartOne,1,q/q0,args=(E1,tau,a1,q0,v,gamma,sigma0,H))[0] 21 b=((q0*h0)**2)*H 22 c=a*b/(16*np.pi) 23 return c 24 25 def PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v,zeta): 26 q=q0*zeta 27 a=GQ(E1,tau,a1,q0,v,sigma0,gamma,H,h0,q) 28 b=1/(2*np.sqrt(a)) 29 c=erf(b) 30 return c 31 32 def DrawPQ(E1,tau,a1,q0,sigma0,gamma,H,h0,xmin,xmax): 33 p1=[] 34 p2=[] 35 p3=[] 36 p4=[] 37 p5=[] 38 39 v1=0.00005 40 v2=0.013 41 v3=0.316 42 v4=7.85 43 v5=187 44 45 l = np.linspace(xmin, xmax, 50) 46 for o in l: 47 p1.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v1,o)) 48 p2.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v2,o)) 49 p3.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v3,o)) 50 p4.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v4,o)) 51 p5.append(PQ(E1,tau,a1,q0,sigma0,gamma,H,h0,v5,o)) 52 53 plt.figure() 54 plt.semilogy(l,p1,"b--",label="$v=0.00005m/s$",linewidth=5) 55 plt.semilogy(l,p2,label="$v=0.013m/s$",linewidth=5) 56 plt.semilogy(l,p3,"r--",label="$v=0.316m/s$",linewidth=1) 57 plt.semilogy(l,p4,"b--",label="$v=7.85m/s$",linewidth=1) 58 plt.semilogy(l,p5,label="$v=187m/s$",linewidth=1) 59 plt.xlabel("zeta[-]") 60 plt.ylabel("P[-]") 61 plt.title("Apparent contact area") 62 plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) 63 plt.show() 64 65 E1=10.0**9 66 a1=1000.0 67 tau=1.0/10**3 68 q0=2000.0 69 sigma0=0.2*(10**6) 70 gamma=0.5 71 H=0.85 72 h0=5.0/10**4 73 xmin=1 74 xmax=100 75 DrawPQ(E1,tau,a1,q0,sigma0,gamma,H,h0,xmin,xmax)
图：
![Python画图[通俗易懂]](http://images0.cnblogs.com/blog/120731/201302/28132813-c59b9e733a31465bbe81a15b912103ef.png)
步骤三：
公式：![Python画图[通俗易懂]](http://images0.cnblogs.com/blog/120731/201302/28133021-baabee6fa93a4aabbdffd8492bc150a7.png)
这一步的代码没有什么难度，只是注意一下坐标的问题即可。
目前还有一点问题，即第二张图中的图例不知道为什么总是无法显示完全，需要再查一下官方文档。
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/155983.html原文链接：https://javaforall.net
