
[Sponsors] 
January 21, 2008, 02:39 
File: backwardDdtScheme.c
Fun

#1 
New Member
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 10 
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() 

January 21, 2008, 04:20 
Yes you are quite right rho.bo

#2 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 15 
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 

January 21, 2008, 20:05 
Dr Henry Weller,
Thanks for

#3 
New Member
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 10 
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 bugfree in the future. Another thing. Could you kindly explain something I don't understand. See my post, http://www.cfdonline.com/cgibin/Op...cus/discus.cgi Liu Huafei 

January 23, 2008, 20:15 
Dr Henry Weller,
Thanks for y

#4 
New Member
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 10 
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 Φf =ΦC +(xf – xC) *ψ(rf) * (dΦ/ dx)f Or Φf =ΦC +(xf – xC) / (xD– xC) *ψ(rf) * (ΦDΦC) f ,the face C,central node D downwind node ψ(rf) the flux limiter rf is based on the following formula rf=(ΦCΦU)/( ΦDΦC) = 2*grad(ΦC)*(xDxC) /( ΦDΦC) 1 so for the faceflux mf mf>0 Φf =ΦP+fx*ψ(rf) * (ΦNΦP) =(1fx*ψ(rf)) ΦP +f x*ψ(rf)* ΦN mf<0 Φf =ΦN+(1fx)*ψ(rf) * (ΦPΦN) =(1fx) *ψ(rf)* ΦP +[1(1fx) *ψ(rf)]* * ΦN fx is the geometric interpolation factor, defined by fx=(xNxf)/(xNxP) 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ψ(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? 

January 24, 2008, 04:04 
I don't see the problem. The

#5 
Senior Member
Join Date: Mar 2009
Posts: 854
Rep Power: 15 
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 faceweight 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 

January 24, 2008, 21:07 
Thanks for your answer.

#6 
New Member
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 10 
Thanks for your answer.


January 27, 2008, 20:41 
If convection scheme is done b

#7 
New Member
Liu Huafei
Join Date: Mar 2009
Location: Shanghai, China
Posts: 20
Rep Power: 10 
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, Φf =max(mf,0)/mf * ΦP +min(mf,0)/mf * Φ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 Φf =ΦP+fx*ψ(rf) * (ΦNΦP) = (1fx*ψ(rf)) ΦP +f x*ψ(rf)* ΦN mf<0 Φf =ΦN+(1fx)*ψ(rf) * (ΦPΦN) = (1fx) *ψ(rf)* ΦP +[1(1fx) *ψ(rf)]* Φ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 

Thread Tools  
Display Modes  


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