2-D Boundary Layer CFD Finite Difference Implicit Solver with Arduino

1// Code to Solve the 2D Boundary Layer problem with Finite Difference Implicit Solver for Arduino Project Hub
2//
3//Written by Adrianos  Botis: adrianosbotis@gmail.com
4//
5//Use at your own and change the parameters:  nx,ny,fh,Rel to change the fluid and grid properties
6// Keep in mind that the space discretization is kept low to avoid memory problems on the arduino board
7//
8
9#include "Plotter.h"
10
11const int nx = 25;
12const int ny = 25;
13float nxc = 25;
14float nyc = 25;
15float fh = 0.097;
16float Rel = 6062.0;
17float v[nx][ny];
18float u[nx][ny];
19float Su[ny];
20float CN[ny];
21float CP[ny];
22float CS[ny];
23float x[nx];
24float y[ny];
25float Rex[nx];
26float dx = 1.0;
27float dy = 1.0;
28float uvel = 0.0;
29float xcord = 0.0;
30Plotter p; // create plotter
31void setup() {
32
33   p.Begin(); // start plotter
34  
35   p.AddXYGraph( "U velocity profile at different X positions", 1500, "Uvel [m/s]", uvel, "y coordinate", xcord); // Add XY plot
36}
37
38void loop() {
39  //-----------------------------------------
40  // Defitiniton of Grid Parameters
41  
42  dx = 1.0/(nxc-1.0);
43  x[1] = 0.0;
44  Rex[1] = 0.0;
45  for (int i = 2; i <= nx; i++) {
46    x[i] = x[i-1]+dx;
47    Rex[i]=Rel*x[i];
48  }
49
50  dy=fh/(nyc-1.0);
51  y[1]=0.0;
52   for (int i = 2; i <= ny; i++) {
53       y[i]=y[i-1]+dy;
54   }
55  //---------------------------------------
56
57  // Definition of Boundary Conditions
58  for (int i = 1; i <= nx; i++){
59    for (int j = 1; j <= ny; j++){
60      u[i][j] = 0.0;
61      v[i][j] = 0.0;
62    }
63  }
64
65  for (int j = 1; j <= ny; j++) {
66        u[1][j]=1.0; // Set horizontal Velocity at the grid with y=0 equal to 1 (inflow) 
67  }
68  
69  for (int i = 1; i <= nx; i++) {   // Set No slip boundary condition to the fluid/plate interface
70        v[i][1]=0.0;
71        u[i][1]=0.0;
72  }
73  //---------------------------------------
74
75  // Iterations to Solve the 2D Equations
76  float val;
77  for (int i = 2; i <= nx; i++) { // Iterating through all X/Y grid points
78    for (int j = 2; j <= nx-1; j++) {   // Calculate coefficients
79          CN[j] = -1.0/(Rel*(dy*dy)) - v[i-1][j]/(2*dy);
80          CP[j] = (u[i-1][j]/dx)+(2.0/(Rel*(dy*dy)));
81          CS[j] = (v[i-1][j]/(2*dy)) - (1.0/(Rel*((dy*dy))));
82          val = u[i-1][j];
83          Su[j] = (val*val)/dx;
84      }
85  
86      CP[1] = 1.0;
87      CS[1] = 0.0;
88      Su[1] = 0.0;
89      CN[ny] = 0.0;
90      CP[ny] = 1.0;
91      Su[ny] = 1.0;
92
93     // Solve Tridiagonal Problem
94     float Q;
95     for (int I = 2; I <= ny; I++){
96          Q = CN[I]/CP[I - 1];
97          CP[I] = CP[I] - CS[I - 1]*Q;
98          Su[I] = Su[I] - Su[I - 1]*Q;
99     }
100      
101      Q = Su[ny]/CP[ny];
102      u[i][ny] = Q;
103      for (int I = ny-1; I >=1; I--){
104        Q = (Su[I] - CS[I]*Q)/CP[I];
105        u[i][I] = Q;
106      }
107      // Calculate vertical velocity component
108      for (int j=2; j <= ny; j++){
109          v[i][j]=(dy/dx)*(u[i-1][j]-u[i][j]) + v[i][j-1];
110      }
111  }
112  // -----------------------------------------------------
113
114  // Perform Plotting 
115  p.Plot();
116
117  for (int kk = 1; kk <= nx; kk++){
118    for (int k = ny; k >= 1; k--){
119      uvel = u[kk][k];
120      xcord = y[k];
121      p.Plot();
122      delay(10);
123  }
124}
125}
126
127// END of code!