用TTM方法生成翼型网格（Python & MATLAB）

2016-4-5 14:46| 发布者: aspen| 查看: 3859| 评论: 24|原作者: Rainyboy

摘要: 为了练习Python，把以前拿MATLAB写过的这个算法重写了一遍，算是对Python在数值计算有一些体会了吧。原理见附件，程序入口点是CFDMesh.py，其实整个大流程就是一个迭代而已。代码贴在下面，Python语义中很重要的一 ...

为了练习Python，把以前拿MATLAB写过的这个算法重写了一遍，算是对Python在数值计算有一些体会了吧。原理见附件，程序入口点是CFDMesh.py，其实整个大流程就是一个迭代而已。
代码贴在下面，Python语义中很重要的一部分就是缩进……如果用“代码”功能的话就贴乱了，我还是用分割线吧。
附件里也有每个文件的压缩包。

第一次用Python写这种规模的数值程序，肯定有很多地方的实现过于罗嗦，还请大家不吝指正！

============代码文件分割线：CFDMesh.py=========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
import numpy as np;import matplotlib.pyplot as plt
from fCreatInitMesh import *;
from fPlotMesh import *;
from fAdjustBound import *;
from fPoint_Itera import *;
from fCheckConvergence import *;
if __name__ == '__main__':
#-----------------------------
#参数表
#-----------------------------
#根据结构的对称性，只储存和计算一半的网格
#因此，周向分网数实际上是实际值的一半
IM = 18;#实际周向分网数=IM*2
JM = 21;#径向分网数
L_C = 1;#弦长
R_OUT = 2.5*L_C;#网格外径
#生成初始网格
(x,y)=CreatInitMesh(R_OUT,JM,IM);
#画出初始网格
PlotMesh(x,y,IM,JM,1);
#准备迭代：准备储存当前步和推进步的空间
y_next = np.zeros((JM,IM));
x_next = np.zeros((JM,IM));
y_curr = np.zeros((JM,IM));
x_curr = np.zeros((JM,IM));
y_curr[:,:] = y[:,:];
x_curr[:,:] = x[:,:];
#设置边界条件
x_next[0,:] = x_curr[0,:];
x_next[JM-1,:] = x_curr[JM-1,:];
y_next[0,:] = y_curr[0,:];
y_next[JM-1,:] = y_curr[JM-1,:];
x_next[:,IM-1] = x_curr[:,IM-1];
y_next[:,0] = y_curr[:,0];
#x_next[:,0] = x_curr[:,0];#这个边界条件需要动态修正,在AdjustBound中进行
y_next[:,IM-1] = y_curr[:,IM-1];
#储存误差历程的向量
MaxError = np.zeros((1,5000));
re = 1;
for n in range(0,500):
      #进一步调整边界条件
      AdjustBound(x_curr,y_curr,x_next,y_next,IM,JM);
      #迭代一步
      Point_Itera(x_curr,y_curr,x_next,y_next,IM,JM);
      #检查收敛
      (re , MaxError[0][n]) = CheckConvergence(x_curr,y_curr,x_next,y_next,IM,JM);
      if re == 1:
         break;
      #递推
      x_curr[:,:] = x_next[:,:];
      y_curr[:,:] = y_next[:,:];
#画出最终网格
PlotMesh(x_next,y_next,IM,JM,2);
#画出残差随迭代的变化
fig = plt.figure(3);
ax = fig.add_subplot(1,1,1);
ax.plot(MaxError[0][0:n]);
plt.show();
============代码文件分割线：================================

============代码文件分割线：fCreatInitMesh.py========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
#子模块，用于生成初始网格
import numpy as np
from fFZeros import *
def CreatInitMesh(R_OUT,JM,IM):
#翼型参数 y = p0*x^0.5 + p1*x + p2*x^2 + p3*x^3 + p4*x^4
p0 = 0.2969;
p1 = -0.1260;
p2 = -0.3516;
p3 = 0.2843;
p4 = -0.1015;
#直角坐标表示的点
x = np.zeros((JM,IM));
y = np.zeros((JM,IM));
#极坐标表示的点
r = np.zeros((JM,IM));
thita = np.linspace(0,np.pi,IM);
#得到最内层的点
for i in range(0,IM-1):
      k1 = np.cos(thita);
      k2 = np.sin(thita);
      try:
         r[0,i] = FZeros(y_s,(-0.5,10,0.5),(p0,p1,p2,p3,p4,k1,k2));
      except NoRootException ,e:
         print i;
         print "在所设置的区间" , e.ebound , "及参数表" , e.eparalist , "内没有根" ;
r[0,IM-1] = 0.5;
#得到最外层的点
r[JM-1,:] = R_OUT;
#插值得到中间的点
for i in range(1,JM-1):
      r[i,:] = (JM-i-1)*1.0/(JM-1)*r[0,:] + (-i)*1.0/(1-JM)*r[JM-1,:];
#极坐标转换为直角坐标
for i in range(0,JM):
      x[i,:] = r[i,:] * np.cos(thita) + 0.5;
      y[i,:] = r[i,:] * np.sin(thita);
return (x,y);

def y_s((r,p0,p1,p2,p3,p4,k1,k2)):
return p0*(k1*r+0.5)**0.5 + p1*(k1*r+0.5) + p2*(k1*r+0.5)**2 + p3*(k1*r+0.5)**3 + p4*(k1*r+0.5)**4 - k2*r;

if __name__ == '__main__':
try:
      IM = 18;#实际周向分网数=IM*2
      JM = 21;#径向分网数
      L_C = 1;#弦长
      R_OUT = 2.5*L_C;#网格外径

      CreatInitMesh(R_OUT,JM,IM);
except NoRootException ,e:
      print "在所设置的区间" , e.ebound , "及参数表" , e.eparalist , "内没有根" ;
============代码文件分割线：================================

============代码文件分割线：fPlotMesh.py========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
#子模块，画出网格
import matplotlib.pyplot as plt
import numpy as np
def PlotMesh(x,y,IM,JM,n):
fig = plt.figure(n);
ax = fig.add_subplot(1,1,1);
for i in range(0,JM):
      ax.plot(x,y,color='blue');
      ax.plot(x,-y,color='blue');
for i in range(0,IM):
      ax.plot(x[:,i],y[:,i],color='blue');
      ax.plot(x[:,i],-y[:,i],color='blue');
plt.show();
============代码文件分割线：================================

============代码文件分割线：fZeros.py========================
# -*- coding: cp936 -*-
#牛顿法求一元函数在制定初值附近的根
#范雨 2010/12/5
def FZeros(func,bound,paralist):
"""对指定的一元函数在指定的点求导数值

      func    函数对象，其定义形式应为，（注意，该函数以一个元组为参数）
               test((x,p1,p2,p3,p4,...))
      bound 求根区间及初始值(x_min,x_max,x_init)
      paralist (p1,p2,p3,p4,...)
"""
(x_min,x_max,x_init) = bound;
root_last = x_init;
root_now = x_init;
while True:
      root_now = root_last - func((root_last,) + paralist)/diff(func,root_last,paralist);
      if abs(root_last - root_now) <= 1e-4:
         break;
      if root_now >= x_min and root_now <= x_max:
         root_last = root_now;
      else:
         raise NoRootException(bound,paralist);
return root_now;

def diff(func,x,paralist):
"""对指定的一元函数在指定的点求导数值

      func    函数对象，其定义形式应为，（注意，该函数以一个元组为参数）
               test((x,p1,p2,p3,p4,...))
      x       自变量的值
      paralist (p1,p2,p3,p4,...)
"""
dx = 1e-8;
return (func((x + dx,) + paralist) - func((x - dx,) + paralist))/(2*dx);

class NoRootException(Exception):
"""求根过程自定义的异常"""
def __init__(Self,bound,paralist):
      Exception.__init__(Self);
      Self.ebound = bound;
      Self.eparalist = paralist;

if __name__ == '__main__':
def test((x,p1,p2,p3)):
      y = p1*x**2 + p2*x + p3;
      return y;
try:
      print FZeros(test,(5,10,6),(1,-4,3));
except NoRootException ,e:
      print "在所设置的区间" , e.ebound , "及参数表" , e.eparalist , "内没有根" ;
============代码文件分割线：================================

============代码文件分割线：fAdjustBound.py========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
#子函数用于动态调整边界条件
def AdjustBound(x_curr,y_curr,x_next,y_next,IM,JM):
i=0;
for j in range(1,JM-1):
      pxpI = (x_curr[j,i+1]-x_curr[j,i+1])/2;
      pxpJ = (x_curr[j+1,i]-x_curr[j-1,i])/2;
      pypI = (y_curr[j,i+1]+y_curr[j,i+1])/2;
      pypJ = (y_curr[j+1,i]-y_curr[j-1,i])/2;
      a = pxpI**2 + pypI**2;
      b = pxpI*pxpJ + pypI*pypJ;
      c = pxpJ**2 + pypJ**2;
      x_next[j,i]= (a*(x_curr[j+1,i]+x_curr[j-1,i])-b*(x_curr[j+1,i+1]-x_curr[j-1,i+1]-x_curr[j+1,i+1]+x_curr[j-1,i+1])/2 + c*(x_curr[j,i+1]+x_curr[j,i+1]))/(a+c)/2;
i=IM-1;
for j in range(1,JM-1):
      pxpI = (x_curr[j,i-1]-x_curr[j,i-1])/2;
      pxpJ = (x_curr[j+1,i]-x_curr[j-1,i])/2;
      pypI = (y_curr[j,i-1]+y_curr[j,i-1])/2;
      pypJ = (y_curr[j+1,i]-y_curr[j-1,i])/2;
      a = pxpI**2 + pypI**2;
      b = pxpI*pxpJ + pypI*pypJ;
      c = pxpJ**2 + pypJ**2;
      x_next[j,i]= (a*(x_curr[j+1,i]+x_curr[j-1,i])-b*(x_curr[j+1,i-1]-x_curr[j-1,i-1]-x_curr[j+1,i-1]+x_curr[j-1,i-1])/2 + c*(x_curr[j,i-1]+x_curr[j,i-1]))/(a+c)/2;
============代码文件分割线：================================

============代码文件分割线：fCheckConvergence.py========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
#子函数，检查收敛
import numpy as np;
def CheckConvergence(x_curr,y_curr,x_next,y_next,IM,JM):
lim_max_error = 1e-5;
max_error = np.max(np.abs(x_curr[:,:]-x_next[:,:])) + np.max(np.abs(y_curr[:,:]-y_next[:,:]));
re = 1;
if max_error < lim_max_error:
      re = 1;
else:
      re = 0;
print max_error;
return (re , max_error);
============代码文件分割线：================================

============代码文件分割线：fPoint_Itera.py========================
# -*- coding: cp936 -*-
#用TTM方法生车贴体网格
#范雨 2010/12/5
#子函数，用于迭代一步
def Point_Itera(x_curr,y_curr,x_next,y_next,IM,JM):
for i in range(1,IM-1):
      for j in range(1,JM-1):
         #得到四个偏导数
         pxpI = (x_curr[j,i+1]-x_curr[j,i-1])/2;
         pxpJ = (x_curr[j+1,i]-x_curr[j-1,i])/2;
         pypI = (y_curr[j,i+1]-y_curr[j,i-1])/2;
         pypJ = (y_curr[j+1,i]-y_curr[j-1,i])/2;
         a = pxpI**2 + pypI**2;
         b = pxpI*pxpJ + pypI*pypJ;
         c = pxpJ**2 + pypJ**2;
         x_next[j,i]= (a*(x_curr[j+1,i]+x_curr[j-1,i])-b*(x_curr[j+1,i+1]-x_curr[j-1,i+1]-x_curr[j+1,i-1]+x_curr[j-1,i-1])/2 + c*(x_curr[j,i+1]+x_curr[j,i-1]))/(a+c)/2;
         y_next[j,i]= (a*(y_curr[j+1,i]+y_curr[j-1,i])-b*(y_curr[j+1,i+1]-y_curr[j-1,i+1]-y_curr[j+1,i-1]+y_curr[j-1,i-1])/2 + c*(y_curr[j,i+1]+y_curr[j,i-1]))/(a+c)/2;
============代码文件分割线：================================

初始网格：

最终网格：

误差随迭代变化：

CFDMesh.rar

TTM方法原理.part1.rar

TTM方法原理.part2.rar

TTM方法原理.part3.rar

本帖最后由 Rainyboy 于 2010-12-6 03:08 编辑

12 / 2 页下一页

收藏分享邀请

上一篇：Simply Fortran编译器下一篇：数值模拟偏微分方程的三种方法介绍

发表评论

账号		自动登录	找回密码
密码			我要加入

用TTM方法生成翼型网格（Python & MATLAB）

最新评论

相关分类