Stationary Fluid Flow

Henry Daniels-Koch and Holly Rudel



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.



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