【PY】多项式拟合与样条插值(以及积分)

问题的提出

昨天打算处理叶尖泄漏流动量的,参考GT2017-63655,单位轴向弦长的轴向动量(the tip leakage flow axial momentum per unit axial chord)如下式定义:

image.png

式子里各个参数我都可以通过输出间隙中间平面的参数来得到,然后乘除之后做一个积分。
那么首先我要弄清楚这个平面上的数据结构

面上的数据结构

间隙中间平面的网格,可以认为是一个长的矩形网格,如图:

image.png

然后利用Post中的Export之后,选择这个面,选择要输出的物理量,就得到一个csv文件。

  1. 间隙网格配置
    径向方向上17个点,弦向方向上是44个点,所以一共有17*44=748个数据点。
  2. 数据点的排布
    csv文件中也同样给了748个数据点,那么这些数据点是怎么排布的呢?
    • 是输出相同轴向位置上的17个径向点,然后再移动到下一个轴向位置输出
    • 且先输出尾缘,再逐渐前移
    • 径向方向上则是从tip到casing的顺序,如果用叶高来表达是98%-100%叶高
    • P.S. 有个地方需要注意的是,输出的头34个数据,就是尾缘位置处,两组相近的17个数据重叠在一起,跟后面的规律不一样,所以头34个数据先不用

积分的问题

那么首先我就要对某个轴向位置上的17个径向点进行积分
积分应该怎么做呢,照理说积分其实就是相加,但是有问题:

  • 径向点数少,只有17个
  • 分布很不均匀(近壁面加密)
  • 我还想做不同叶高范围的动量对比,如果按照径向高度来分,中间33%的区域只有三四个点

那么应该怎么做?想了想其实也好解决,而且应该也是更合理的做法:

  • 把我17个离散的数据拟合或插值,得到径向方向上n多数据,然后就可以方便的积分了,结果肯定比我用17个点做出来的要精确的多

多项式拟合

我一开始用的是多项式拟合,用numpy里的polyfit,其实很简单:

import numpy as np
coeffs = np.polyfit(x, y, degree)
p = np.poly1d(coeffs)

说明:

  1. 把离散的x和y数据导入,然后用polyfit一处理,degree是你要用几次多项式来拟合,这样就会得到一个coeffs的list,这是多项式前的系数
  2. 比如f=ax3 + bx2 + cx3 +d,那么coeffts就是[a,b,c,d]
  3. 然后再用poly1d的命令, p就等于ax3 + bx2 + cx3 +d了
  4. 也就是你比如再弄一个xnew = np.arange(0,1,1000)
    那么的ynew=p(xnew),这就是拟合之后0-1范围内的x和y的值了

问题:

  1. 拟合效果不好,于是还是用了样条插值,这样不管能不能用拟合,效果都不错

样条插值

样条插值的做法也很简单:

from scipy import interpolate 
tck = interpolate.splrep(x, y) #需要先使用splrep函数计算出B-Spline曲线的参数
x_new = np.linspace(min(x), max(x), 10000)
y_bspline = interpolate.splev(x_new, tck) #将参数传递给splev函数计算出各个取样点的插值结果

说明:

  1. 就是这样,先import插值的命令
  2. 然后先用splrep计算出b样条曲线的参数tck
  3. 然后再将参数传递给splev函数计算出各取样点的插值结果
    这样得到的结果就很美啦
image.png

这样我在98-100%的间隙范围内就可以去对低、中、高三个范围进行积分了。当然,要先把数据准备好的

对各个轴向位置的径向17个点进行积分,就得到泄漏流动量沿着弦向方向的44个点,如果要积分得到总的动量,就再进行插值和积分的步骤。

积分

  1. trapz
    积分使用的是trapz的命令
    trapz其实是trapezoid的缩写,也就是梯形
    就是梯形法数值积分,也就是一个个矩形块的面积相加吧
    应用起来也非常简单
res = np.trapz(y_bspline,x_new)

在b样条处理中得到的y值和x值,用这个命令已处理,就能得到积分结果了
积分的精确度当然是和x_new的分隔挂钩的,分得区间越多积分越精确。

  1. quad
    此外,还有quad的命令, 这个精度比trapz搞,但是需要精确表达式
    例如:
from scipy import integrate
def half_circle(x): 
  return (1-x**2)**0.5
pi_half, err = integrate.quad(half_circle, -1, 1)

half_circle是圆面积的表达式,后面的两个参数-1和1是积分的范围
quad运行之后得到一个结果和一个误差,所以最后一句里还有个err

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 200,045评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,114评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 147,120评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,902评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,828评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,132评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,590评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,258评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,408评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,335评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,385评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,068评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,660评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,747评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,967评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,406评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,970评论 2 341

推荐阅读更多精彩内容