Kelvin-Helmholtz instability

Numerical simulation of the KH instability with SU2

Geometry, Boundary and Initial conditions

The Euler equations are soved on a square of edge equal to 1. The computational grid is a uniform structured mesh with 1000x1000 points. Periodic boundary conditions are enforced on all sides of the square. The initial condition is set to trigger the instability as:

$$ (\rho, u, v, p) = \begin{cases}(2, -0.5, 0.01 \sin(2 \pi x), 2.5) \quad &for |y| \leq 0.25 \\ (1, 0.5, 0.01 \sin(2 \pi x), 2.5) \quad &for |y| > 0.25 \end{cases} $$

The working fluid is an ideal gas with gamma = 1.4 and R = 1.

Numerical setup

The convective fluxes are discretized using the Approximate Riemann Solver of Roe. The scheme is made second order in space by employing the MUSCL reconstruction. Gradients are computed via the weighted least squares approach, anche the slope of the resconstruction is limited via the Van Albada limiter. A dual time stepping approach is used to evolve the equations in time. The physical time stepping is performed via a BDF2 scheme, making the simulation second order in time. The internal loop of the dual time stepping is integrated with the implicit euler scheme. The internal loop is iterated until the RMS of the residual of the density equation falls below the value of 10^-7. The physical time step used is 0.004.

Imposing initial condition in SU2

The initial condition can be imposed by modifying function CFVMFlowSolverBase::SetInitialCondition in file SU2_CFD/include/solvers/CFVMFlowSolverBase.inl. Ther following code was addedd at line 1027 in version 7.4.0 of SU2:

/*--- Loop over the multigrid levels. ---*/
for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {

    /*--- Loop over all grid points. ---*/
    SU2_OMP_FOR_STAT(omp_chunk_size)
    for (iPoint = 0; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) {

        /* Set the pointers to the coordinates and solution of this DOF. */
        const su2double *coor = geometry[iMesh]->nodes->GetCoord(iPoint);
        su2double *solDOF     = solver_container[iMesh][FLOW_SOL]->GetNodes()->GetSolution(iPoint);

        /* Set KH inital condition. */
	      su2double r0 = 2, r1 = 1, u0 = -0.5, u1 = 0.5, v = 0.01, p = 2.5, v0,v1;
	      if (abs(coor[1]) <= 0.25) {
	          v0 = v*sin(2*M_PI*coor[0]);
	          solDOF[0] = r0;
	          solDOF[1] = r0 * u0;
	          solDOF[2] = r0 * v0;
	          solDOF[3] = (p / 0.4) + 0.5 * r0 * (u0*u0 + v0*v0);
	      } else {
	          v1 = v*sin(2*M_PI*coor[0]);
	          solDOF[0] = r1;
	          solDOF[1] = r1 * u1;
	          solDOF[2] = r1 * v1;
	          solDOF[3] = (p / 0.4) + 0.5 * r1 * (u1*u1 + v1*v1);

	      }
      
    }
    END_SU2_OMP_FOR
}

Simulation result