Partially working Karman Vortexs

BasicSolver/2D-seq/canal.par

@@ -32,7 +32,7 @@ jmax       50	   # number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      100.0   # final time
+te      60.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -62,7 +62,7 @@ y2                  4.0
 # -----------------------
 # Shape 0 disable, 1 Rectangle/Square, 2 Circle
 
-shape               1 
+shape               0 
 xCenter             10.0
 yCenter             2
 xRectLength         6.0

BasicSolver/2D-seq/dcavity.par

@@ -32,7 +32,7 @@ jmax       100		# number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      60.0		# final time
+te      10.0		# final time
 dt      0.02	    # time stepsize
 tau     0.5	    	# safety factor for time stepsize control (<0 constant delt)
 
@@ -42,7 +42,7 @@ tau     0.5	    	# safety factor for time stepsize control (<0 constant delt)
 itermax  1000		# maximal number of pressure iteration in one time step
 eps      0.001		# stopping tolerance for pressure iteration
 rho      0.5
-omg      1.7		# relaxation parameter for SOR iteration
+omg      1.8		# relaxation parameter for SOR iteration
 gamma    0.9		# upwind differencing factor gamma
 
 # Particle Tracing Data:

BasicSolver/2D-seq/karman.par

@@ -0,0 +1,72 @@
+#==============================================================================
+#                            Laminar Canal Flow
+#==============================================================================
+
+# Problem specific Data:
+# ---------------------
+
+name karman             # name of flow setup
+
+bcTop      1			#  flags for boundary conditions
+bcBottom   1			#  1 = no-slip      3 = outflow
+bcLeft     3			#  2 = free-slip    4 = periodic
+bcRight    3			#
+
+gx    0.0      # Body forces (e.g. gravity)
+gy    0.0      #
+
+re    10000.0	   # Reynolds number
+
+u_init    1.0      # initial value for velocity in x-direction
+v_init    0.0      # initial value for velocity in y-direction
+p_init    0.0      # initial value for pressure
+
+# Geometry Data:
+# -------------
+
+xlength    30.0     # domain size in x-direction
+ylength    4.0	   # domain size in y-direction
+imax       200      # number of interior cells in x-direction
+jmax       50	   # number of interior cells in y-direction
+
+# Time Data:
+# ---------
+
+te      100.0   # final time
+dt      0.02    # time stepsize
+tau     0.5     # safety factor for time stepsize control (<0 constant delt)
+
+# Pressure Iteration Data:
+# -----------------------
+
+itermax  200       # maximal number of pressure iteration in one time step
+eps      0.001   # stopping tolerance for pressure iteration
+rho      0.52
+omg      1.75      # relaxation parameter for SOR iteration
+gamma    0.9       # upwind differencing factor gamma
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   60
+startTime           0
+injectTimePeriod    0.5
+writeTimePeriod     0.5
+
+x1                  0.0
+y1                  1.9
+x2                  0.0
+y2                  2.1
+
+# Obstacle Geometry Data:
+# -----------------------
+# Shape 0 disable, 1 Rectangle/Square, 2 Circle
+
+shape               2 
+xCenter             4.0
+yCenter             2.0
+xRectLength         2.0
+yRectLength         1.0
+circleRadius        1.0
+
+#===============================================================================

BasicSolver/2D-seq/src/main.c

@@ -78,7 +78,7 @@ int main (int argc, char** argv)
         
 
         /* Added function for particle tracing. Will inject and advance particles as per timePeriod */
-        //trace(&particletracer, solver.u, solver.v, t);
+        trace(&particletracer, solver.u, solver.v, t);
 
         t += solver.dt;
         nt++;

BasicSolver/2D-seq/src/solver.c

@@ -23,7 +23,7 @@
 #define S(i, j)   s[(j) * (imax + 2) + (i)]
 #define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
 
-static double distance(i, j, iCenter, jCenter)
+static double distance(double i, double j, double iCenter, double jCenter)
 {
     return sqrt(pow(iCenter - i, 2)
                 + pow(jCenter - j, 2) * 1.0);
@@ -133,11 +133,12 @@ void initSolver(Solver* solver, Parameter* params)
     double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
     solver->dtBound  = 0.5 * solver->re * 1.0 / invSqrSum;
 
-    int iCenter = 0, jCenter = 0, iRectLength = 0, jRectLength = 0, radius = 0;
+    double xCenter = 0, yCenter = 0, radius = 0;
     double x1=0, x2=0, y1=0, y2=0;
 
     int* s = solver->s;
     
+
     switch(params->shape)
     {
         case NOSHAPE:
@@ -161,15 +162,15 @@ void initSolver(Solver* solver, Parameter* params)
 
             break;
         case CIRCLE:
-            iCenter = params->xCenter/dx;
-            jCenter = params->yCenter/dy;
+            xCenter = params->xCenter;
+            yCenter = params->yCenter;
             radius  = params->circleRadius;
 
             for (int j = 1; j < jmax + 1; j++) 
             {
                 for (int i = 1; i < imax + 1; i++) 
                 {
-                    if(distance(i, j, iCenter, jCenter) <= radius)
+                    if(distance((i*dx), (j*dy), xCenter, yCenter) <= radius)
                     {
                         S(i, j) = LOCAL;
                     }
@@ -605,6 +606,13 @@ void setSpecialBoundaryCondition(Solver* solver)
             if(S(0,j) == NONE) U(0, j) = 1.0;
         }
     }
+    else if (strcmp(solver->problem, "karman") == 0) 
+    {
+        for (int j = 1; j < jmax + 1; j++) 
+        {
+            U(0, j) = 1.0;
+        }
+    }
 }
 
 void setObjectBoundaryCondition(Solver* solver)

BasicSolver/2D-seq/vector.plot

@@ -1,7 +1,6 @@
-set terminal png size 1800,768 enhanced font ,12
+set terminal png size 3600,768 enhanced font ,28
 set output 'velocity.png'
+set size ratio -1
 set datafile separator whitespace
-set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 5
-
 
 plot 'velocity.dat' using 1:2:3:4:5 with vectors filled head size 0.01,20,60 lc palette

BasicSolver/2D-seq/vis_files/karman_animate.plot

@@ -0,0 +1,12 @@
+unset border; unset tics; unset key;
+set term gif animate delay 20
+set output "trace.gif"
+set xrange [0:30]
+set yrange [0:4]
+set object 1 circle front at 4.0,2.0 size 1.0 fillcolor rgb "black" lw 2
+
+
+do for [ts=0:300] {
+    plot "particles_".ts.".dat" with points pointtype 7
+}
+unset output