Debug
时间步长 5e-7，50个步长左右就发散了
- 改成udmi
- source 收敛性
- message
为了试验，
A omega = 40e-6 * 2*pi* 2e+4 = 5.0265
A omega * rho/c = 3.3833
发现 v 分量的边界条件，可能有问题
DEFINE_ADJUST(axis_BC, d)
{
Thread *t;
face_t f;
int zone_ID = 3;
cell_t c;
//Thread *t0 ;
t = Lookup_Thread(d, zone_ID);
if(NULL != THREAD_STORAGE(t,SV_UDS_I(2))
{
begin_f_loop(f,t)
{
c0 = F_C0(f,t);
C_UDSI(c0,t->t0,2) = 0.;
}
end_f_loop(f,t)
}
}
既然轴对称边界上没有 F_UDSI,
Note that F_UDSI is available for wall and flow boundary faces, only. If a UDS attempts to access any other face zone, then an error will result.
那么
DEFINE_PROFILE(axis_boundary, t, i)
{
face_t f0;
Domain *d = Get_Domain(1);
Thread *t0 = Lookup_Thread(d, 3);
begin_f_loop(f0, t0)
{
F_PROFILE(f0, t0, i) = 0.;
}
end_f_loop(f0, t0)
}
1. 标量波动方程（无流场）
$$\frac{\partial \rho}{\partial t} + \rho_0 \nabla \cdot v^\prime = 0$$
$$\rho_0 \frac{\partial v^\prime}{\partial t} = -\nabla p$$
速度分解
$$v^\prime = \nabla \phi + \nabla \times \mathbf{A}$$
$$\begin{cases}\nabla \cdot v^\prime = \zeta\ \nabla \times v^\prime = \vec{\omega} \end{cases}$$
轴对称
$$v^\prime = v_r^\prime\hat{r} + v_z^\prime \hat{z}$$
$$\begin{cases} \nabla\cdot v^\prime = \frac{1}{r} \frac{\partial r v_r^\prime}{\partial r} + \frac{\partial v_z^\prime}{\partial z}\ \nabla\times v^\prime = \frac{\partial v^\prime_r}{\partial z} - \frac{\partial v_z^\prime}{\partial r} \end{cases} $$
波动方程的标量形式
$$\frac{\partial\rho}{\partial t} + \rho_0 \nabla\cdot v^\prime = 0$$
$$\rho_0 \frac{\partial \nabla \cdot v^\prime}{\partial t} = - \nabla^2 p$$
$$\rho_0 \frac{\partial \nabla \times v^\prime}{\partial t} = 0$$
也即
$$\begin{cases} \frac{\partial \rho^\prime}{\partial t} + \rho_0 \zeta = 0 \ \rho_0\frac{\partial \zeta}{\partial t} = -c^2 \nabla^2 \rho^\prime \ \frac{\partial \omega}{\partial t} = 0 \end{cases}$$
边界条件
# include "udf.h"
# include "sg.h"
enum
{
RHO,
ZETA,
N_REQURIED_UDS
}
DEFINE UDS_UNSTEADY(uds_time, c, t, i, apu, su)
{
real physical_dt, vol, rho0 = 998.2, phi_old;
physical_dt = RP_Get_Real("physical-time-step");
vol = C_VOLUME(c,t);
if(i == RHO) {
*apu = -vol / physical_dt;
*su = vol*C_STORAGE_R(c,t,SV_UDSI_M1(i))/physical_dt;
} else if(i == ZETA){
*apu = -rho0*vol / physical_dt;
*su = rho0*vol*C_STORAGE_R(c,t,SV_UDSI_M1(i))/physical_dt;
}
}
DEFINE UDS_FLUX(uds_flux, f, t, i)
{
real k = 1483.0*1483.0;
cell_t c0, c1;
THREAD *t1, *t0;
real A[ND_ND], dRHO[ND_ND], dr0[ND_ND], dr1[ND_ND], es[ND_ND], A_by_ds, ds;
c0 = F_C0(f, t);
t0 = THREAD_T0(t);
if (i != ZETA) return 0.;
if (BOUNDARY_FACE_THREAD_P(t))
{
if(NULL != THREAD_STORAGE(t, SV_UDS_I(RHO)) ) {
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
source = k * (F_UDSI(f, t, RHO) - C_UDSI(c0, t0, RHO)) * A_by_es / ds;
if( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO)) {
BOUNDARY_SECONDARY_GRADIENT_SOURCE(GSource, SV_UDSI_G(RHO), dRHO, es, A_by_es, k);
source += GSource;
}
} else {
source = 0.;
}
} else {
c1 = F_C1(f, t);
t1 = THREAD_T1(t);
INTERIOR(f, t, A, ds, es, A_by_es, dr0, dr1);
source = k * (C_UDSI(c1, t1, rho) - C_UDSI(c0, t0, rho)) * A_by_es / ds;
if( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO)) {
BOUNDARY_SECONDARY_GRADIENT_SOURCE(GSource, SV_UDSI_G(RHO), dRHO, es, A_by_es, k);
source += GSource;
}
}
return source;
}
DEFINE SOURCE(uds0, c, t, dS, eqn)
{
real rho0 = 998.2;
dS[eqn] = 0.;
return - rho0 * C_UDSI(zeta);
}
// Initial Condition
//----- Boundary Condition -------//
DEFINE PROFILE(f, t, i)
{
//real A[ND_ND], dRHO[ND_ND], dr0[ND_ND], dr1[ND_ND], es[ND_ND], A_by_ds, ds;
cell_t c0 = F_C0(f, t);
THREAD *t0 = THREAD_T0;
begin_f_loop(f, t)
{
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
f_profile(f,t,i) = C_UDSI(c0, t0, ZETA)/ ds;
}
end_f_loop(f, t)
}
DEFINE PROFILE(f, t, i)
{
//
begin_f_loop(f, t);
THREAD *t0 = THREAD_T0;
begin_f_loop(f, t)
{
f_profile(f, t, i) = 5.026
}
end_f_loop(f, t)
2. 标量波动方程
$$(\frac{\partial }{\partial t} + \bar{v} \cdot \nabla) \rho + \rho_0 \nabla \cdot v^\prime = 0$$
$$\rho_0 [\frac{\partial v^\prime}{\partial t} + (v^\prime \cdot \nabla ) \bar{v} + (\bar{v} \cdot \nabla) v^\prime + (\bar{v}\cdot\nabla)\bar{v} ] = \rho_0\nabla (\bar{v^\prime v^\prime}) - \nabla p^\prime + \mu \nabla^2 v^\prime $$
所以用速度势
$$\rho_0 {\frac{\partial \nabla \cdot v^\prime}{\partial t} + \nabla\cdot [(v^\prime \cdot \nabla) \bar{v}] + \nabla\cdot[(\bar{v}\cdot\nabla)v^\prime] + \nabla \cdot [(\bar{v}\cdot\nabla)\bar{v}]} = \nabla^2 (\bar{v^\prime v^\prime}) - \nabla^2 p^\prime + \mu \nabla^2 \nabla \cdot v^\prime$$
$$\rho_0 {\frac{\partial \nabla\times v^\prime}{\partial t} + \nabla\times[(v^\prime \cdot \nabla ) \bar{v}] + \nabla\times[ (\bar{v} \cdot \nabla) v^\prime] + \nabla\times[ (\bar{v}\cdot\nabla)\bar{v} ]} = \mu \nabla^2 \nabla\times v^\prime $$
利用
$$\nabla \cdot[(\bar{v}\cdot \nabla)v^\prime] = (\bar{v}\cdot \nabla) (\nabla \cdot v^\prime) + \nabla v^\prime : \nabla\bar{v}$$
$$\nabla\times[(v^\prime \cdot \nabla ) \bar{v}]=(v^\prime \cdot \nabla ) (\nabla \times \bar{v}) + (\nabla \times v^\prime)\cdot \nabla \bar{v}$$
综合
$$\begin{cases}(\frac{\partial}{\partial t} + \bar{v} \cdot \nabla) \rho + \rho_0 \zeta= 0 \ \rho_0 [\frac{\partial \zeta}{\partial t} + \bar{v} \cdot \nabla \zeta + 2 \nabla v^\prime : \nabla \bar{v} + \nabla \bar{v} : \nabla \bar{v}] = - \nabla^2 p^\prime +\mu \nabla^2 \zeta \ \rho_0 [\frac{\partial \omega}{\partial t} + (v^\prime \cdot\nabla )\bar{\omega} + (\bar{v}\cdot \nabla) \omega] = \mu \nabla^2 \omega \end{cases}$$
其中。速度$$v^\prime = \nabla \phi - \nabla \times \vec{\psi}$$ 根据这两个方程，
$$\begin{cases} \nabla^2 \phi = \zeta\ \nabla^2 \psi = \omega\ \end{cases}$$
边界条件
#include "udf.h"
enum {
RHO,
ZETA,
OMEGA,
N_REQUIRED_UDS
}
DEFINE_FLUX(zeta, f, t, i)
{
real k = 1483.0*1483.0;
cell_t c0, c1;
THREAD *t1, *t0;
real A[ND_ND], dRHO[ND_ND], dr0[ND_ND], dr1[ND_ND], es[ND_ND], A_by_ds, ds;
c0 = F_C0(f, t);
t0 = THREAD_T0(t);
if (i != ZETA) return 0.;
if (BOUNDARY_FACE_THREAD_P(t)) {
if(NULL != THREAD_STORAGE(t, SV_UDS_I(RHO)) ) {
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
source = k * (F_UDSI(f, t, RHO) - C_UDSI(c0, t0, RHO)) * A_by_es / ds;
if( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO)) {
BOUNDARY_SECONDARY_GRADIENT_SOURCE(GSource, SV_UDSI_G(RHO), dRHO, es, A_by_es, k);
source += GSource;
}
} else {
source = 0.;
}
if(NULL != THREAD_STORAGE(t, SV_UDS_I(ZETA))
zeta =F_UDSI(f, t, ZETA);
} else if(NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(ZETA)) {
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
zeta =C_UDSI(c0, t0, ZETA) + NV_DOT(C_UDSI_G(c0, t0, ZETA), dr0);
} else {
zeta =C_UDSI(c0, t0, ZETA);
}
NV_DS(V, =, F_U(f,t), F_V(c,t), 0, *, zeta);
return rho0*NV_DOT(V, A) + source;
} else {
c1 = F_C1(f, t);
t1 = THREAD_T1(t);
INTERIOR_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0, dr1);
source = k * (C_UDSI(c1, t1, rho) - C_UDSI(c0, t0, rho)) * A_by_es / ds;
if( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO)) {
SECONDARY_GRADIENT_SOURCE(GSource, SV_UDSI_G(RHO), dRHO, es, A_by_es, k);
source += GSource;
}
if(NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(ZETA)) {
//INTERIOR_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0, dr1);
ZETA = C_UDSI(c0, t0, ZETA) + NV_DOT(C_UDSI_G(c0, t0, ZETA), dr0);
} else {
ZETA = (C_UDSI(c0, t0, ZETA) + C_UDSI(c1, t1, ZETA))/2.0;
}
NV_D(V, =, C_U(c0, t0), C_V(c0, t0), 0);
NV_D(V,+=, C_U(c1, t1), C_V(c1, t1), 0);
return rho0*NV_DOT(V, A)/2.0*zeta + source;
}
}
DEFINE_FLUX(omega, f, t, i)
{
c0 = F_C0(f, t);
t0 = THREAD_T0(t);
if (i != OMEGA) return 0.;
if (BOUNDARY_FACE_THREAD_P(t)) {
if(NULL != THREAD_STORAGE(t, SV_UDS_I(OMEGA)) ) {
omega = F_UDSI(f, t);
} else if ( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(OMEGA) ) {
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
omega = C_UDSI(c0, t0) + NV_DOT(C_UDSI_G(c0, t0, OMEGA), dr0);
} else {
omega = C_UDSI(c0, t0);
}
NV_D(V, =, F_U(f, t), F_V(f, t), 0, omega);
return rho0*NV_DOT(V, A);
} else {
c1 = F_C1(f, t);
t1 = THREAD_T1(t);
if(NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(OMEGA)) {
INTERIOR_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0, dr1);
omega = C_UDSI(c0, t0, OMEGA) + NV_DOT(C_UDSI_G(c0, t0, OMEGA), dr0);
} else {
omega = (C_UDSI(c0, t0, OMEGA) + C_UDSI(c1, t1, OMEGA))/2.0;
}
NV_D(V, =, C_U(c0, t0), C_V(c0, t0), 0);
NV_D(V,+=, C_U(c1, t1), C_V(c1, t1), 0);
return rho0*NV_DOT(V, A)/2.0*omega;
}
}
DEFINE_FLUX(rho, f, t, i)
{
c0 = F_C0(f, t);
t0 = THREAD_T0(t);
if (i != RHO) return 0.;
if (BOUNDARY_FACE_THREAD_P(t)) {
if(NULL != THREAD_STORAGE(t, SV_UDS_I(RHO)) ) {
rho = F_UDSI(f, t);
} else if ( NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO) ) {
BOUNDARY_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0);
rho = C_UDSI(c0, t0) + NV_DOT(C_UDSI_G(c0, t0, RHO), dr0);
} else {
rho = C_UDSI(c0, t0);
}
NV_D(V, =, F_U(f, t), F_V(f, t), 0, RHO);
return rho0*NV_DOT(V, A);
} else {
c1 = F_C1(f, t);
t1 = THREAD_T1(t);
if(NULL != T_STORAGE_R_NV(t0, SV_UDSI_G(RHO)) {
INTERIOR_FACE_GEOMETRY(f, t, A, ds, es, A_by_es, dr0, dr1);
rho = C_UDSI(c0, t0, RHO) + NV_DOT(C_UDSI_G(c0, t0, RHO), dr0);
} else {
rho = (C_UDSI(c0, t0, RHO) + C_UDSI(c1, t1, RHO))/2.0;
}
NV_D(V, =, C_U(c0, t0), C_V(c0, t0), 0);
NV_D(V,+=, C_U(c1, t1), C_V(c1, t1), 0);
return rho0*NV_DOT(V, A)/2.0*rho;
}
}
DEFINE_SOURCE(source_rho, c, t, i)
{
real rho0 = 998.2;
return - rho0*
}