杭州锐达数字技术有限公司
查看: 8594|回复: 24
打印 上一主题 下一主题

[Python] 用TTM方法生成翼型网格(Python & MATLAB)

? [复制链接]
跳转到指定楼层
楼主
发表于 2010-12-6 02:09 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

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

为了练习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

4.11 KB, 阅读权限: 20, 下载次数: 7

售价: 10 点体能 ?[记录] ?[购买]

所有程序代码(Python版)

TTM方法原理.part1.rar

190 KB, 下载次数: 32

原理摘要(Python版)

TTM方法原理.part2.rar

190 KB, 下载次数: 35

原理摘要(Python版)

TTM方法原理.part3.rar

142.15 KB, 下载次数: 29

原理摘要(Python版)

评分

1

查看全部评分

沙发
?楼主| 发表于 2010-12-6 03:06 | 只看该作者
本帖最后由 Rainyboy 于 2010-12-6 03:08 编辑

MATLAB版的实现更全,一共实现了8种迭代算法,并比较了各种算法的性能。
%方法
nMethod = 7;
% nMethod = 0;%Jacobi点迭代
% nMethod = 1;%Jacobi线迭代
% nMethod = 2; %Jacobi点迭代 + 超松弛 (不可行,计算结果会发散,此处作为对比)
% nMethod = 3;%Jacobi线迭代 + 超松弛??(不可行,计算结果会发散,此处作为对比)
% nMethod = 4;%Guass-Siedel点迭代
% nMethod = 5;%Guass-Siedel线迭代
% nMethod = 6;%Guass-Siedel点迭代 + 超松弛
% nMethod = 7;%Guass-Siedel线迭代 + 超松弛


TTM原理.part1.rar

190 KB, 下载次数: 34

原理及算法对比

TTM原理.part2.rar

190 KB, 下载次数: 25

原理及算法对比

TTM原理.part3.rar

190 KB, 下载次数: 27

原理及算法对比

TTM原理.part4.rar

68.09 KB, 下载次数: 32

原理及算法对比

MatLab代码.rar

8.95 KB, 阅读权限: 20, 下载次数: 8

售价: 20 点体能 ?[记录]

全部MATLAB版本的代码

点评

很不错!学习了? 发表于 2010-12-8 11:24
板凳
发表于 2011-4-21 15:27 | 只看该作者
本帖最后由 woai3181 于 2011-4-21 15:28 编辑

好像不错,但下载不了。谢谢分享
地板
发表于 2011-4-21 15:32 | 只看该作者
能不能让我下载下你的MATLAB版TTM算法,我的权限不够下载。万分感谢!
5
发表于 2011-9-26 23:32 | 只看该作者
值得学习的好贴~~
6
发表于 2012-3-27 12:34 | 只看该作者
学习了~~
头像被屏蔽
7
发表于 2012-9-18 15:13 | 只看该作者
向楼主致敬












? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?
登录/注册后可看大图

济宁:wangzhan678.com??天泰电子:aizk007.com
8
发表于 2012-11-3 13:12 | 只看该作者
9
发表于 2012-12-19 14:47 | 只看该作者
非常有用啊,感谢分享!
10
发表于 2013-5-6 07:51 | 只看该作者
好东西,来看看
11
发表于 2014-3-17 09:07 | 只看该作者
没权限,不能下载怎么办
12
发表于 2014-4-17 21:22 | 只看该作者
这是什么,翼型?
13
发表于 2015-7-11 12:50 | 只看该作者
干货啊,楼主好人
14
发表于 2015-7-11 14:37 | 只看该作者
本帖最后由 eggtor 于 2015-7-11 14:40 编辑

你好,请问能否给我发一下您的TTM法划分翼型外流场网格的matlab代码,我没有权限下载764261405@qq.com
15
发表于 2015-7-12 12:27 | 只看该作者
没权限,不能下载
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表