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 repres