SOD 激波管数值模拟 VL 分裂 weno 5阶精度 以及ROE 格式 python 源码

最新推荐文章于 2024-04-26 16:55:17 发布
最新推荐文章于 2024-04-26 16:55:17 发布
阅读量4.6k 收藏 23
点赞数 5
分类专栏： python3--numpy
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/yhf_naive/article/details/69396674
版权 
#!/usr/bin/python3
################ZJU#####
#@auther: Hongfei Ye  
#date: 2016.4.6
#class: Computational aerodynamics
#project name:One-dimensional SOD shock tube 
#teacher:Shi Xing
#platform : python3.6.1 matplotlib 2.0.0 numpy 1.12.1
########################
from numpy import *
import matplotlib.pyplot as  pl
##########initial#########
delta_t=0.0001
delta_x=0.005#you must choose [0.0001 0.005] because of the sick problem
Nt=int(0.20/delta_t)
Nx=int(1/delta_x)
gama=1.4
P=array([1.0]*int(Nx/2)+[0.1]*int(Nx/2))
F_p=array([[0.0]*Nx,[0.0]*Nx,[0.0]*Nx])
F_n=array([[0.0]*Nx,[0.0]*Nx,[0.0]*Nx])
Den=array([1.0]*int(Nx/2)+[0.125]*int(Nx/2))
u=array([0.0]*Nx)
Den_u=Den*u
FF_p=array([[0.0]*Nx,[0.0]*Nx,[0.0]*Nx])
FF_n=array([[0.0]*Nx,[0.0]*Nx,[0.0]*Nx])
E=P/(gama-1)+0.5*u**2*Den
#########displacement########
def show(s):
    x=[[i/Nx-0.5] for i in range(Nx)]
    pl.figure(1)
    pl.plot(x,P,label="$Pressure$")
    pl.plot(x,u,label="$velocity$")
    pl.plot(x,Den,label="$Density$")
    pl.legend()
    pl.title(s)
    pl.xlabel("Positon x")
    pl.show()
#########compute#######
def VLsplit(level=1):
    global P,u,C,F_p,F_n,Den,E
    for j in range(Nt):
        C=sqrt(abs(gama*P/Den))
        ma=u/C
        lamta=array([u,u-C,u+C])
        for i in range(Nx):
            if(ma[i]>=1):
                F_p[0][i]=Den[i]*u[i]
                F_p[1][i]=u[i]*(Den[i]*u[i]**2+P[i])
                F_p[2][i]=u[i]*(E[i]+P[i])
                F_n[0][i]=0
                F_n[1][i]=0
                F_n[2][i]=0
            if(ma[i]<=-1):
                F_n[0][i]=Den[i]*u[i]
                F_n[1][i]=u[i]*(Den[i]*u[i]**2+P[i])
                F_n[2][i]=u[i]*(E[i]+P[i])
                F_p[0][i]=0
                F_p[1][i]=0
                F_p[2][i]=0
            else:
                fp=Den[i]*C[i]*((ma[i]+1)/2)**2
                fn=-1*Den[i]*C[i]*((ma[i]-1)/2)**2
                F_p[0][i]=fp
                F_p[1][i]=(fp/gama)*((gama-1)*u[i]+2*C[i])
                F_p[2][i]=(fp/(2*gama**2-2))*((gama-1)*u[i]+2*C[i])**2
                F_n[0][i]=fn
                F_n[1][i]=(fn/gama)*((gama-1)*u[i]-2*C[i])
                F_n[2][i]=(fn/(2*gama**2-2))*((gama-1)*u[i]-2*C[i])**2
        if(level==2):##choose the order of accuracy
            for i in range(Nx-4):
                Den[i+2]   -=(0.5*delta_t/delta_x)*(3*F_p[0][i+2]-4*F_p[0][i+1]+F_p[0][i]-F_n[0][i+4]+4*F_n[0][i+3]-3*F_n[0][i+2])
                Den_u[i+2] -=(0.5*delta_t/delta_x)*(3*F_p[1][i+2]-4*F_p[1][i+1]+F_p[1][i]-F_n[1][i+4]+4*F_n[1][i+3]-3*F_n[1][i+2])
                E[i+2]     -=(0.5*delta_t/delta_x)*(3*F_p[2][i+2]-4*F_p[2][i+1]+F_p[2][i]-F_n[2][i+4]+4*F_n[2][i+3]-3*F_n[2][i+2])
        else:
            for i in range(Nx-2):
                Den[i+1]-=   (delta_t/delta_x)*(F_p[0][i+1]+F_n[0][i+2]-F_p[0][i]-F_n[0][i+1])
                Den_u[i+1]-= (delta_t/delta_x)*(F_p[1][i+1]+F_n[1][i+2]-F_p[1][i]-F_n[1][i+1])
                E[i+1]-=     (delta_t/delta_x)*(F_p[2][i+1]+F_n[2][i+2]-F_p[2][i]-F_n[2][i+1])
        u=Den_u/Den
        P=(gama-1)*(E-0.5*Den*u**2)
    show("Van Leer split")
        
def ROE():
    global P,u,C,F_p,F_n,Den,E
    F1=[Den*u,Den*u**2+P,u*(E+P)]
    for j in range(Nt):
        G=[Den,Den_u,E]
        F=[Den_u,Den*u**2+P,u*(E+P)]
        H=(E+P)/Den
        for i in range(Nx-1):
            aveu=(sqrt(Den[i])*u[i]+sqrt(Den[i+1])*u[i+1])/(sqrt(Den[i])+sqrt(Den[i+1]))
            aveDen=((sqrt(Den[i])+sqrt(Den[i+1]))/2)**2
            aveH=(sqrt(Den[i])*H[i]+sqrt(Den[i+1])*H[i+1])/(sqrt(Den[i])+sqrt(Den[i+1]))
            aveP=((gama-1)/gama)*(aveDen*aveH-0.5*aveDen*aveu**2)
            aveE=aveDen*aveH-aveP
            aveDen_u=aveDen*aveu
            avec=((gama-1)*(aveH-0.5*aveu**2))
            aveA=array([[0,                                       1,                                  0],\
                        [((gama-3)/2)*aveu**2,                    (3-gama)*aveu,                      gama-1],\
                        [(gama-2)*aveu**3/2-aveu*avec/(gama-1),avec/(gama-1)+(3-gama)*aveu**2/2,gama*aveu]])
            aveA_1,aveS=linalg.eig(aveA)
            aveA_1=abs(diag(aveA_1))
            aveS_1=linalg.inv(aveS)
            D=dot(aveS,aveA_1)
            K=0.5*dot(aveS_1,array([x[i+1] for x in G])-array([x[i] for x in G]))
            K=dot(D,K)
            for k in range(3):
                F1[k][i+1]=0.5*(F[k][i]+F[k][i+1])-K[k]
            if(i!=0):
                Den[i]-=   (delta_t/delta_x)*(F1[0][i+1]-F1[0][i])
                Den_u[i]-= (delta_t/delta_x)*(F1[1][i+1]-F1[1][i])
                E[i]-=     (delta_t/delta_x)*(F1[2][i+1]-F1[2][i])
        u=Den_u/Den
        P=(gama-1)*(E-0.5*Den*u**2)
    show("ROE")

    
def WENO():
    global P,u,C,F_p,F_n,Den,E,FF_p,FF_n,IS1_p,IS1_n,IS2_p,IS2_n,IS3_p,IS3_n
    for j in range(Nt):
        F=[Den_u,Den*u**2+P,u*(E+P)]
        C=sqrt(abs(gama*P/Den))
        lamta=array([u,u-C,u+C])
        lamta_p=(abs(lamta)+lamta)/2
        lamta_n=lamta-lamta_p
        for i in range(Nx):
            F_p[0][i]=Den[i]/(2*gama)*(2*(gama-1)*lamta_p[0][i]+lamta_p[1][i]+lamta_p[2][i])
            F_p[1][i]=Den[i]/(2*gama)*(2*(gama-1)*lamta_p[0][i]*u[i]+lamta_p[1][i]*(u[i]-C[i])+lamta_p[2][i]*(u[i]+C[i]))
            F_p[2][i]=Den[i]/(2*gama)*((gama-1)*lamta_p[0][i]*u[i]**2+1/2*lamta_p[1][i]*(u[i]-C[i])**2+1/2*lamta_p[2][i]*(u[i]+C[i])**2+(3-gama)/(2*gama-2)*(lamta_p[1][i]+lamta_p[2][i])*(C[i])**2)
            F_n[0][i]=Den[i]/(2*gama)*(2*(gama-1)*lamta_n[0][i]+lamta_n[1][i]+lamta_n[2][i])
            F_n[1][i]=Den[i]/(2*gama)*(2*(gama-1)*lamta_n[0][i]*u[i]+lamta_n[1][i]*(u[i]-C[i])+lamta_n[2][i]*(u[i]+C[i]))
            F_n[2][i]=Den[i]/(2*gama)*((gama-1)*lamta_n[0][i]*u[i]**2+1/2*lamta_n[1][i]*(u[i]-C[i])**2+1/2*lamta_n[2][i]*(u[i]+C[i])**2+(3-gama)/(2*gama-2)*(lamta_n[1][i]+lamta_n[2][i])*(C[i])**2)
        for i in range(Nx-4):
            IS1_p=IS2_p=IS3_p=IS1_n=IS2_n=IS3_n=0
            for k in range(3):
                IS1_p+=1/4*((F_p[k][i]-4*F_p[k][i+1]+3*F_p[k][i+2])**2)+     13/12*((F_p[k][i]  -2*F_p[k][i+1]+F_p[k][i+2])**2)
                IS2_p+=1/4*((F_p[k][i+1]-F_p[k][i+3])**2)+                 13/12*((F_p[k][i+1]-2*F_p[k][i+2]+F_p[k][i+3])**2)
                IS3_p+=1/4*((3*F_p[k][i+2]-4*F_p[k][i+3]+F_p[k][i+4])**2)+   13/12*((F_p[k][i+2]-2*F_p[k][i+3]+F_p[k][i+4])**2)
                IS1_n+=1/4*((F_n[k][i+4]-4*F_n[k][i+3]+3*F_n[k][i+2])**2)+   13/12*((F_n[k][i+4]-2*F_n[k][i+3]+F_n[k][i+2])**2)
                IS2_n+=1/4*((F_n[k][i+3]-F_n[k][i+1])**2)+                 13/12*((F_n[k][i+3]-2*F_n[k][i+2]+F_n[k][i+1])**2)
                IS3_n+=1/4*((3*F_n[k][i+2]-4*F_n[k][i+1]+F_n[k][i])**2)+     13/12*((F_n[k][i+2]-2*F_n[k][i+1]+F_n[k][i])  **2)
            oumiga1_p=(0.1/(IS1_p+1e-6)**2/(0.1/(IS1_p+1e-6)**2+0.6/(IS2_p+1e-6)**2+0.3/(IS3_p+1e-6)**2))#attention,  oumiga is a num instead of an vector
            oumiga2_p=(0.6/(IS2_p+1e-6)**2/(0.1/(IS1_p+1e-6)**2+0.6/(IS2_p+1e-6)**2+0.3/(IS3_p+1e-6)**2))
            oumiga3_p=1-oumiga2_p-oumiga1_p
            oumiga1_n=(0.1/(IS1_n+1e-6)**2/(0.1/(IS1_n+1e-6)**2+0.6/(IS2_n+1e-6)**2+0.3/(IS3_n+1e-6)**2))
            oumiga2_n=(0.6/(IS2_n+1e-6)**2/(0.1/(IS1_n+1e-6)**2+0.6/(IS2_n+1e-6)**2+0.3/(IS3_n+1e-6)**2))
            oumiga3_n=1-oumiga2_n-oumiga1_n
            for k in range(3):
                FF_p[k][i+2]=oumiga1_p*(2*F_p[k][i]/6    -7*F_p[k][i+1]/6   +11*F_p[k][i+2]/6)+oumiga2_p*(-1*F_p[k][i+1]/6 +5*F_p[k][i+2]/6   +2*F_p[k][i+3]/6)+oumiga3_p*(2*F_p[k][i+2]/6  +5*F_p[k][i+3]/6   -1*F_p[k][i+4]/6)
                FF_n[k][i+2]=oumiga1_n*(2*F_n[k][i+4]/6  -7*F_n[k][i+3]/6   +11*F_n[k][i+2]/6)+oumiga2_n*(-1*F_n[k][i+3]/6 +5*F_n[k][i+2]/6   +2*F_n[k][i+1]/6)+oumiga3_n*(2*F_n[k][i+2]/6  +5*F_n[k][i+1]/6     -1*F_n[k][i]/6)
        for i in range(Nx-4):
            if(i>=2):
                Den[i+1]-=   (delta_t/delta_x)*(FF_p[0][i+1]-FF_p[0][i]+FF_n[0][i+2]-FF_n[0][i+1])
                Den_u[i+1]-= (delta_t/delta_x)*(FF_p[1][i+1]-FF_p[1][i]+FF_n[1][i+2]-FF_n[1][i+1])
                E[i+1]-=     (delta_t/delta_x)*(FF_p[2][i+1]-FF_p[2][i]+FF_n[2][i+2]-FF_n[2][i+1])
        u=Den_u/Den
        P=(gama-1)*(E-0.5*Den*u**2)
    show("WENO")
############choose function##########
#VLsplit()
#VLsplit(2)
ROE()
#WENO()

yhfnaive
关注
  • 5
    点赞
  • 踩
  • 23
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
利用5阶精度的WENO方法求解欧拉方程 cfd
10-31
利用5阶精度的WENO方法求解欧拉方程，里面有源代码和说明文档，具体可参考论文Efficient Implementation of Weighted ENO Schemes。-The use of 5-order accuracy of the WENO method for solving Euler equations, which source code and documentation, specific reference papers Efficient Implementation of Weighted ENO Schemes.
一维sod激波管问题weno格式
02-17
一维sod激波管问题weno格式
WENO5_1DB.zip_BURGER_WENO格式_weno matlab_weno方程
07-15
使用5阶WENO格式，3阶RK格式，LF通量分裂，对一维Burger s方程进行求解。所得结果可以与《Weighted Essentially Non-Oscillatory Schemes》1994中结果进行对比
minmod.rar_TVD matlab_TVD 格式_TVD 计算_minmod函数_计算流体力学
09-14
minmod限制器函数，计算流体力学TVD格式中的开关函数。
计算流体力学系列（二）：针对如下Sod 激波管问题（Matlab版）
Nidebear的博客
03-27 9268
原创性声明：本文属于作者原创性文章，属于课程作业内容，转载请注明出处，仅为同学提供学习交流之用。如需程序请＋企鹅：2659692290。
一维Sod激波管问题Riemann间断解算例及python代码
03-30
对一维Sod激波管问题的分析和数值求解，内附python代码及注释。本算例中初值已经给定，首先通过理论推导和代入数值分析确定了算例中属于左膨胀波右侧激波情形，而后根据已知信息进行编程，用二分法求解未知压强，进而计算所有参数并绘制图像。
五阶weno格式的m文件
06-07
五阶精度weno格式matlab代码 时间推进三阶rk格式 求解激波稀疏波问题
Sod激波管 五阶WENO 流通矢量分裂 三阶RK
最新发布
tokisaygg的博客
04-26 549
以λ+ =（（λ² + epsl²）^0 λ + .5）/2 ，λ- = （（λ² - epsl²）^0 λ + .5）/2 将三个特征值进行分裂，求h j+1/2 在（a，b，c，d，e）下的差分格式，联立三个三阶差分模板（a，b，c），（b，c，d），（c，d，e），求解出权值c1，c2，c3。4.使用特征值以及对应矩阵求出 F+=。1.构造模板 a=f j-2 ，b=f j-1 ，c= f j ，d= f j+1 ， e=f j+2。2.求出特征值 λ1=u ，λ2=u+c ，λ3=u-c。
matlab-流体Sod激波管问题L-W格式、Roe格式、Van Leer格式（迎风格式）、5阶WENO格式计算其数值解
a_zxswer的博客
07-04 6963
目录Sod 激波管问题1、LW：2、Roe3、Van Leer4、五阶weno源码+文档+演示视频 Sod 激波管问题 针对如下Sod 激波管问题： 分别用L-W格式（搞定）、Roe格式（搞定）、Van Leer格式（迎风格式）、5阶WENO格式（）计算其数值解。 要求：1） 简述求解原则，注明时间推进格式，时间步长； 2） 列出全部计算公式和步骤； 3） 使用三套空间网格数，分别是100,200,400； 4） 画出t=0.14时刻密度、速度及压力的分布，并与精确解进行比较（要求同一套网格四种格式的数
一维激波管问题——计算流体力学（Roe格式）
03-31
本文用python3.6对一维激波管lax激波管问题进行了计算，使用了Roe一阶迎风格式，激波管问题的初值采用了lax激波管的初值（具体见代码），推进时间为0.14s，网格数为1001，Roe 格式在激波处与解析解符合的较好, 而在接触间断处存在有耗散。
weno格式源代码
10-06
weno代码
sod_激波管精确解_激波管_精确解_一维python_一维激波管_
09-30
用python写的一维激波管精确解，作为参考
加权基本非振荡 (WENO) 方案：使用 5 阶 WENO 方案求解非线性双曲方程。-matlab开发
05-30
C.-W. 审查的五阶 WENO 方案的一维实现。 Shu 在“对流主导问题的高阶加权基本非振荡方案”，SIAM 评论，51：82-126，（2009 年）。 本代码旨在作为该方法实现的指南。 因此，当前代码背后的哲学是可读而不是高效。 医学博士 (2015)
sod_SOD_Sod激波管_RIEMANN_激波管_
10-04
利用黎曼解计算sod激波管问题，输出t=0.14时的压强速度密度。
一维激波管（Lax shock tube）问题的数值求解
热门推荐
LSEC小陆的博客
05-28 1万+
一维激波管（Lax shock tube）问题的数值求解 文章目录一维激波管（Lax shock tube）问题的数值求解@[toc] 问题描述一般的守恒算法格式精确解常用守恒格式（先介绍标量）Lax-Friedrichs格式Roe格式Lax-Wendroff格式ENO格式简介（先介绍标量）5阶WENO格式空间离散Lax-Friedrichs通量分裂原函数方法重构左值右值空间算子和时间三阶Rung...
Sod激波管 五阶WENO格式+流通矢量分裂+三阶RK推进 （Python）
04-26
Sod激波管 五阶WENO格式+流通矢量分裂+三阶RK推进 （Python）
五阶精度weno格式代码.zip_5阶WENO格式_5阶weno_partlygmd_五阶weno格式的m文件_激波
09-23
五阶精度的有限区线性重构Weighted Essentially Non-Oscillatory（WENO）格式是一种高级的数值方法，常用于解决偏微分方程，特别是流体力学中的守恒律方程组，如Euler方程或Navier-Stokes方程。该方法在保持高分辨率...
SOD激波管代码_一维SOD激波管求解_源码
10-01
总的来说，Python求解一维SOD激波管问题结合ROE平均法，提供了一种直观且有效的数值模拟手段，能够帮助理解和验证各种数值方法在处理激波问题时的行为。通过阅读和理解这段代码，读者不仅可以深入学习流体动力学的...
sod激波管中的sod
01-10
SOD是指"Superoxide Dismutase"，中文称为超氧化物歧化酶。 超氧化物歧化酶是一种重要的抗氧化酶，存在于生物体内。它能够催化超氧阴离子（O2-）的还原反应，将其转化为较为稳定的氧气（O2）和过氧化氢（H2O2）。这个反应对于细胞内氧化还原平衡的维持以及减少由活性氧自由基引起的细胞损伤具有重要作用。 SOD激波管中的"SOD"是对超氧化物歧化酶的简称。激波管是一种常用于研究气体动力学和激波传递等问题的实验装置。在SOD激波管中，我们通常把超氧化物歧化酶作为一种反应催化剂，并将氧气与超氧类为初始气体。当激发源激发通过激波管时，气体会高速逃逸，而在反应催化剂的作用下，超氧类与超氧化物歧化酶发生反应生成较为稳定的氧气和过氧化氢物质。这一反应在激波管中的过程被称为"SOD反应"。 通过在实验装置中使用SOD激波管，我们能够研究超氧自由基的特性以及反应动力学，从而更好地了解细胞内氧化还原平衡的调控机制，以及超氧化物歧化酶在细胞保护中的重要作用。这对于进一步研究氧化应激以及相关疾病的发生机制有着重要的意义。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值