[Scabbard]

Originally Posted by aerosaravanan

Hi Scabbard,
I tried to simulate for Re 3900 with low reynolds model(LaunderSharmakepsilon). This time I got reasonable strouhal Number but high drag value.It's confusing me lot.
Then for your last post same thing happening for me.I got few papers on circular cylinder analysis in OpenFOAM for wide range of Re.And they got reasonable values like CD, strouhal number with komegaSST model.

Hi aerosavanvan,

I am really sick with this, I almost follow all the settings in the paper, but can not get even one correct result. This situation last for a month... Hope some one can help me solve this

Best regards,
Scabbard

[davibarreira]

Hey guys, did you manage to get a solution for this problem? I did simulations for a wide range of Reynolds numbers and compared the Drag Coefficiente (Cd), the lift (Cl) and the Strouhal number (St) with values from the literature.
The standard k-epsilon breaks in the region between 1e4 and 1e6 (the results for Re = 1e6 are ok, the Cd is spot on, also for Re = 3900 and 1000 the results are pretty good). I also tried k-omega SST and kklOmega. The k-omega gave me values of Cd around 0.9 for Re 1e4 and 1e5, but the experimental is around 1.2... Also the Cl values were very high (1.12 compared to 0.5 expected from literature). The kklOmega showed to be quite sensitive to turbulence initial boundary values.
Any thoughts?

[arieljeds]

Hi guys, I know this thread has been quiet for a while but I'm running a similar case for a wide variety of Re numbers but I haven't been able to get a drag coefficient even close to what I'm expecting. So, a couple of questions:

1) How many iterations were necessary to even get a drag coefficient of 1.65? Mine seems to be around 80!

2) Can anyone give some advice on where I might be going wrong? I've attached my forces object and controlDict files and my blockMestDict because from looking at @Scabbard and @aerosaravanan 's cases, my initial and boundary conditions and discretisation schemes look the same (or almost the same)

Please help! I'm getting very confused and have been spending too long trying to understand where the problem is!

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

D	     6;
depth       -10;
air          10;	        	 // Distance above free surface in computational domain  


Lx1     #calc "40.0*$D";	     	 // Unscaled length of the tank (downstream)
Lx2     #calc "-20.0*$D"; 	     	 // Unscaled length of tank (upstream)
Ly1     #calc "10*$D";           	 // Unscaled width of tank/2 (positive)
Ly2     #calc "-10*$D";        	 // Unscaled width of tank/2 (negative)


r1      #calc "$D/2";      	  	 // Positive radius location 
r2      #calc "-1.0*($D/2)";         	 // Negative radius location 
o_r1     7;	        		 // Outer circle radius location 1
o_r2     -7	;	         // Outer circle radius location 2

// Locations on the inner cylinder and outer boundary layer finer mesh 
ang1  	   45;		// Corresponds to pi/4
ang2  	  135;		// Corresponds to 3*pi/4
ang3  	  225;		// Corresponds to 5*pi/4
ang4  	  315;		// Corresponds to 7*pi/4
ang1_rad  #calc "degToRad($ang1)";
ang2_rad  #calc "degToRad($ang2)";
ang3_rad  #calc "degToRad($ang3)";
ang4_rad  #calc "degToRad($ang4)";
p1y       #calc "$o_r1*sin($ang2_rad)";		// Location on outer circle, positive y. Corresponds to both pi/4 and 3*pi/4
p2y	  #calc "$o_r1*sin($ang4_rad)";         // Location on outer circle, corresponds to 5*pi/4 and 7*pi/4
p1x 	  #calc "$o_r1*cos($ang1_rad)";         // Location on outer circle, corresponds to pi/4 and 7*pi/4
p2x	  #calc "$o_r1*cos($ang3_rad)";         // Location on outer circle, corresponds to 3*pi/4 and 5*pi/4
p3y       #calc "$r1*sin($ang1_rad)";		// Location on inner circle, y value. Corresponds to pi/4 and 3*pi/4
p4y	  #calc "$r1*sin($ang4_rad)";		// Location on inner circle, y value. Corresponds to 5*pi/4 and 7*pi/4
p3x	  #calc "$r1*cos($ang1_rad)";		// Location on inner circle, x value. Corresponds to 3*pi/4 and 7*pi/4
p4x	  #calc "$r1*cos($ang3_rad)"; 		// Location on inner circle, x value. Corresponds to 3*pi/4 and 5*pi/4

// Locations on the circles to create the arc
ang5       20;					// Angle to pick a location on the circle to create arcs
ang5_rad   #calc "degToRad($ang5)";
arc1	   #calc "$r1*cos($ang5_rad)";		// Point one, inner circle (cylinder), positive values
arc2       #calc "$r1*sin($ang5_rad)";		// Point two, inner circle (cylinder), positive values
arc1_n	   #calc "$r2*cos($ang5_rad)";		// Point one, inner circle (cylinder), negative values
arc2_n     #calc "$r2*sin($ang5_rad)";		// Point two, inner circle (cylinder), negative values 
arc3	   #calc "$o_r1*cos($ang5_rad)";	// Point one, outer circle, positive values 
arc4       #calc "$o_r1*sin($ang5_rad)";	// Point two, outer circle, positive values 
arc3_n     #calc "$o_r2*cos($ang5_rad)";        // Point one, outer circle, negative values 
arc4_n     #calc "$o_r2*sin($ang5_rad)";        // Point two, outer circle, negative values  
 
 // Values for number of cells
 zn	   30;					// Change this number to change the  number of cells in the vertical direction 
 bl	   2;					// Change this to increase the grading in the boundary layer (1 = no grading)
 ySep	   10; 
 nUp	   #calc "std::ceil(((-1.0)*$Lx2+$o_r2)/(3*($D)/16))";		// Cell length upstream, to make upstream mesh finer in oscillatory flow, divide by 3D/16 instead of D/4
 nDown     #calc "std::ceil(($Lx1-$o_r1)/((3*$D)/16))";			// Cell length downstream of the cylinder
 nBL	   #calc "std::ceil(($o_r1-$r1)/((atan(1.000)*4*$D)/76))";		// Cell length in the boundary layer. NOTE: replace atan(1.000)*4 with pi. To make the mesh finer, increase 76
 nROI	   #calc "std::ceil(($Ly1-$o_r1)/((3.0*$D)/20))";			// Cell length 

 	
vertices
(
   ( $r1    0    $depth )		// 0
   ( $o_r1  0    $depth )		// 1
   ( $Lx1   0    $depth )		// 2
   ( $Lx1  $p1y  $depth )		// 3		
   ( $p1x  $p1y  $depth )		// 4
   ( $p3x  $p3y  $depth )		// 5
   ( $Lx1  $Ly1  $depth )		// 6 	
   ( $p1x  $Ly1  $depth )		// 7
   (   0   $Ly1  $depth )		// 8
   (   0   $o_r1 $depth )		// 9
   (   0   $r1   $depth )		// 10
 
   ( $r2    0    $depth )	        // 11
   ( $o_r2  0    $depth )		// 12
   ( $Lx2   0    $depth )		// 13
   ( $Lx2  $p1y  $depth )	   	// 14
   ( $p2x  $p1y  $depth )		// 15
   ( $p4x  $p3y  $depth )		// 16
   ( $Lx2  $Ly1  $depth )		// 17 
   ( $p2x  $Ly1  $depth )		// 18

   ( $r1    0    $air )			// 19 
   ( $o_r1  0    $air ) 		// 20
   ( $Lx1   0    $air )	  		// 21
   ( $Lx1  $p1y  $air )			// 22
   ( $p1x  $p1y  $air )			// 23
   ( $p3x  $p3y  $air )			// 24
   ( $Lx1  $Ly1  $air )			// 25 
   ( $p1x  $Ly1  $air )			// 26
   (   0   $Ly1  $air )			// 27
   (   0   $o_r1 $air )			// 28
   (   0   $r1   $air )			// 29

   ( $r2    0    $air )			// 30
   ( $o_r2  0    $air )			// 31
   ( $Lx2   0    $air )			// 32
   ( $Lx2  $p1y  $air )			// 33
   ( $p2x  $p1y  $air )			// 34
   ( $p4x  $p3y  $air )			// 35
   ( $Lx2  $Ly1  $air )			// 36
   ( $p2x  $Ly1  $air )			// 37
	
   ( $Lx1  $p2y  $depth )		// 38 
   ( $p1x  $p2y  $depth )		// 39 
   ( $p3x  $p4y  $depth )		// 40 
   ( $Lx1  $Ly2  $depth )		// 41
   ( $p1x  $Ly2  $depth )		// 42
   (   0   $Ly2  $depth )		// 43
   (   0   $o_r2 $depth )		// 44
   (   0   $r2   $depth )		// 45
 	
   ( $Lx2  $p2y  $depth )	        // 46 
   ( $p2x  $p2y  $depth )		// 47 
   ( $p4x  $p4y  $depth )		// 48 
   ( $Lx2  $Ly2  $depth )		// 49 
   ( $p2x  $Ly2  $depth )		// 50 

   ( $Lx1  $p2y   $air )		// 51
   ( $p1x  $p2y   $air )		// 52
   ( $p3x  $p4y   $air )		// 53
   ( $Lx1  $Ly2   $air )		// 54
   ( $p1x  $Ly2   $air )		// 55
   (  0    $Ly2   $air )		// 56
   (  0    $o_r2  $air )		// 57 
   (  0    $r2    $air )		// 58
	
   ( $Lx2  $p2y   $air )		// 59
   ( $p2x  $p2y   $air )		// 60 
   ( $p4x  $p4y   $air )		// 61
   ( $Lx2  $Ly2   $air )		// 62
   ( $p2x  $Ly2   $air )		// 63
);

blocks
(  
   hex (5 4 9 10 24 23 28 29)    ( $nBL   $ySep  $zn ) simpleGrading ( $bl 1  1 )	// 0
   hex (0 1 4 5 19 20 23 24)     ( $nBL   $ySep  $zn ) simpleGrading ( $bl 1  1 )	// 1 
   hex (1 2 3 4 20 21 22 23)     ( $nDown $ySep  $zn ) simpleGrading (  1  1  1 )	// 2
   hex (4 3 6 7 23 22 25 26)     ( $nDown $nROI  $zn ) simpleGrading (  1  1  1 )	// 3
   hex (9 4 7 8 28 23 26 27)     ( $ySep  $nROI  $zn ) simpleGrading (  1  1  1 )	// 4	
   hex (16 10 9 15 35 29 28 34)  ( $ySep  $nBL   $zn ) simpleGrading (  1 $bl 1 )	// 5 
   hex (11 16 15 12 30 35 34 31) ( $ySep  $nBL   $zn ) simpleGrading (  1 $bl 1 )	// 6
   hex (12 15 14 13 31 34 33 32) ( $ySep  $nUp   $zn ) simpleGrading (  1  1  1 )	// 7
   hex (15 18 17 14 34 37 36 33) ( $nROI  $nUp   $zn ) simpleGrading (  1  1  1 )	// 8
   hex (9 8 18 15 28 27 37 34)   ( $nROI  $ySep  $zn ) simpleGrading (  1  1  1 )	// 9


   hex (40 45 44 39 53 58 57 52)  ( $ySep  $nBL   $zn ) simpleGrading (  1 $bl 1 )	// 10
   hex (0 40 39 1 19 53 52 20)    ( $ySep  $nBL   $zn ) simpleGrading (  1 $bl 1 )	// 11
   hex (1 39 38 2 20 52 51 21)    ( $ySep  $nDown $zn ) simpleGrading (  1  1  1 )	// 12 
   hex (39 42 41 38 52 55 54 51)  ( $nROI  $nDown $zn ) simpleGrading (  1  1  1 )	// 13
   hex (44 43 42 39 57 56 55 52)  ( $nROI  $ySep  $zn ) simpleGrading (  1  1  1 )	// 14
   hex (48  47 44 45 61 60 57 58) ( $nBL   $ySep  $zn ) simpleGrading ( $bl 1  1 )	// 15
   hex (11 12 47 48 30 31 60 61)  ( $nBL   $ySep  $zn ) simpleGrading ( $bl 1  1 )	// 16
   hex (12 13 46 47 31 32 59 60)  ( $nUp   $ySep  $zn ) simpleGrading (  1  1  1 )	// 17
   hex (47 46 49 50 60 59 62 63)  ( $nUp   $nROI  $zn ) simpleGrading (  1  1  1 )	// 18 
   hex (44 47 50 43 57 60 63 56)  ( $ySep  $nROI  $zn ) simpleGrading (  1  1  1 )	// 19 
);

edges
(  
   arc 0 5   ( $arc1   $arc2  $depth )	// 0       --> This is at 20 degrees
   arc 5 10  ( $arc2   $arc1  $depth )	// 1	   --> 
   arc 1 4   ( $arc3   $arc4  $depth )	// 2
   arc 4 9   ( $arc4   $arc3  $depth )	// 3
   arc 19 24 ( $arc1   $arc2  $air )	// 4
   arc 24 29 ( $arc2   $arc1  $air )	// 5
   arc 20 23 ( $arc3   $arc4  $air )	// 6 
   arc 23 28 ( $arc4   $arc3  $air )	// 7 
   arc 11 16 ( $arc1_n $arc2  $depth )	// 8 
   arc 16 10 ( $arc2_n $arc1  $depth )	// 9 
   arc 12 15 ( $arc3_n $arc4  $depth )	// 10 
   arc 15 9  ( $arc4_n $arc3  $depth )	// 11
   arc 30 35 ( $arc1_n $arc2  $air )	// 12
   arc 35 29 ( $arc2_n $arc1  $air )	// 13
   arc 31 34 ( $arc3_n $arc4  $air )	// 14
   arc 34 28 ( $arc4_n $arc3  $air )    // 15

   // completing the circle. Commented values correspond to the first half, but with negative y values 	
   arc 0 40  ( $arc1   $arc2_n  $depth ) // 0
   arc 40 45 ( $arc2   $arc1_n  $depth ) // 1 
   arc 1 39  ( $arc3   $arc4_n  $depth ) // 2
   arc 39 44 ( $arc4   $arc3_n  $depth ) // 3 
   arc 19 53 ( $arc1   $arc2_n  $air )	 // 4
   arc 53 58 ( $arc2   $arc1_n  $air ) 	 // 5
   arc 20 52 ( $arc3   $arc4_n  $air )	 // 6
   arc 52 57 ( $arc4   $arc3_n  $air )	 // 7
   arc 11 48 ( $arc1_n $arc2_n  $depth ) // 8
   arc 48 45 ( $arc2_n $arc1_n  $depth ) // 9
   arc 12 47 ( $arc3_n $arc4_n  $depth ) // 10 
   arc 47 44 ( $arc4_n $arc3_n  $depth ) // 11
   arc 30 61 ( $arc1_n $arc2_n  $air )	 // 12
   arc 61 58 ( $arc2_n $arc1_n  $air )	 // 13
   arc 31 60 ( $arc3_n  $arc4_n  $air )	 // 14
   arc 60 57 ( $arc4_n $arc3_n  $air ) 	 // 15

);

boundary
(
   
   top
   {
       type slip;
       faces
       (
           (7 8 27 26)
           (6 7 26 25)
           (8 18 37 27)
           (18 17 36 37)
       );
   }   
       
   bottom
   {
       type slip;
       faces
       (
          (49 50 63 62)
          (50 43 56 63)
          (43 42 55 56)
          (42 41 54 55)
       );
   }

   inlet
   {
       type patch;
       faces
       (
           (14 13 32 33)
           (17 14 33 36)

           (46 13 32 59)
           (46 49 62 59)
       );
   }

   outlet
   {
       type patch;
       faces
       ( 
           (2 3 22 21)
           (3 6 25 22)

           (38 51 21 2)
           (41 54 51 38)
       );
   }

   cylinder
   {
       type wall;
       faces
       (
           (10 5 24 29)
           (5 0 19 24)
           (16 10 29 35)
           (11 16 35 30)

           (48 11 30 61)
           (45 48 61 58)
           (40 45 58 53)
           (0 40 53 19) 
       );
   }
   
   
   atmosphere
   {
       type patch;
       faces
       (
       	  (24 23 28 29) // f
       	  (19 20 23 24) // f
       	  (20 21 22 23) // f
       	  (23 22 25 26) // f
       	  (28 23 26 27) // f
       	  (35 29 28 34) // f
       	  (31 30 35 34) // f
       	  (32 31 34 33) // f
       	  (33 34 37 36) // f
       	  (34 28 27 37) // f 
       	  
      
          (58 57 52 53) // f
          (53 52 20 19)	// f
          (52 51 21 20)	// f
          (52 55 54 51)	// f
          (57 56 55 52)	// f
          (60 57 58 61)	// f
          (31 60 61 30)	// f
          (32 59 60 31) // f
          (62 63 60 59)	// f
          (63 56 57 60)	// f
       );
   }
   
   seaFloor
   {
       type patch; 
       faces
       (
          (5 10 9 4)  	// b          
          (0 5 4 1) 	// b        
          (1 4 3 2) 	// b         
          (4 7 6 3)  	// b          
          (4 9 8 7) 	// b      
          (16 15 9 10) 	// b         
          (12 15 16 11) // b     
          (13 14 15 12) // b     
          (14 17 18 15) // b        
          (15 18 8 9) 	// b
          
          (45 40 39 44) // b         
          (40 0 1 39) 	// b    
          (39 1 2 38)	// b         
          (39 38 41 42)	// b         
          (44 39 42 43)	// b          
          (47 48 45 44)	// b         
          (12 11 48 47)	// b        
          (13 12 47 46)	// b	        
          (49 46 47 50)	// b
          (50 47 44 43)	// b
       );
   } 
   
);

mergePatchPairs
(
);

// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    application		waveFoam;
    
    startFrom		latestTime; //startTime; 
     
    startTime		0;
    
    stopAt		endTime;
    
    endTime		730;
    
    deltaT		0.001;
    
    writeControl	adjustableRunTime;
    
    writeInterval	0.2; 
    
    purgeWrite		0;
    
    writeFormat		ascii;
    
    writePrecision	6;
    
    writeCompression	off; 
    
    timeFormat		general;
    
    timePrecision	6;
    
    runTimeModifiable	true;
    
    adjustTimeStep	yes; 
    
    maxCo		0.8; // increase the Co as much as possible 
    
    maxAlphaCo		0.8;
    
    maxDeltaT		1;
    
    
   functions
   {
      #includeIfPresent "../waveGaugesNProbes/surfaceElevationAnyName_controlDict";
      #include "forces_object";
      #include "forceCoeffs_object";
      minmaxdomain
      {
	type fieldMinMax;
	//type banana;

	functionObjectLibs ("libfieldFunctionObjects.so");

	enabled true;

	mode component;

	outputControl timeStep;
	outputInterval 1;

	log true;

	fields (p p_rgh U);
       }

   }

// ************************************************************************** //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      forceCoeffsIncompressible;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


forceCoeffs_object
{
	type 	forceCoeffs;
	functionObjectLibs ("libforces.so");
	patches (cylinder);

	pName 	p;
	Uname 	U;
	rhoName rhoInf;				// Reference Density 
	rhoInf  1.0;				// Density/1000
	
	log 	true;				// Dump to file -- change to false to not save a file 

	CofR 	(0.0 0 0);			// Centre of rotation 
	liftDir (0 1 0);			// Direction of lift coefficient  ---> y-direction (transverse to flow)   ?????
	dragDir (1 0 0);			// Direction of drag coefficient  ---> drag acts in opposite direction to flow (-x) 
	pitchAxis (0 0 1);			// Pitch moment axis  --> changed from pitchAxis (0 0 1)
	magUInf 1.1556;				// free stream velocity magnitude 
	lRef 6.0;         			// reference length (cylinder diameter)
	Aref 6.0;         			// reference area: D

        outputControl   timeStep;
        outputInterval  1;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      forcesIncompressible;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

forces_object
{
    type            forces;  
    functionObjectLibs ( "libforces.so" );

    enabled         true;
    outputControl   timeStep;
    outputInterval  2;

    pName           p;
    UName           U;	
    rhoName         rhoInf; 	// Incompressible solver
    log             on;

    
    rhoInf      1.0;    	// Fluid density
    patches     ( cylinder); 

    CofR        (0 0 0);
    pitchAxis   (0 0 1);	// Gives moment coefficient

   // #includeEtc "caseDicts/postProcessing/forces/forces.cfg"
}

// ************************************************************************* //

Any advice anyone has would be seriously appreciated, I've been looking everywhere and have not made progress on this at all and am getting very frustrated!

[davibarreira]

What are you simulating exaclty? In this thread, the simulations are 2D cylinders, using solvers such as pisoFoam and pimpleFoam. I see you are using waveFoam, so it's going to be hard to comment your results.

Post a picture of your simulation. Also, what is the Reynolds number and turbulence model?

[arieljeds]

Hi Davidberreira,

Yes, admittedly I realise this isn't quite the right place to post it but I've been posting this question in several different threads that I thought might be more relevant but have had no replies, so thought I'd chance my luck with this one..I hope that's ok!

EDIT: I didn't answer this question: I am trying to simulate the forces in the very near wake and directly on the cylinder and capture any features that might be present there (i.e. vorticity) as well as compute the drag and lift coefficients for comparisons to other results.

The Reynolds number is 6.93e6 and I am trying to use a k-omega SST turbulence model.

I'm posting a picture of what the p_rgh and U contour plots look like, although this is after only 328 iterations (265.4 seconds, I was not able to run it for longer due to space).

The simulations "looks" correct and I also seem to be getting appropriate velocity values. I have probes set up in the wave field and around the cylinder. Additionally, I am seeing what I was expecting for this case, which is strong vorticity on the cylinder but no vortex shedding (which makes sense because KC = 0.5 only).

Please let me know if you want me to provide you with any more information that might help..

Thanks a lot for the response!!

[davibarreira]

Ok. So it's hard to say what is causing the bad results. But from my expirience on the 2D case with pisofoam, here is what may be the reason.
First: Look in the literature if people have tried this before, and if they get a close drag coeff. It might be that the model is just not good for your case (for example, for Re 1e5, kwsst gives drag 0.9, but the experimental is 1.2; and the simulation is correct, but the model is that is just not good enough)
Second: mesh refinement. Can you post a picture of your mesh? For external flow the cells close to the surface are very important to get good estimates for drag and lift. Check the value of y+, use yplus application in OpenFoam and report what it gives.
Third: Domain size. Lookin at your picture, it seems large enough, but it might be a source of error. By the way, what is the drag you are getting and what is the experimental?
Fourth: Turbulence parameters in the inlet. Try to verify if your simulation is sensible to them and if the results improve with you perhaps increase or decrease Turbulence intensity.

[arieljeds]

Hi Davibarreira,

Again, thank you for your detailed response. Ok so I've attached three images of my mesh to show you the full layout, an up close of the cylinder/boundary layer, and the vertical cells as well.

Yes in the literature people say they are getting the right drag coefficient, but in one paper, they do not use a turbulence model at all. Also, in the related literature, the focus is on loading on the front of the turbine, so the features in the wake are not as important (given the Reynolds number). So the main problem so far is that I can't really tell you what drag coefficient I'm getting because the case hasn't converged but if I calculate the mean drag coefficient, I am getting crazy answers that are almost not worth mentioning (777.5 is the current average!!)

I'm also attaching my plot of forces and of the drag coefficient. Is it possible that it's correct but just needs far more iterations?? Please note: I have not been able to run it for as long as I'd like due to storage issues, which I am currently solving.

[davibarreira]

I've never user waveFoam, so don't know exaclty how it works. Having said that...

A drag coefficient of 700 seems too high (again, what is the expected value?). My guess would be this, your using a solver that uses density (since you are solving p_rgh), so you might have a value for rho in your constant directory, which should be around 1e-6... So your drag coefficient is being divided by such rho, which is causing the value to be this high. (in the formula for drag coefficient you divide by rho, so...).

To check if the problem is in the simulation or in some mistake in the estimation of the drag coefficient, you might want to try to compute the value from Paraview. You can do a visual check by looking at the values of pressure around the cylinder to see if they are consistent. If your pressure has indeed a very high value, this explains the high drag coefficient, so it might be the density problem or some error in the simulation. If the pressure value looks correct, then it might be some mistake in the dictionary for calculating the forces.

Also, your mesh seems a bit coarse (not enough to get the 700), so after you solve the really high drag values, you might want to do a mesh independence analysis.

Cheers.

[arieljeds]

Yes, the drag coefficient is far too high. Using an analytical expression to calculate the drag (drag for Re = 4e6 drops down significantly), I'm finding that I should only expect a drag value somewhere near 0.0086.

Your explanation though is very clear... I'm trying to understand how the pressures are calculated because I'm having sort of the same issue in trying to get the pressures to be values that I expect (and so I was suspecting that rho is being divided through somewhere). Just to understand your explanation, given that I have p_rgh, my pressures are calculated as just dphi/dt (given that it usually would be p = rho*dphi/dt, and we are dividing by rho). Am I understanding this, as far as you know?

Also, I do suspect that the issue has to do with my forces and forceCoeffs objects because visually, the pressures look correct. Could you have a look at the forceCoeffs_object and forces_object (I included them in my original post).

Lastly, I'm not sure how to do a mesh independence analysis...


Thank you again for your help, I really appreciate it

[davibarreira]

So, your force_coeff dictionary is using rhoinf = 1, which means you are not considering the actual rho for the fluid... This means that if your pressure is the absolute pressure and not the relative pressure, this might be the cause for your really high value. So just change rhoinf to the rho you have in your transportProperties file.

For the mesh independence, to be very brief, just increase the number of cells and plot the results to see if they converge as the mesh gets more refined. I'm putting an example of mesh independence analysis.

P.S: in the figure, "malha" means mesh (it's in portugues)

[arieljeds]

Hi again,

Ok that's very helpful I think... again, thank you for your help, really really appreciate your explanations.

So at the moment, I'm using pressure boundary condition that assumes zerogradient at the surface and am setting the internal field to 0 (in my 0/p_rgh file), and then, as you see I have my rhoInf at 1.0.

Are you suggesting that I either: set rhoInf to 1028 (density in salt water) in the forces function object, or else change the internal field to 1028 in my 0/p_rgh file?

I am not at my work computer but will see if I can do a mesh dependence analysis based on your plots.

[arieljeds]

Whoops, I guess I meant leave it in the transport properties as 1028 and adjust in rhoInf... One question though -- isn't rhoInf irrelevant for incompressible flow?

[davibarreira]

It is, this is why you usually leave at 1, because your pressure should be already divided by rho. But since your solver uses p_rgh, it needs the value of rho, and I would guess the pressure used to solve for the drag is not divided by rho yet.

Your solver probably solves both for p_rgh and p. Your forcecoeffs uses p, not p_rgh. I'm saying for you to change rhoinf to 1028.

Just try it, and see if it improves the results. As I said, I haven't worked with this solver, so I might be wrong.

[arieljeds]

Ah ok, yes I understand what you're saying. Thanks, I've changed it now so we'll see if that makes a difference. (Also I get now why it would be redundant for incompressible if it was solving for p... once again, thank you for a clear explanation). What you're saying I think does make sense and I'll let you know if it works out with this solver.

Otherwise though, do you have any comment on the Aref value? I'm setting it to Aref = 2*pi*r*h, where h is the water depth = 25. However, in my computational domain, I have only set the depth to reach down to -10 m to save time. Is this a problem, as far as you know? Or should the h value for Aref be the cell height, and not the entire cylinder?

Again, thank you so much for your replies, you've been a great help

[arieljeds]

Daviberreira - I have another quick query that you may be able to answer.. In both function objects, you sete the pName. Should I be setting my p_name to p_rgh instead of p then, if I'm solving for p_rgh?

I changed the rhoInf and am still seeing fluctuating Cd that goes as high as 400. I'm playing around with Aref during run-time to see if I see significant changes in the Cd.

[davibarreira]

Did anything change when you altered the rhoInf?

Before running the simulation you should understand the models and its variables. I cant say which of the variables you should use, as I havent worked with this solver and this kind of simulation. P_rgh is just Pressure - rho*g*h. The solver probably solves for p and p_rgh, but you just need to set p_rgh (the purpose of this variable is to make easier to prescribe boundary conditions).

Did you integrate te pressure values yourself to see if the results are right or also too high?

To adress the question about the Aref. What forcecoeffs dictionary does is to integrate your pressure values over the surface you defined, then apply the drag coefficient formula with the Aref you prescribed. So it should be clear to you whether it makes sense or not to procede as you said. Are you sure you are using this solver correclty? This wavefoam is solving for h? If Yes, then what is the Z axis for?

[arieljeds]

Hi Daviberreira,

Apologies for not writing back to you sooner but yes, it did alter it. So I still have an oscillating drag coefficient but when I calculate the average, it is now down to -1. However, the drag is showing no sign of converging to one value (newest plot is attached).

To be honest I can't be 100% sure that I am using it correctly as I'm still just learning but I am fairly positive that I am. I am able to find the surface elevation, which is as expected. I'm int he process of running another case and then I will integrate pressure myself. Do you use the patchIntegrate command? I am going to try this and then compare...

Also, I messed around with the Aref value while it was running to see if I could see any differences but I'm not sure that I can do that and see any changes... do you know anything about changing the values at run time?

[davibarreira]

Sorry, dont knoe about changing at run. Check if you are using the solver properly