1. 网格

gambit 网格

-----------------
|               |
|               | 6.6 cm
|               | 5 mm
-----------------
     12 cm

20170807/default_id4568.msh

2. fluent 声流的计算

$$\nabla \cdot \bar{v} = 0 \ \rho_0 [\bar{v} \cdot \nabla \bar{v}] = - \rho_0\nabla (\bar{v^\prime v^\prime}) - \nabla \bar{p} + \mu \nabla^2 \bar{v}$$

Rayleigh 应力项，近似为 force source

流程

1. 导入网格 File - Read - Case
2. 网格单位 Grid - Scale - mm, Change Length Units - Scale
3. 定义方程 Define - Models - Solver - Axisymmetric
4. Define - Materials - Fluent Database 找到 water-liquid(h2o<l>) , copy-修改粘性系数($$\mu = 0.016 \sqrt{K}$$) Change/Create
5. Define - Viscous - Laminar (?????)
6. Define - User - Defined Functions - Add - source_term.c, velocity_profile.c - compile - load
7. Define - Boundary Condition - fluid - Set - Source Terms - "Axial Momentum"
8. Define - Boundary Condition - inlet- velocity inlet- set
9. Initialize
10. Solve
11. File - Write

#include "udf.h"
DEFINE_SOURCE(axial_force, C0, t, dS, eqn)
{
    real K = 20.27; real c = 1482; real rho = 998.2;
    real W0 = K*c/rho; 
    real beta = 1/0.03; // length of beam: 1/beta
    real source;

    real x[ND_ND];
    Domain *d = Get_Domain(1); 
    Thread *tw = Lookup_Thread(d,2); // lookup thread of boundary axis; 
    face_t f;

    begin_f_loop(f, tw)
    {
        if(C0 == F_C0(f, tw)) {
            F_CENTROID(x, f, thread);
            source = W0/c * exp (- beta * abs(x[0]));
            break;     // jump out of f_loop;
        } else {
            source = 0;
        }
    }
    end_f_loop(f, tw)

    dS[eqn] = 0;
}

DEFINE_PROFILE(velocity_inlet, thread, position)
{
    real K = 20.27;
    real S = 0.004; // beam width 4mm
    real rho = 998.2;

    real x[ND_ND]; real y;
    face_t f;

    begin_f_loop(f, thread)
    {
        F_CENTROID(x, f, thread);
        y = x[1];
        F_PROFILE(f, thread, position) = sqrt(2*K/pow(rho,2)/M_PI/pow(S,2)) * exp( - y / S);
    }
    end_f_loop(f, thread)
}

3. fluent 声波的计算

$$(\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} ] = \nabla (\bar{v^\prime v^\prime}) - \nabla p^\prime + \mu \nabla^2 v^\prime$$

利用 fluent 的 uds 自定义标量传输方程解声波方程，二维轴对称方程化成标准形式(忽略粘滞， 雷诺应力)

$$\nabla\cdot(\rho_0 v^\prime + \rho \bar{v}) = - \frac{\partial \rho}{\partial t}$$

$$\rho_0 \nabla\cdot(\bar{v}\bar{v}_r + v_r^\prime\bar{v} ) = -\rho_0 \frac{\partial v^\prime_r}{\partial t} - \frac{\partial p^\prime}{\partial r} - \rho_0(v_r^\prime \frac{\partial \bar{v}_r}{\partial r} + v_z^\prime \frac{\partial \bar{v}_r}{\partial z}) $$

$$\rho_0\nabla\cdot(\bar{v}\bar{v}_z + v_z^\prime\bar{v} ) = -\rho_0 \frac{\partial v^\prime_z}{\partial t} - \frac{\partial p^\prime}{\partial z} - \rho_0(v_r^\prime \frac{\partial \bar{v}_z}{\partial r} + v_z^\prime \frac{\partial \bar{v}_z}{\partial z}) $$

结果，

#include "udf.h"

DEFINE_UDS_UNSTEADY(name, c, t, i, apu, su)
{
    real physical_dt, vol, rho;
    rho = 998.2;
    physical_dt = RP_Get_Real("physical-time-step");
    vol = C_VOLUME(c,t);
    if(i == 0) {
        *apu = - vol / physical_dt ;
        *su = vol * C_STORAGE_R(c, t, SV_UDSI_M1(i))/ physical_dt;
    } 

    else if (i == 1) {
       *apu = - rho * vol / physical_dt;
       *su = rho * vol * C_STORAGE_R(c, t, SV_UDSI_M1(i))/ physical_dt;  
    }

    else if (i == 2) {
       *apu = - rho * vol / physical_dt;
       *su = rho * vol * C_STORAGE_R(c, t, SV_UDSI_M1(i))/ physical_dt;
    }
}

#include "udf.h"
DEFINE_UDS_FLUX(name, f, t, i)
{
    Thread *t0, * t1;
    cell_t c0, c1;
    real NV_VEC(flux), NV_VEC(A);

    t0 = THREAD_T0(f,t);
    c0 = F_C0(f,t);
    F_AREA(A, f, t);
    if(i == 0) {
        if(BOUNDARY_FACE_THREAD_P(t)) {
            if (NNULLP(THREAD_STORAGE(t, SV_UDSI(i))
                rho_p = F_UDSI(f, t, 0);
            else
                rho_p = C_UDSI(f, t, 0);
            NV_DS(flux, =, F_U(f,t), F_V(f,t), *, rho_p);
            NV_DS(flux, +=, F_UDSI(f, t, 1), F_UDSI(f, t, 2), * rho);
            return NV_DOT(flux, A);
        } else {
            c1 = F_C1(f,t);
            t1 = F_C1_THREAD(f,t);
            NV_DS(flux, =, C_U(c0,t0), C_V(c0, t0, * C_UDSI(c0, t0, 0));
            NV_DS(flux, +=, C_U(c1,t1), C_V(c1, t1, * C_UDSI(c1, t1, 0));
            NV_DS(flux, +=, C_UDSI(c0, t0, 1), C_UDSI(c0, t0, 2), * rho);
            NV_DS(flux, +=, C_UDSI(c1, t1, 1), C_UDSI(c1, t1, 2), * rho);
            return NV_DOT(flux, A)/ 2.0;
        }
    } 
    else if (i==1) {
        if(BOUNDARY_FACE_THREAD_P(t)) {
            if (NNULLP(THREAD_STORAGE(t, SV_UDSI(i))
                u_p = F_UDSI(f, t, 1);
            else
                u_p = C_UDSI(f, t, 1);
            NV_DS(flux, =, F_U(f,t), F_V(f,t), *, (F_U(f, t)+u_p);
            return NV_DOT(flux, A);
        } else {
            c1 = F_C1(f,t);
            t1 = F_C1_THREAD(f,t);
            NV_DS(flux, =, C_U(c0,t0), C_V(c0,t0), *, (C_U(c0, t0)+C_UDSI(c0, t0, 1));
            NV_DS(flux, +=, C_U(c1,t1), C_V(c1,t1), *, (C_U(c1, t1)+C_UDSI(c1, t1, 1));
            return rho*NV_DOT(flux, A)/2.0;
        }
    } 
    else if (i==2) {
        if(BOUNDARY_FACE_THREAD_P(t)) {
            if (NNULLP(THREAD_STORAGE(t, SV_UDSI(i))
                v_p = F_UDSI(f, t, 2);
            else
                v_p = C_UDSI(f, t, 2);
            NV_DS(flux, =, F_U(f,t), F_V(f,t), *, (F_V(f, t)+v_p);
            return NV_DOT(flux, A);
        } else {
            c1 = F_C1(f,t);
            t1 = F_C1_THREAD(f,t);
            NV_DS(flux, =, C_U(c0,t0), C_V(c0,t0), *, (C_V(c0, t0)+C_UDSI(c0, t0, 2));
            NV_DS(flux, +=, C_U(c1,t1), C_V(c1,t1), *, (C_V(c1, t1)+C_UDSI(c1, t1, 2));
            return rho*NV_DOT(flux, A)/2.0;
        }
    }
}

#include "udf.h"

DEFINE_SOURCE(my_source1, c, t, dS, eqn)
{
    real source;
    real rho = 998.2, c2 = 1483.0 * 1483.0; 
    source = - c2 * C_UDSI_G(c, t, 0)[0] - rho * (C_U_G(c, t)[0]*C_UDSI(c, t, 1) + C_U_G(c, t)[1]*C_UDSI(c, t, 2)); 

    dS[eqn] = 0;
    return source;
}

DEFINE_SOURCE(my_source2, c, t, dS, eqn)
{
    real source;
    real rho = 998.2, c2 = 1483.0 * 1483.0;
    source = - c2 * C_UDSI_G(c, t, 0)[1] - rho * (C_V_G(c, t)[0]*C_UDSI(c, t, 1) + C_V_G(c, t)[1]*C_UDSI(c, t, 2));
    dS[eqn] = 0;
    return source;
}

1. 无扰动声波方程

$$\frac{\partial \rho}{\partial t} + \rho_0 \nabla \cdot v^\prime = 0$$

$$\rho_0 \frac{\partial v^\prime}{\partial t} = - \nabla p$$

#include "udf.h"

DEFINE_UDS_UNSTEADY(my_unsteady, c, t, i, apu, su)
{
    real physical_dt, vol, rho;
    rho = 998.2;
    physical_dt = RP_Get_Real("physcial-time-step");
    vol = C_VOLUME(c,t);
    if(i==0) {
        *apu = - vol / physcial_dt;
        *su = vol * C_STORAGE_R(c, t, SV_UDSI_M1(i)) /physcial_dt;
    } else if (i == 1) {
        *apu = - rho * vol / physical_dt;
        *su = vol * rho * C_STORAGE_R(c, t, SV_UDSI_M1(i)) /physcial_dt;
    } else if (i == 2) {
        *apu = - rho * vol / physical_dt;
        *su = vol * rho * C_STORAGE_R(c, t, SV_UDSI_M1(i)) /physcial_dt;
    }
}

DEFINE_UDS_FLUX(my_flux, f, t, i)
{
    real NV_VEC(V), NV_VEC(A);
    if(i == 0)
        F_AREA(A, f, t);
        NV_D(V, =, F_UDSI(f, t, 1), F_UDSI(f, t, 2));
        return rho * NV_DOT(V, A); 
    else    
        return 0.0;
}

DEFINE_SOURCE(axial_source, c, t, dS, eqn)
{
    real c2 = 1481.0*1481.0;
    dS[eqn] = 0.;
    return = - c2 * C_UDSI_G(c, t, 0)[1];
}

DEFINE_SOURCE(radial_source, c, t, dS, eqn)
{    
    real c2 = 1481.0*1481.0;
    dS[eqn] = 0.;
    return -c2 * C_UDSI_G(c, t, 0)[0];
}