|
[Sponsors] |
|
August 13, 2013, 08:10 |
|
#1 | |
New Member
Richard
Join Date: Aug 2013
Location: Zilina, Slovakia
Posts: 20
Rep Power: 12 |
Quote:
I am Ph.D. student, I would like use this program for simulations in my thesis. I wrote similar algorithm in 1D in Visual Basic, it works without any problem. OpenFOAM and C++ is new for me... Richard |
||
August 13, 2013, 16:45 |
Dont forget the chain rule
|
#2 |
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Germany
Posts: 213
Rep Power: 18 |
Hi,
when psat is a function of temperature, you are not allowed to pull it out of the gradient in the laplacian function. But you could use the chain rule to get an implicit term of F and an explicit term of psat. Here we go: Eqn (1): This is your starting point Code:
fvm::ddt(dw,F) + fvm::laplacian(delta,psat*F) Code:
fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat) + delta*psat*fvc::grad(F)) Code:
fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat)) + fvm::div(delta*psat*fvc::grad(F)) Code:
fvm::ddt(dw,F) + fvc::laplacian(delta*F,psat) + fvm::laplacian(delta*psat,F) Regards Fabian |
|
August 13, 2013, 18:02 |
|
#3 |
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,208
Rep Power: 26 |
Code:
dF/dt=laplacian(delta,P) P is partial pressure and could be rewrite F = P / Psat Psat is partial saturated pressure (function of Temperature) I tried this, but without succes... fvm::ddt(dw,F)=fvm::laplacian(delta*Psat,F) fvm::ddt(dw,F)=fvc::laplacian(delta,Psat*F) Code:
Eqn (1): This is your starting point Code: fvm::ddt(dw,F) + fvm::laplacian(delta,psat*F) Eqn (2): Split the gradient term Code: fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat) + delta*psat*fvc::grad(F)) Eqn (3): Split the divergence term Code: fvm::ddt(dw,F) + fvc::div(delta*F*fvc::grad(psat)) + fvm::div(delta*psat*fvc::grad(F))
__________________
Injustice Anywhere is a Threat for Justice Everywhere.Martin Luther King. To Be or Not To Be,Thats the Question! The Only Stupid Question Is the One that Goes Unasked. |
|
August 14, 2013, 02:43 |
|
#4 |
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Germany
Posts: 213
Rep Power: 18 |
Hi
The first laplacian Code:
+ fvc::laplacian(delta*F,psat) Code:
+ fvm::laplacian(delta*psat,F) Regards Fabian |
|
August 14, 2013, 05:07 |
|
#5 | |
Senior Member
Ehsan
Join Date: Oct 2012
Location: Iran
Posts: 2,208
Rep Power: 26 |
Quote:
__________________
Injustice Anywhere is a Threat for Justice Everywhere.Martin Luther King. To Be or Not To Be,Thats the Question! The Only Stupid Question Is the One that Goes Unasked. |
||
August 14, 2013, 09:56 |
|
#6 | |
New Member
Richard
Join Date: Aug 2013
Location: Zilina, Slovakia
Posts: 20
Rep Power: 12 |
Quote:
Thank you very much, I am going to try it. Yesterday I used other approach (similar which I using in Visual Basic) for implicit scheme and I think that results was correct. - generate fvscalarmatrix "mat" from laplacian(delta,F) - multiply each cell in fvscalarmatrix mat by Psat in point of this cell (using lduAddressing) - solve (fvm::ddt(dw,F)==mat) But this code is more complicated and your approach seems better. dw is aproximation function of sorption curve (relation between water content and relative humidity of air) it is nonlinear function, I forgot write in equation Richard |
||
August 14, 2013, 10:47 |
|
#7 |
New Member
Richard
Join Date: Aug 2013
Location: Zilina, Slovakia
Posts: 20
Rep Power: 12 |
The part "fvc::laplacian(delta*F,Psat)" causes stability problem,...
volScalarField F internalField nonuniform List<scalar> 4000 ( 0.503748 0.500001 0.5 0.5 0.5 0.5 0.49999 0.44179 -361.167 -2.31903e+06 0.500001 0.5 0.5 0.5 0.5 0.5 0.49999 0.44179 -361.165 -2.31902e+06 |
|
August 14, 2013, 12:22 |
|
#8 |
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 |
you could try to implement it implicitly
fvc::laplacian(delta*F,Psat) = fvc::div( delta*F*grad(Psat) ) volVectorField dps = delta * fvc::grad(Psat); surfaceScalarField dpsf = fvc::interpolate(dps) & mesh.Sf(); then you can view the dps field as a 'convective' change of F due to the Psat gradient and rewrite it like this fvc::laplacian(delta*F,Psat) = fvm::div(dpsf, F) ...should work |
|
August 14, 2013, 13:35 |
|
#9 |
Senior Member
Fabian Roesler
Join Date: Mar 2009
Location: Germany
Posts: 213
Rep Power: 18 |
Hi
That's a pitty. But the code of niklas seems to be promissing. Ricardo, I didn't get your approach in total. Could you explain it in a little more detail? Especially the part of the lduAddressing. Regards Fabian |
|
August 14, 2013, 14:21 |
|
#10 | |
New Member
Richard
Join Date: Aug 2013
Location: Zilina, Slovakia
Posts: 20
Rep Power: 12 |
implementation to solver >> solve (fvm::ddt(dw,F)==fvm::laplacian(lambdaD*Psat,F)+fv m::div(dpsf, F)); doen't works, too
Quote:
a part of code ... Code:
float tdelta = 1.0; float tdeltaN = 1.0; float fdelta = 1.0; float fdeltaN = 1.0; double relFaktor = 1.0; int iter = 0; #include "boundaryCalculation.H" do { Tt = T; #include "heatParameters.H" solve ( fvm::ddt(cap,T)==fvm::laplacian(lambda, T)//+hv*fvc::laplacian(lambdaD, Psat*F) ); T=(relFaktor)*T+(1-relFaktor)*Tt; tdeltaN=diffFields(T,Tt); #include "saturatedpressure.H" #include "moistureParameters.H" fvScalarMatrix mat1 ( fvm::laplacian(lambdaD,F) ); label Nc = mesh.nCells(); //Total number of cells // Diagonal coefficients calculation for(label i=0; i<Nc; i++) { mat1.diag()[i]*=Psat[i]; } //Off-diagonal coefficients calculation for(label faceI=0; faceI<mat1.lduAddr().lowerAddr().size(); faceI++) { mat1.lower()[faceI]*=Psat[mesh.faceOwner()[faceI]]; mat1.upper()[faceI]*=Psat[mesh.faceNeighbour()[faceI]]; } dimensionedScalar Pa ("Pa",dimensionSet(1,-1,-2,0,0,0,0), 1.0 ); mat1*=Pa; solve ( fvm::ddt(dw,F)==mat1 //fvm::laplacian(lambdaD*Psat,F)//+fvc::laplacian(dT,T)// // ); F=(relFaktor)*F + (1-relFaktor)*Ft; fdeltaN=diffFields(F,Ft); #include "w_prepocet.H" iter++; if (iter==1) { tdelta=tdeltaN; fdelta=fdeltaN; } if ((tdeltaN>tdelta) || (fdeltaN>fdelta)) { relFaktor=relFaktor/2; T=(relFaktor)*T+(1-relFaktor)*Tt; F=(relFaktor)*F + (1-relFaktor)*Ft; } } while (((tdeltaN>maxDeltaT) || (fdeltaN>maxDeltaH)) && iter < maxIter); I tried next approach, rewrite partial pressure by gradient of temperature, it seems, that gives similar values like lduaddressing script I don't know where is problem, I have to check all functions, material parameters and other properties again,... thanks to all for support and advices Richard |
||
August 14, 2013, 15:11 |
|
#11 |
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29 |
you REALLY need to be more specific when you say 'it doesnt work'.
|
|
August 14, 2013, 15:39 |
|
#12 |
New Member
Richard
Join Date: Aug 2013
Location: Zilina, Slovakia
Posts: 20
Rep Power: 12 |
I was wrong , It works, but obtained values ...
dimensions [0 0 0 0 0 0 0]; internalField nonuniform List<scalar> 4000 ( -1.35814e-10 4.58031e-11 4.58031e-11 4.58031e-11 4.5803e-11 4.58031e-11 4.58031e-11 4.58031e-11 -1.82592e-10 -1.77054e-12 -2.23032e-11 So it seems (two scheme unstable), that a problem could be somewhere else, Is right using parameters (lambda), stored as volScalaField and (ununiform field) calculated in each step equation by this way? (fvm::laplacian(lambda,T)) Is it interpreted correctly to scheme or is used average values? Thanks Richard |
|
April 14, 2020, 01:13 |
|
#13 | |
Senior Member
TWB
Join Date: Mar 2009
Posts: 402
Rep Power: 19 |
Quote:
I face similar problem whereby I need to convert Code:
fvm::laplacian((alpha*rho*DkEff(F1)), k_) Code:
fvm::laplacian((alpha*DkEff(F1)), rho*k_) Code:
fvm::laplacian((alpha*rho*DkEff(F1)), k_) + fvc::laplacian((alpha*k_*DkEff(F1)), rho) Code:
fvc::laplacian((alpha*k_*DkEff(F1)), rho) Another error is "note: template argument deduction/substitution failed" Any idea how to solver it? Must been trying to solve it for a month. Thanks |
||
|
|