[liuhuafei]

File: backwardDdtScheme.c
Function:
template<class>
tmp<geometricfield<type,> >
backwardDdtScheme<type>::fvcDdt
(
const volScalarField& rho,
const GeometricField<type,>& vf
)
{
......

rDeltaT.value()*( coefft*vf.boundaryField() -
(coefft0*rho.oldTime().boundaryField() *vf.oldTime().boundaryField()
- coefft00*rho.oldTime().oldTime().boundaryField()
*vf.oldTime().oldTime().boundaryField()
)
)
......

coefft*vf.boundaryField()
where is rho.boundaryField()

[henry]

Yes you are quite right rho.boundaryField() should be in coefft*vf.boundaryField() , i.e.

coefft*rho.boundaryField()*vf.boundaryField()

Thanks for reporting it

Henry

[liuhuafei]

Dr Henry Weller,
Thanks for your answer. I am checking the source code. I am happy that I will find a bug in the foam. Maybe the foam will be bug-free in the future.

Another thing. Could you kindly explain something I don't understand. See my post,

http://www.cfd-online.com/cgi-bin/Op...cus/discus.cgi

Liu Huafei

[liuhuafei]

Dr Henry Weller,
Thanks for your answer. Could you kindly explain something I don't understand.
post,

Some questions about the implementation of the convection scheme
--------------------------------------------------------------------------------
By Liu Huafei on Sunday, January 20, 2008 - 05:18 pm: Edit Post
It seems that something wrong with the implementation of the convection scheme. I guess the face value is based on the following formula

&Phi;f =&Phi;C +(xf &ndash; xC) *&psi;(rf) * (d&Phi;/ dx)f
Or
&Phi;f =&Phi;C +(xf &ndash; xC) / (xD&ndash; xC) *&psi;(rf) * (&Phi;D-&Phi;C)

f ,the face
C,central node
D downwind node
&psi;(rf) the flux limiter

rf is based on the following formula
rf=(&Phi;C-&Phi;U)/( &Phi;D-&Phi;C) = 2*grad(&Phi;C)*(xD-xC) /( &Phi;D-&Phi;C) -1

so
for the faceflux mf
mf>0 &Phi;f =&Phi;P+fx*&psi;(rf) * (&Phi;N-&Phi;P)
=(1-fx*&psi;(rf)) &Phi;P +f x*&psi;(rf)* &Phi;N
mf<0 &Phi;f =&Phi;N+(1-fx)*&psi;(rf) * (&Phi;P-&Phi;N)
=(1-fx) *&psi;(rf)* &Phi;P +[1-(1-fx) *&psi;(rf)]* * &Phi;N

fx is the geometric interpolation factor, defined by fx=(xN-xf)/(xN-xP)
P is owner cell
N is neighbour cell

but in the limitedSurfaceInterpolationScheme.H and limitedSurfaceInterpolationScheme.c
template<class>
tmp<surfacescalarfield> limitedSurfaceInterpolationScheme<type>::weights
(
const GeometricField<type,>& phi
) const
{
const fvMesh& mesh = this->mesh();
tmp<surfacescalarfield> tWeights(this->limiter(phi)); this is&psi;(rf)!!!!!!!
surfaceScalarField& Weights = tWeights();
const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights(); fx !!!! scalarField& pWeights = Weights.internalField();
forAll(pWeights, face){
pWeights[face] = pWeights[face]*CDweights[face]
+ (1.0 - pWeights[face])*pos(faceFlux_[face]); It is right??????
}
Maybe openfoam is based on other formula. Could you kindly explain clearly it?

[henry]

I don't see the problem. The only thing in the code which may be confusing is the reuse of Weights, it is first used to store the limiter field and then becomes the weights field because after each face-weight is calculated the corresponding limiter value is no longer needed. I can add a comment on this if you think it would be useful.

Henry

[liuhuafei]

Thanks for your answer.

[liuhuafei]

If convection scheme is done by such way,some contradiction arises.

LDUMatrix contains the coefficients, lowerPtr,upperPtr,diagPtr ,and
the address reference, uPtr and lPtr. The upperPtr stores the contribution
of the neighbour cell(N) to the owner cell(P) and the lowerPtr store
the contribution of the owner(P) to the neighbour(N).

For upwind scheme, &Phi;f =max(mf,0)/mf * &Phi;P +min(mf,0)/mf * &Phi;N, so in upwind.h
the weight factor is calculated by
tmp<surfacescalarfield> weights() const
{
return pos(this->faceFlux_);
}.

It is OK. From this code,I guess the weight is corresponding to the following code

fvm.lower() = -weights.internalField()*faceFlux.internalField();
fvm.upper() = fvm.lower() + faceFlux.internalField();
fvm.negSumDiag();

fvm.lower is the contribution of the owner(P) to the neighbour(N),fvm.upper
the contribution of the neighbour(N) to the owner(P).

But for the high order convection scheme,
mf>0 &Phi;f =&Phi;P+fx*&psi;(rf) * (&Phi;N-&Phi;P) = (1-fx*&psi;(rf)) &Phi;P +f x*&psi;(rf)* &Phi;N
mf<0 &Phi;f =&Phi;N+(1-fx)*&psi;(rf) * (&Phi;P-&Phi;N) = (1-fx) *&psi;(rf)* &Phi;P +[1-(1-fx) *&psi;(rf)]* &Phi;N

the weight factor is calculated by
pWeights[face] = pWeights[face]*CDweights[face]
+ (1.0 - pWeights[face])*pos(faceFlux_[face]);

It is right? It is corresponding to the fvm.upper rather than fvm.lower,why?
fvm.lower is not the contribution of P to N?

Another question: Diag[i] should be substract the contribution from
all its neighbouring cell.
In fvMatrix, the implementation of function negSumDiag() is as follows
negSumDiag:
for (register label face=0; face<l.size(); face++)
{
Diag[l[face]] -= Lower[face]; // AP[IP]=AP[IP]-AL[IFF] is correct??????
Diag[u[face]] -= Upper[face]; // AP[IN]=AP[IN]-AU[IFF] is correct??????
}

why not the following code?
for (register label face=0; face<l.size(); face++)
{
Diag[l[face]] -= Upper[face]; (Upper is the contribution of N to P)
Diag[u[face]] -= Lower[face]; (Lower is the contribution of P to N)
}

Liu Huafei