[JBUNSW]

Dear Foamers,

I am seeing some weird behaviors from fvc::grad(phi) operator under periodic boundary conditions.

To make sure that fvc::grad(phi) is working correctly, I setup a 1D domain with length 2mm divided into 400 equal length cells. so

delX = 2e-3 / 400

phi is a volScalarField obtained solving a simple transport equation. You can set some random numbers to build phi, doesn't matter actually.

Then I printed out the fvc::grad(phi) along with the phi at all 400 cells. Since the problem is 1D we don't have gradients in y or z direction, therefore, the gradient at internal points can be simply calculated by a central difference scheme as

grad(phi_i) = ( phi_i+1 - phi_i-1 ) / (2 * delX)

I copied the data into a spreadsheet and did the above. Below are the results I'm getting:

Cell 	   phi	        mag(grad(phi))	     abs(d(phi)/dx) 	Err %
				 (OpenFOAM)	      (manual calc)
----------------------------------------------------------------------------
1	244676	6.31092e-317			-		     -	
2	246466	6.38511e-317		358100000		#DIV/0!
3	248257	0				358300000		#DIV/0!
4	250049	0				358400000		#DIV/0!
5	251841	3.58454e+08		358400000		0.015
6	253633	3.5851e+08			358500000		0.003
7	255426	3.58566e+08		358600000		-0.010
8	257219	3.58622e+08		358600000		0.006
9	259012	3.58564e+08		358600000		-0.010
10	260805	3.58506e+08		358500000		0.002
11	262597	3.58447e+08		358400000		0.013
12	264389	3.58503e+08		358500000		0.001
13	266182	3.58444e+08		358500000		-0.016
14	267974	3.58271e+08		358300000		-0.008
15	269765	3.58098e+08		358100000		-0.001
16	271555	3.57924e+08			-		     -	
.	.		.					.		     .
.	.		.					.		     .
.	.		.					.		     .
392	259012	3.58564e+08			-		     -	
393	257219	3.58622e+08		358600000		0.006
394	255426	3.58566e+08		358600000		-0.009
395	253633	3.5851e+08			358500000		0.003
396	251841	3.58454e+08		358400000		0.015
397	250049	3.58398e+08		358400000		-0.001
398	248257	3.58228e+08		358300000		-0.020
399	246466	3.58058e+08		358100000		-0.012
400	244676	1.79043e+08			-		    -

I am using periodic boundary condition in OpenFOAM 2.0.x. Cell 1 is residing at x=0 and cell 400 ends at x = 2mm. As you can see for all internal nodes safely away from the boundaries the error is negligible. But when you get close to the boundary nodes the periodic boundary condition is affecting the way gradient of phi is being calculated inside the domain. The interesting thing is that it is only happening in x=0 end of the domain and the other end apparently is fine.

I have found some similar threads regarding fvc::grad, for instance:
Gradient Operator
Or some threads on periodic boundary condition in general, for example:
How to setup cyclic BCs in simpleFOAM
Or threads like this on gradient at boundaries in general (read the last post):
Gradient at the boundary



My questions:

1. Am I right?! Is this a bug?!

2. Any body knows how to fix this?! Any thoughts?! Because I need to use the periodic boundary condition in a 2D turbulent combustion context, so I need a reliable gradient operator, otherwise I have to hardcode my own implementation of gradient operator, which is not an elegant solution.

3. Speaking of manually coding my custom version of gradient operator, does anybody know how the gradient is calculated at the boundaries?!

4. How the periodic boundary condition is implemented in OpenFOAM?! This wikipedia article gives an idea, but I want to know how and where it is implemented in OpenFOAM. Anybody knows a good starting point?! Which class?!



Appreciate your time for reading this and perhaps trying to answer either of the above...

Jalal