Suggest to me for calculation channel flow, please

[Natky]

Hi, I'm trying to test channel flow case on HiFiLES-code, when I calculated periodic channel flow 3D for lamina flow follow set up my input file, eles.cpp and my domain , at calculation round massflux = -60.7985xxx , Ubulk = -9.67758xxx and bodyforce(1) = -1, the code have show “Fatal error” There is a Fatal error:

Fatal error 'NaN in residual, exiting.' at ../src/eles.cpp:1638

For set up and edit “eles.cpp” at line

6401 // case-specific parameters 6402 // periodic channel: 6403 area = 2_pi; 6404 vol = 4_pi*pi;

6406 // periodic hill (HIOCFD3 Case 3.4): 6407 //area = 9.162; 6408 //vol = 114.34; 6409 //mdot0 = 9.126; initial mass flux \ i dont know that 's set up comment "mdot0" , you know? or suggest to me, ****

and

6438 // set body force for streamwise momentum and energy 6439 body_force(0) = 0.; 6440 //bodyforce(1) = alpha/area/dt(mdot0 - mass_flux); // modified SD3D version 6441 //bodyforce(1) = 1.0/area/dt(mdot0 - 2.0_mass_flux + mdot_old); // HIOCFD C3.4 version 6442 body_force(1) = -1.0; //Hack for bodyforce constant 6443 body_force(2) = 0.; 6444 body_force(3) = 0.; 6445 body_force(4) = body_force(1)_ubulk; // energy forcing

Here is my set up domain periodic channel flow and grid element Lx = 0 : 2pi
Ly = 0 : pi Lz = -1 : 1 and cell element : nonuniform grid = ( 24 * 12 \ 28 )

and here is my input file

---

Solver parameters

// 0: Euler/Navier-Stokes, 1:Advection/Adv-Diffusion equation 0 viscous 1

// 0: Rusanov, 1: Lax-Friedrich, 2: Roe riemann_solve_type 0 vis_riemann_solve_type 0

// 0: Isentropic Vortex, 1: Uniform flow, 2: Sine Wave ic_form 6

// 0: Normal (doesn't have an analytical solution), 1:Isentropic Vortex, 2: Advection-Equation test_case 0 order 5 // Order of basis polynomials dt_type 0 // 0: User-supplied, 1: Global, 2: Local dt 0.0001 CFL 1 n_steps 1000000 adv_type 3 // 0: Forward Euler, 3: RK45 tau 1.0 pen_fact 0.5

---

Mesh options

mesh_file ch3d_pi.msh dx_cyclic 6.283185307179586 dy_cyclic 3.141592653589793 dz_cyclic 2.0

---

Monitoring, plotting parameters

p_res 4 // Plotting resolution, # of nodes per direction write_type 0 // 0: Paraview, 1: Tecplot // Choose extra fields to be written to file: u v w energy pressure mach vorticity q_criterion. // Set to 0 or comment out for no diagnostic fields n_diagnostic_fields 5 u v w pressure mach inters_cub_order 4 // Order of cubature rule for integrating over element interfaces volume_cub_order 4 // Order of cubature rule for integrating over element volumes

monitor_res_freq 100 plot_freq 500 data_file_name ch3dv01zy restart_dump_freq 10000 res_norm_type 2 // 1:L1 norm, 2:L2 norm, 0:infinity norm error_norm_type 2 // 1:L1 norm, 2:L2 norm, 0:infinity norm res_norm_field 0 // Density

---

Element parameters

==== Tris ==== upts_type_tri 0 fpts_type_tri 1 vcjh_scheme_tri 3 c_tri 0.0 sparse_tri 0

==== Quads ==== upts_type_quad 0 // 0: Gauss, 1: Gauss-Lobatto vcjh_schemequad 1 // 0: VCJH, 1: DG, 2: SD, 3: Hu, 4: c+ eta_quad 0.0 sparse_quad 0

==== Hexs ==== upts_type_hexa 0 vcjh_scheme_hexa 1 eta_hexa 0. sparse_hexa 0

==== Tets ==== upts_type_tet 1 fpts_type_tet 0 vcjh_scheme_tet 1 eta_tet 0.0 sparse_tet 0

==== Prisms ==== upts_type_pri_tri 0 upts_type_pri_1d 0 vcjh_scheme_pri_1d 1 eta_pri 0.0 sparse_pri 0

---

Fluid Parameters

gamma 1.4 prandtl 0.72 S_gas 120. T_gas 291.15 R_gas 286.9 mu_gas 1.827E-05

---

Boundary conditions

==== Viscous ==== fix_vis 0 // 0: Sutherland's law, 1: Constant viscosity Mach_free_stream 0.2 nx_free_stream 1. ny_free_stream 0. nz_free_stream 0. Re_free_stream 180. L_free_stream 1. T_free_stream 300 T_wall 300 //Mach_wall 0. //nx_wall 1. //ny_wall 0. //nz_wall 0. //T_wall 300

==== Inviscid ==== rho_bound 1. u_bound 1.0 v_bound 1.0 w_bound 0.0 p_bound 17.857142857142854098

---

Forcing

body_forcing 1 // Used for the periodic channel test case. 0 = no forcing, 1 = forcing

Initial Profile: If ic_form=6, use these parameters to describe polynomial initial condition. Useful for cyclic problems. Specify coeffs of poly of form a + bx + cx^2 + dx^3 + ex^4 + by + cy^2 + dy^3 + ey^4 + bz + cz^2 + dz^3 + ez^4

// Coeffs of x-momentum [need the 13 so reader know how many quantities] x_coeffs 13 1.875 0 0 0 0 0 0 0 0 0 -3.75 0 1.875
// Coeffs of y-momentum y_coeffs 13 0 0 0 0 0 0 0 0 0 0 0 0 0 // Coeffs of y-momentum z_coeffs 13 0 0 0 0 0 0 0 0 0 0 0 0 0 perturb_ic 1 // If 1, add velocity perturbation to kickstart turbulence in periodic channel. 0 otherwise.

Please, suggest to me thank you,

suggest to me or contact to me : nattakit_ch@cmu.ac.th

[mlopez14]

The timestep (1e-4) seems too large for the order (5) of the scheme. I would suggest you try a lower order and a smaller timestep. Let us know if this error happens at order = 1 and dt less than 1e-6.

[Natky]

Thank you for your suggestion and I have set timestep = 1e-6 and set to use the order = 1 of the scheme, but The code can't calculate periodic channel flow. Could you please check the setting in eles.cpp and my input file for me again?

And there is a fatal error:

---

iter, mdot0, mdot_old, mass_flux, body_force(1): 4538, 0.00000000, -60.39439814, -60.42285522, -1.00000000 iter, mdot0, mdot_old, mass_flux, body_force(1): 4539, 0.00000000, -60.42285522, -60.45131858, -1.00000000 iter, mdot0, mdot_old, mass_flux, body_force(1): 4540, 0.00000000, -60.45131858, -60.47978820, -1.00000000

Fatal error 'NaN in residual, exiting.' at ../src/eles.cpp:1638

---

[mlopez14]

Correct me if I'm wrong, but it seems the flow is being initialized going in the +x direction (nx_free_stream 1.) and the forcing function is returning a negative value, so the forcing is going in the -x direction. I would try initializing the flow to go in the -x direction (nx_free_stream -1.).

[Natky]

Thank you for suggestion. But i found the problem during code calculation same as before. I think the problem causes from uncomplete set up at 'eles.cpp' file.Could you please give me a suggestion about setting up the code for calculation periodic channel flow and then i will use the code to develop the application for more calulation on periodic turbulent channel flow domain. if you have any manual or the code user guide, could you please let me know or give it to me?

Thank you very much in advance.

[mlopez14]

The code user guide is the wiki we have here in github. The forcing function is experimental and was used only for one instance that is not published yet. In what part of the domain is the code blowing up? Is it close to the boundary, right at the middle? As long as the boundary is labeled periodic in the mesh and the correct dx_cyclic, etc. variables are in the input file (they seem to be correct), the code would let you know if there is a mismatch.

Do the elements on the periodic boundaries match as if they were adjacent? If they do not, the code would be blowing up at the boundaries.

[Natky]

Thank you for suggestion.