1
$\begingroup$

Problem

I am using finite difference method to solve classic problem of flow around cylinder, for validation of my group's immersive boundary method.

The common way to validate numerical method is calculating drag and lift coefficient. But my result is far form other's results. So it must be error in my force calculation, and I don't know what is wrong.

Two resolution are tested:

  1. Nx = Ny = 512, H = 1/64, dt = 1.8e-3, Nib = 100

  2. Nx = Ny = 1024, H = 1/128, dt = 0.9e-3, Nib = 200

where Nx, Ny is grid size of staggered grid, it means that the pressure field's size if Nx * Ny, u component of velocity field is (Nx + 1) * Ny, v component of velocity field is Nx * (Ny + 1), Nib is the IB points number in cylinder virtual surface.

I have two method to calculate force, the result is shown by coefficient. The first method result:

enter image description here

enter image description here

the second method:

enter image description here

enter image description here

As you can see, almost all results have violent oscillations. I don't know why.

The correct result from

Lai, Ming-Chih, and Charles S. Peskin. "An immersed boundary method with formal second-order accuracy and reduced numerical viscosity." Journal of computational Physics 160.2 (2000): 705-719.

is

enter image description here

So you can see my result have grid independence, but method 1 result is several times the normal coefficient, method 2 have something wrong since its lift coefficient' s origin isn't 0.

I will show my solving pipeline below, thanks for any advice!

Drag and lift coefficient

Formula of lift coefficient:

$$ C_l = \dfrac{F_L}{qS} = \dfrac{F_L}{\frac{1}{2}\rho u^2 S} = \dfrac{2F_L}{\rho u^2 S}, $$

where $F_L$ is lift force, $\rho$ is fluid density, $u$ is inlet velocity, $S$ is reference surface.

Formula of drag coefficient:

$$ C_d = \dfrac{2F_D}{\rho u^2 S}, $$

where $F_D$ is lift force.

Calculating force

All my test here is the same setting: fluid domain size = [0, 8], cylinder center position = (1.85, 4), cylinder diameter D = 0.3, Re = 100, fluid density $\rho$ = 1, fluid inlet velocity $U$= 1, fluid dynamic viscosity $\mu$ is calculate according to Re.

Two resolution are tested:

  1. Nx = Ny = 512, H = 1/64, dt = 1.8e-3, Nib = 100

  2. Nx = Ny = 1024, H = 1/128, dt = 0.9e-3, Nib = 200

In code, it is configed by definition.

#define H          1.0 / 128.0
#define Nx         8.0 / (H)
#define Ny         8.0 / (H)
#define T_STEP     0.9e-3
#define MAX_STEP   100000
#define OUTPUT_GAP 500

#define CENTER_X 1.85
#define CENTER_Y 4
#define CY_R     0.15
#define CY_N     200

#define INLET_FLUID_VELOCITY 1.0
#define FLUID_DENSITY        1.0

#define REYNOLDS_NUM 100.0

#define DYNAMIC_VISCOSITY   FLUID_DENSITY* INLET_FLUID_VELOCITY* CY_R * 2.0 / REYNOLDS_NUM
#define KINEMATIC_VISCOSITY DYNAMIC_VISCOSITY / FLUID_DENSITY

To know the drag and lift force equals to calculating force of solid applied by fluid.

My methods are obtained from this paper:

Lai, Ming-Chih, and Charles S. Peskin. "An immersed boundary method with formal second-order accuracy and reduced numerical viscosity." Journal of computational Physics 160.2 (2000): 705-719.

I have tested two methods in this paper: One is directly sum up the IB force in all IB points, the other is using momentum conservation of NS equation.

The first method

$$ F_D = -\int_{\Omega}f_1 \mathrm{d} \mathbf{x} = -\int_0^{L_b}F_1 \mathrm{d}s, $$

$$ F_L = -\int_{\Omega}f_2 \mathrm{d} \mathbf{x} = -\int_0^{L_b}F_2 \mathrm{d}s, $$

where $F_D,F_L$ is drag force and lift force, $f_1, f_2$ is x and y component of IB force distributed to Euler grid, $F_1, F_2$ is x and y component of IB force at Lagrangian IB points. $\mathrm{d}\mathbf{x}$ seems like position of Euler grid, but it is annotation of the paper. I may think $\mathrm{d}\sigma$ looks nice, because it implies delta surface.

In code, to get the force, I sum up the IB force and then multiply it with h^2, h is interval length of IB points.

The second method:

NS equation in x component

$$ \dfrac{\partial u_1}{\partial t} + (\mathbf u \cdot \nabla)u_1 = -\dfrac{1}{\rho} \dfrac{\partial p}{\partial x} + \dfrac{\mu}{\rho}\nabla \cdot \nabla u_1 + f_1 $$

Integration form

$$ \dfrac{\partial}{\partial t}\int_{\Omega_0}\rho u_1 \mathrm{d} \mathbf x + \int_{\partial \Omega_0} \rho u_1 \mathbf u \cdot \mathbf n \mathrm{d} s = - \int_{\Omega_0} p n_1 \mathrm{d} \mathbf x + \int_{\partial \Omega_0} \mu (\dfrac{\partial u_1}{\partial x_j} + \dfrac{\partial u_j}{\partial x_1}) n_j \mathrm{d} s + \int_{\Omega_0} f_1 \mathrm{d} \mathbf x $$

When fluid become steady, the first term = 0, then

$$ \int_{\partial \Omega_0} \rho u_1 \mathbf u \cdot \mathbf n \mathrm{d} s = - \int_{\Omega_0} p n_1 \mathrm{d} \mathbf x + \int_{\partial \Omega_0} \mu (\dfrac{\partial u_1}{\partial x_j} + \dfrac{\partial u_j}{\partial x_1}) n_j \mathrm{d} s + \int_{\Omega_0} f_1 \mathrm{d} \mathbf x $$

Then drag force

$$ F_D = -\int_{\Omega_0} f_1 \mathrm{d} \mathbf x = -\int_{\partial \Omega_0} \rho u_1 \mathbf u \cdot \mathbf n \mathrm{d} s - \int_{\partial \Omega_0} p n_1 \mathrm{d} s + \int_{\partial \Omega_0} \mu (\dfrac{\partial u_1}{\partial x_j} + \dfrac{\partial u_j}{\partial x_1}) n_j \mathrm{d} s $$

lift force

$$ F_L = -\int_{\Omega_0} f_2 \mathrm{d} \mathbf x = -\int_{\partial \Omega_0} \rho u_2 \mathbf u \cdot \mathbf n \mathrm{d} s - \int_{\partial \Omega_0} p n_2 \mathrm{d} s + \int_{\partial \Omega_0} \mu (\dfrac{\partial u_2}{\partial x_j} + \dfrac{\partial u_j}{\partial x_2}) n_j \mathrm{d} s $$

Therefore, we simply pick any square domain $\Omega_0$ enclosing the cylinder and compute the line integrals of the above equation to obtain the drag force.

My code reproduction is:

#include "calc_solid_force_from_fluid.h"

double dudx(var& u, int i, int j, double h)
{
    double dudx = (u(i + 1, j) - u(i, j)) / h;
    return dudx;
}

double dudy(var& u, int i, int j, double h)
{
    double u_i_j_add_half         = 0.5 * (u(i, j) + u(i, j + 1));         // u(i, j + 1/2)
    double u_i_add_1_j_add_half   = 0.5 * (u(i + 1, j) + u(i + 1, j + 1)); // u(i + 1, j + 1/2)
    double u_i_j_minus_half       = 0.5 * (u(i, j) + u(i, j - 1));         // u(i, j - 1/2)
    double u_i_add_1_j_minus_half = 0.5 * (u(i + 1, j) + u(i + 1, j - 1)); // u(i + 1, j - 1/2)

    double u_i_add_half_j_add_half   = 0.5 * (u_i_j_add_half + u_i_add_1_j_add_half);     // u(i + 1/2, j + 1/2)
    double u_i_add_half_j_minus_half = 0.5 * (u_i_j_minus_half + u_i_add_1_j_minus_half); // u(i + 1/2, j - 1/2)

    double dudy = (u_i_add_half_j_add_half - u_i_add_half_j_minus_half) / h;

    return dudy;
}

double dvdy(var& v, int i, int j, double h)
{
    double dvdy = (v(i, j + 1) - v(i, j)) / h;
    return dvdy;
}

double dvdx(var& v, int i, int j, double h)
{
    double v_i_add_half_j         = 0.5 * (v(i, j) + v(i + 1, j));         // v(i + 1/2, j)
    double v_i_add_half_j_add_1   = 0.5 * (v(i, j + 1) + v(i + 1, j + 1)); // v(i + 1/2, j + 1)
    double v_i_minus_half_j       = 0.5 * (v(i, j) + v(i - 1, j));         // v(i - 1/2, j)
    double v_i_minus_half_j_add_1 = 0.5 * (v(i, j + 1) + v(i - 1, j + 1)); // v(i - 1/2, j + 1)

    double v_i_add_half_j_add_half   = 0.5 * (v_i_add_half_j + v_i_add_half_j_add_1);     // v(i + 1/2, j + 1/2)
    double v_i_minus_half_j_add_half = 0.5 * (v_i_minus_half_j + v_i_minus_half_j_add_1); // v(i - 1/2, j + 1/2)

    double dvdx = (v_i_add_half_j_add_half - v_i_minus_half_j_add_half) / h;

    return dvdx;
}

double integrand_x(var&           u,
                   var&           v,
                   Neumann_field& p,
                   int            i,
                   int            j,
                   double         rho,
                   double         dynamic_viscosity,
                   double         h,
                   double         n1,
                   double         n2)
{
    double Fx = 0.0;

    Fx += -rho * u(i, j) * (u(i, j) * n1 + v(i, j) * n2);
    Fx += -p(i, j) * n1;
    Fx += dynamic_viscosity * (2.0 * dudx(u, i, j, h) * n1 + (dudy(u, i, j, h) + dvdx(v, i, j, h)) * n2);

    return Fx;
}

double integrand_y(var&           u,
                   var&           v,
                   Neumann_field& p,
                   int            i,
                   int            j,
                   double         rho,
                   double         dynamic_viscosity,
                   double         h,
                   double         n1,
                   double         n2)
{
    double Fy = 0.0;

    Fy += -rho * v(i, j) * (u(i, j) * n1 + v(i, j) * n2);
    Fy += -p(i, j) * n2;
    Fy += dynamic_viscosity * (2.0 * dvdy(v, i, j, h) * n2 + (dudy(u, i, j, h) + dvdx(v, i, j, h)) * n1);

    return Fy;
}

std::tuple<double, double> calc_solid_force_from_fluid(NSSolver2D& ns_solver,
                                                       double      rho,
                                                       double      dynamic_viscosity,
                                                       double      xmin,
                                                       double      xmax,
                                                       double      ymin,
                                                       double      ymax)
{
    var&           u = ns_solver.u;
    var&           v = ns_solver.v;
    Neumann_field& p = ns_solver.p;

    double h = ns_solver.h;

    // u

    int ixmin = int(xmin / ns_solver.h);
    int ixmax = int(xmax / ns_solver.h);
    int iymin = int(xmin / ns_solver.h - 0.5);
    int iymax = int(xmax / ns_solver.h - 0.5);

    double Fx = 0.0, Fy = 0.0;

    // Fx

    // left
    for (int i = ixmin, j = iymin; j < iymax; j++)
    {
        Fx += integrand_x(u, v, p, i, j, rho, dynamic_viscosity, h, -1.0, 0.0);
    }

    // right
    for (int i = ixmax, j = iymin; j < iymax; j++)
    {
        Fx += integrand_x(u, v, p, i, j, rho, dynamic_viscosity, h, 1.0, 0.0);
    }

    // down
    for (int i = ixmin, j = iymin; i < ixmax; i++)
    {
        Fx += integrand_x(u, v, p, i, j, rho, dynamic_viscosity, h, 0.0, -1.0);
    }

    // up
    for (int i = ixmin, j = iymax; i < ixmax; i++)
    {
        Fx += integrand_x(u, v, p, i, j, rho, dynamic_viscosity, h, 0.0, 1.0);
    }

    // Fy

    ixmin = int(xmin / ns_solver.h - 0.5);
    ixmax = int(xmax / ns_solver.h - 0.5);
    iymin = int(xmin / ns_solver.h);
    iymax = int(xmax / ns_solver.h);

    // left
    for (int i = ixmin, j = iymin; j < iymax; j++)
    {
        Fy += integrand_y(u, v, p, i, j, rho, dynamic_viscosity, h, -1.0, 0.0);
    }

    // right
    for (int i = ixmax, j = iymin; j < iymax; j++)
    {
        Fy += integrand_y(u, v, p, i, j, rho, dynamic_viscosity, h, 1.0, 0.0);
    }

    // down
    for (int i = ixmin, j = iymin; i < ixmax; i++)
    {
        Fy += integrand_y(u, v, p, i, j, rho, dynamic_viscosity, h, 0.0, -1.0);
    }

    // up
    for (int i = ixmin, j = iymax; i < ixmax; i++)
    {
        Fy += integrand_y(u, v, p, i, j, rho, dynamic_viscosity, h, 0.0, 1.0);
    }

    Fx *= h; // ds in line integrals
    Fy *= h; // ds in line integrals

    return {Fx, Fy};
}

The code directly return what the force I need.

The code consider two main aspects:

  1. There is transform from world space xmin, xmax to index space ixmin, ixmax, which is related with staggered grid.

  2. To get the dudy, dvdx in staggered grid, I interpolate the u, v as figure below.

enter image description here

Other details about solver implementation

SIMPLE algorithm

The solving algorithm of NS equation is SIMPLE algorithm.

for (int time_iter = 0; time_iter < MAX_STEP; time_iter++)
{
    _NS_solver.apply_boundary();
    _NS_solver.cal_conv();
    _NS_solver.apply_pressure_gredient();
    _NS_solver.u2F(cylinder);
    _NS_solver.apply_ib_force(cylinder);
    _NS_solver.apply_boundary();
    _NS_solver.refresh_vars();

    std::cout << "PPE" << std::endl;
    for (int corr_iter = 0; corr_iter < 1; corr_iter++)
    {
        _NS_solver.cal_divu();
        _PE_solver.solve(_NS_solver.Div_u, _NS_solver.p);
        _NS_solver.apply_pressure_gredient();
    }

    ...
}

More specific details can't be provided because they are contributions of our groups, but I think the pipeline is enough.

Immersive boundary model

There are many models of immersive boundary method, what we used is "direct forcing" model, basically it means that the IB force is evaluated by the difference between the actual velocity of solid and the fluid velocity in solid position.

$$ \mathbf{F}(\mathbf{X}_l^{(m)}) = \dfrac{\mathbf{U}^{(d)}(\mathbf{X}_l) - \widetilde{\mathbf{U}}(\mathbf{X}_l)}{\Delta t} \forall l;m, $$

More information see:

Uhlmann, Markus. "An immersed boundary method with direct forcing for the simulation of particulate flows." Journal of computational physics 209.2 (2005): 448-476.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

The question is solved.

I use wrong coefficient formula. I should use 2d formula but actually I use 3d.

And my calculation method using momentum conversation may be wrong, but that is fine, I can use other.

In 2d:

$$ C_L = \dfrac{F_L}{(1/2)\rho u^2_{\infty}D} $$

$$ C_D = \dfrac{F_D}{(1/2)\rho u^2_{\infty}D} $$

where $D$ is feature length.

In 3d:

$$ C_L = \dfrac{F_L}{(1/2)\rho u^2_{\infty}S} $$

$$ C_D = \dfrac{F_D}{(1/2)\rho u^2_{\infty}S} $$

where $S$ is section surface.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.