CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Bugs

Bug of backwardDdtScheme

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   January 21, 2008, 02:39
Default File: backwardDdtScheme.c Fun
  #1
New Member
 
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 8
liuhuafei is on a distinguished road
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()
liuhuafei is offline   Reply With Quote

Old   January 21, 2008, 04:20
Default Yes you are quite right rho.bo
  #2
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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
henry is offline   Reply With Quote

Old   January 21, 2008, 20:05
Default Dr Henry Weller, Thanks for
  #3
New Member
 
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 8
liuhuafei is on a distinguished road
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 is offline   Reply With Quote

Old   January 23, 2008, 20:15
Default Dr Henry Weller, Thanks for y
  #4
New Member
 
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 8
liuhuafei is on a distinguished road
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?
liuhuafei is offline   Reply With Quote

Old   January 24, 2008, 04:04
Default I don't see the problem. The
  #5
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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
henry is offline   Reply With Quote

Old   January 24, 2008, 21:07
Default Thanks for your answer.
  #6
New Member
 
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 8
liuhuafei is on a distinguished road
Thanks for your answer.
liuhuafei is offline   Reply With Quote

Old   January 27, 2008, 20:41
Default If convection scheme is done b
  #7
New Member
 
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 8
liuhuafei is on a distinguished road
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
liuhuafei is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Serious bug in backwardDdtScheme henry OpenFOAM Bugs 1 June 2, 2007 05:37


All times are GMT -4. The time now is 17:43.