\label{Intro}

The laws governing fluid flow around an obstruction can be modeled after a set of two coupled partial differential equations for \(\psi\) and \(\xi\): \[\nabla^{2}\psi = \xi\] \[\nabla^{2}\xi = \frac{1}{\nu}(\frac{\partial \psi}{\partial y}\frac{\partial \xi}{\partial x} - \frac{\partial \psi}{\partial x}\frac{\partial \xi}{\partial y})\] \(\psi\), or the stream function, represents the path of the fluid when looking at its contour lines. \(\xi\), or the vorticity, represents the curl of the fluid flow, or how much the fluid curves while travelling around this obstruction.

In this project, we are interested in computationally modeling the path of a fluid around an obstruction using a two dimensional grid of points. We assume our fluid flowing is symmetric around the obstruction and thus we only consider the top half of our boundary. By modeling our fluid flow on a grid, we are able to use relaxation methods, the general idea of which is to split the matrix into easily invertible and harder to invert points. We can then construct the solution iteratively by using the averages of the four nearest neighbors in order to compute the solution of a grid point in the next iteration. We use a linear combination of updated values using a factor \(\omega\) which we set to values between 1 and 2. in order to anticipate future corrections. This variant of the Gauss-Sidel method, which allows our solution to converge faster, is called successive overrelaxation.

We first model our fluid so that it is in a free flow pattern. We then use a combination of implementation of our boundary conditions and overrelaxation to calculate \(\psi\) and \(\xi\) values iteratively. This process is repeated until our residuals, which we also compute each time we sweep through the psi and xi matrices, are reasonably small.

By plotting the contour lines of \(\psi\), we will be able to see how the fluid flows around the obstruction. We can also plot \(\xi\) to see how the fluid swirls around the obstruction.

In addition to looking at just how the fluid flows around the obstruction, we are also interested in the parameters that change this fluid flow, particularly its viscosity and initial velocity, encompassed in the Reynolds number. Also, we are interested in how moving the obstruction closer to the edges of the boundaries changes the fluid flow pattern.

\label{Methods}

We first create a structure, called “fluid”, containing matrices for \(\psi\), \(\xi\) and residual values for both \(\xi\) and \(\psi\). We initialize values for our psi and xi matrices to model free flow as though there is no obstruction.

The sweep method first updates \(\psi\) on the interior of its matrix, including the surface of the plate. To do this, we use the method of overrelaxation, or a linear combination of the nearest neighbors using centered differencing. However, along the plate, because we have a viscous fluid, the particles ‘stick’ to the plate, so we set the \(\psi\) values equal to zero.

Once the interior has been updated, we use these values to implement the exterior boundaries of our grid. To implement our boundary conditions, we iterate through each point along the boundary of the matrix. For the left boundary, the fluid has not been affected by the obstruction, and therefore has only the x component of the velocity equal to a constant. Therefore values on the boundary are equal to the values directly to the right of the boundary. The same boundary condition applies to the upper boundary, and therefore the values on this boundary are equal to the values directly below it. On the right boundary, nothing changes due to the obstruction, therefore the values on this boundary are equal to the values to the left of them. Due to the viscosity of the fluid, the \(\psi\) values on the edges of the plate are 0. On the streamlines where our line of symmetry lies, we choose \(\psi\) to be equal to zero.

Because the \(\psi\) and \(\xi\) values are coupled, we use these updated \(\psi\) values in order to calculate \(\xi\). We first update \(\xi\) values along the left, top, and right side of the plate using our boundary conditions. We can then use those updated values to update \(\xi\) on the interior using overrelaxation. Finally, we update \(\xi\) on the right boundary. The other boundaries are not affected by the plate and thus do not have any vorticity, and we can leave the \(\xi\) values as they are. For the right boundary, however, we know that the vorticity is not changing in the horizontal direction, so we can set values along the right-hand boundary equal to the points to their immediate left.

Finally, we compute residuals for both our \(\xi\) and \(\psi\) matrices. This represents the error from the successive overrelaxation method. For both \(\psi\) and \(\xi\), we calculate the norm of these residuals. We then loop over this whole process of successive overrelaxation and implementation of boundary conditions until the norm of \(\psi\) is below a certain a reasonably low value, which we determined to be \(10^{-13}\), around 9000 sweeps. We used the norm of the residual of \(\psi\) instead of the norm of the residual of \(\xi\) because it was higher at higher numbers of sweeps and therefore more indicative of when we had reached a suitable number of sweeps.

As we iterated through our grid of points representing \(\psi\) and \(\xi\) using our while loop, we saw that the contour lines of our \(\psi\) function began to model what we would expect of a fluid going around an obstruction: namely that the fluid did in fact flow around the obstruction. While the contour values are now shown on our graphs the \(\psi\) contour values went from 0 along the boundary of the plate (because the viscous fluid does not flow along the edges of the plate but instead gets stuck), to 0.6131 at the edge of the boundary, theoretically not affected by the plate obstruction. The shape of the contours clearly demonstrates that our fluid flows around the obstruction, and can be seen in our graphs (Figure 1). Additionally, we can see a very weak eddy forming directly behind the plate with a \(\psi\) contour value of 0.01. Similarly, our vorticity contours clearly show that there is no curl of the vector field far enough away from the obstruction, as we would expect, but very high values around the boundaries of the plate, where the fluid is being diverted away (Figure 2). We can also see this in the surface plot of \(\xi\) (Figure 3).

For all of our plots, in order to compare them to each other, the contour spacing is the same. While the Some of the contours plots look very black, which is partially due to larger ranges on contour values as well as large changes in \(\psi\) or \(\xi\).

We continued executing “sweeps” until the norm of the residual of \(\psi\) reached a value of \(10^{-13}\), at which point the model clearly showed the fluid flowing around the obstruction. We plotted the residual of both \(\psi\) and \(\xi\) on a surface plot over the grid and saw that for low tolerances such as \(10^{-13}\), the norm for \(\psi\) was 0 everywhere (Figure 4). The \(\xi\) residual grid was also zero at most places. However, close to the surface of the plate, where the fluid is swirling much more because it is interacting with the obstruction, the residuals were notably higher (Figure 5). We can explain this by looking at the surface plot of \(\xi\), (Figure 3) where we can see that there is a high vorticity around the plates, then a sharp drop off when inside of the plate. Because our residuals are calculated using neighboring points, this sharp drop off causes some computational error which results in higher residuals (though still relatively small).

When we saw the norm of the residuals was 0 almost everywhere, we were a little alarmed, leading us to look at the norms of the residuals with a higher minimum norm to terminate the program. This would allow us to look at our residuals with fewer sweeps. When we set the norm of the residual of \(\psi\) to \(10^{-4}\), we saw higher residuals of \(\xi\) near the boundary of the plate and higher residuals of \(\psi\) throughout the grid. This is because our values have not run through enough sweeps for the successive overrelaxation method to be completely accurate, which validated our code because we knew the process of overrelaxation was occuring.

When initial parameters are set to regular suggested values and we do enough sweeps so that the norm converges to zero, we were interested in seeing how our relaxation term \(\omega\) affected the accuracy of our code. We plotted the norm versus the number of sweeps (Figure 6), then did a semi-log plot (Figure 7) in order to see the resolution between the \(\omega\) values more carefully. We noticed a general pattern that for \(\omega\) values of 1.5 or less, the norm converges to zero very quickly, but for values of 1.6 to 1.64, the norm takes longer to converge. The initial drop off in the graph comes from when we arbitrarily set the norm to a large number (100) so that our while loop would initiate. We also found that the residual of the norm of \(\psi\) is smallest around 1.64 but at values above 1.65, the norm diverges, indicating that this is a tipping point.

We were also interested in how the number of grid points affects our model. We varied the resolution of the grid from 16 to 128 cells on a side and looked at the affects on \(\psi\). We observed how the contour lines of \(\psi\) are better defined on the plate and appear smoother as they arch around the plate with more grid cells. We also noticed that the contour lines of \(\xi\) look similar in that the vorticity is zero nearly everywhere except around the surface of the plate, where it is quite large. However with 128 grid cells there is much more resolution around the surface of the plate, to the point that you can even see the outline of the obstruction (Figures 10, 11).

Additionally, we played around with moving the plate and changing its size to see how that affected our model. We first changed the boundaries of the plate such that top of the plate was halfway up the grid at y value of 0.5 (stretching it vertically). Initially, the residual of the norm of \(\psi\) did not converge, indicating that our plate was too close to our boundaries. We then doubled the viscosity to 0.2 and found that the norm converged. The contour graph of \(\psi\) looks as one might predict with stream lines around the plate that flatten out after the plate (Figure 12). This flattening (instead of coming down the right side of the plate) might also be due to the high viscosity. We also plotted the contour of \(\xi\), and noticed it had a similar pattern as the graphs prior (Figure 13) Additionally, we moved the right end of the plate closer (stretching it horizontally) to the right boundary at an x value of 0.75. We see relatively flat streamlines over the plate that continue to the edge of the grid (Figure 14). The \(\xi\) contour lines also appear to swirl in a similar pattern to the graphs above (Figure 15).

Finally, we increased the velocity of the fluid from 1 to 10 to see how the initial velocity changes the fluid flow. Initially the norm did not converge. We used the final \(\psi\) and \(\xi\) grids with a velocity of 1 to then sweep through the grid at a velocity of 10, which allowed the norm to converge. We observe streamline contours behind the plate that circle around resembling eddies (Figure 16). Because the fluid is flowing so quickly, the \(\psi\) contours are close together and appear black. When the initial velocity was at its max of 10, our contour plot of vorticty showed a lot of curl, making the entire plot look black. To rectify this, when the contours were the closest together, we chopped off the top in order to see the full resolution of the eddies (Figure 17). Our surface plot of \(\xi\) shows the full range of \(\xi\)(Figure 18).