I've recently got interested in CFD and started a small project by solving the radial Reynoldsequation. Why the Reynoldsequation? I recently encountered it through my studies and somehow got stuck :)
I would greatly appreciate it, if someone were so kind to have a read through my steps and point out my thought mistakes or any other errors.
In a first approach, I would like to solve it without the source term, which reduces it basically to a two dimensional diffusive process.
$\frac{\partial}{\partial r}(rh^3\frac{\partial p}{\partial r}) + \frac{1}{r}\frac{\partial}{\partial \theta}(h^3\frac{\partial p}{\partial \theta}) = 0$
If one were to imaging a ring, with two Diriclet Boundary Conditions: $p(r_I) = 2$ and $p(r_O) = 250$, so I don't need to worry about the singularity at $r=0$. I also have a periodic condition along the angle at $\theta = 0$ and $\theta = 2\pi$ The actual pressure values don't matter, as long as they are larger as zero. Negative pressures don't really work :D In the end, I shall have a unsymmetrical problem with several boundary conditions on the ring, that's why I don't wish to reduce this to a one dimensional one.
Following the approach by Patanker or any other book on Finite Volume Methods (such as The Finite Volume Method in Computational Fluid Dynamics by F. Moukalled), the semi discretization with the midpoint rule will give me a second order accurate approximation:
$r_eh_e^3\Delta \theta(\frac{\partial p}{\partial r})_e - r_wh_w^3\Delta \theta(\frac{\partial p}{\partial r})_w + \frac{\Delta r}{r_i}(h_s^3(\frac{\partial p}{\partial \theta})_s - h_n^3(\frac{\partial p}{\partial \theta})_n)$
That can be confirmed using the Taylor analysis for example.
For convenience, I wrote the code in Matlab, so the East direction is $j+1$ and the South direction is $i+1$ for a matrix given by $M(i,j)$ Lowerbound letters, such as e and s denote the faces midway between two points.
I use a structured grid and equidistant in r and $\theta$ ($\Delta r \neq \Delta \theta$)
Now I'm getting to the part, where I've been scratching my head for some time now.
If one were to approximate the gradient at the faces of the control volume like: $(\frac{\partial p}{\partial r})_e = \frac{p_E - p_i}{\Delta r}$ and $(\frac{\partial p}{\partial r})_w = \frac{p_i - p_W}{\Delta r}$ and similar for the south and north faces, the formal order of accuracy is what?
Taylor analysis yields a formal order of accuracy of $1$ for the gradient approximation as above, while my books and other sources (internet etc.) refer to the combined discretization (see below) as of order $2$ on a uniform grid and others as of order $1$. They contradict themselves, which is kind of annoying.
$\frac{r_eh_e^3(p_E-p_i)\Delta \theta}{\Delta r} - \frac{r_wh_w^3(p_i - p_W)\Delta \theta}{\Delta r}$
I understand, that when I choose $d_1 = r_eh_e^3\Delta \theta$ and $d_2 = r_wh_w^3\Delta \theta$ (similar for South and North) and in the special case of $d_1 = d_2$, the above discretization would become a second order central difference scheme look a like. But I don't have that case.
$\frac{d (p_E - 2p_i + p_W)}{\Delta r}$
I implemented the discretization and performed an order of accuracy study yielding an order of accuracy of $1$ by comparing the solution of the scheme in radial coordinates to the exact solution. I did that by integrating the resulting pressure field to obtain the hydrostatic lift force. I can now compare the numerical one to the exact one.
Is this actually allowed? The order of accuracy I get is basically 1 (+/- 0.003) for fine grids and starts to oscillate for coarser grids as expected. Would it be better to directly compare the resultant pressure fields?
That's basically it. I sifted through so many texts that I'm now confused to a point, where every explanation I try to come up with is clouded in self doubt.
I would like to find out, which formal order of accuracy I should obtain and how I can adapt the discretization to achieve an order of accuracy of $2$.
The above discretization's of the gradients are basically a linear profile between two points, but I'm not sure how a higher order profile would look like.
Edit 1:
I use a staggered grid, where the pressure is stored at the centre of the cell and the height at the cell faces. The radii, for example $r_e$, are not approximated but direclty computed at the position of the faces by $(j+c \pm 0.5)*\Delta r$, where c is a constant to account for the inner radius, as my domain does not extend to $r=0$ but starts at $r_I$.
To be able to use the exact solution, the height is currently uniform and could be taken out of the equation, but will be later evaluted through the interpolation of a matrix. Additionally, the density and viscosity could also be added, currently assumed constant, and would also be evaluated at the cell faces.
I also updated the above equations, such that capital letters denote a cell centre and the small letters the cell faces, where i is also at the centre of a cell and conincides with the centre of my stencil.
Edit 2:
I'm currently using these three books:
1: Patanker: Numerical Heat Transfer and Fluid Flow Link
2: F. Moukalled: The Finite Volume Method in Computational Fluid Dynamics Link
3: J.H. Ferziger: Computational Methods for Fluid Dynamics Link
Patanker doesn't go much into detail regarding the order of accuracy, while in the second book they note (page: 216), that the semi discretization is of order two, but the linear profile used to approximate the gradient is of first order, which would make the complete approximation be of order $1$ from my point of view.
The third book goes more into detail regarding the linear profile used for the gradient on pages: 72-74. They note, that the linear gradient approximation is of formal order $1$, but behaves as a second order approximation on a uniform grid, when one chooses the faces midway between two neighbouring points, as I have. On Page 75 it states, that the scheme converges asymptotically in a second order manner.
Several lecture notes on the internet just mentioned, it would be second order accurate without going much into detail.
\enclose{horizontalstrike}{After going through the texts again, from my understanding, I should obtain a order of accuracy of one, which I do. Yet, I would like to be certain of it^^}
Answer May I answer my own question? Well, anyway.
I've now been able to obtain an order of accuracy of $2$. I had the error norm computed wrong or in other words, I did not scale the error of the obtained lift force to the grid.
