
[Sponsors] 
January 10, 2012, 11:31 
Manipulating source fvc::grad(vsf1+vsf2)

#1 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
Hello,
I need to manipulate the source term (fvc::grad(vsf1+vsf2)) of the momentum equation. The physical behaviour should be, that every identified internal face should be treated as a zero gradient condition: Code:
int icounterinternalfaces = 0; forAll(phi,phii) { label P = mesh.faceOwner()[phii]; label N = mesh.faceNeighbour()[phii]; if(vsf1[P]<(vsf2[N]vsf2[P])) { if(vsf1[N]<vsf3[N]){//faces which should be temporarily treated with a zero gradient condition for this timestep .... } } } I'm using the gauss gradient scheme ...==fvc::grad(vsf1+vsf2) First way I: Should I copy the gaussGrad.H/C and manipulate the code? But then I need to use something like ...==fvc::grad(vsf1,vsf2,vsf3) to be able to distinguish between normal internal faces and the ones which should have a zero gradient condition. So I need to create also a new gradSchemeclass. Second way: I use the fixedInternalValueFvPatchField, create temporarily an internal patch consisting of the specified faces and treating them with a zero gradient condition. But I'm not sure, if this is a good way. I would be happy about every comment on this. I tend to use the first way, or do you think, there is a more easy way like manipulating directly the outcome of the gradcalculation like: Code:
volVectorField sourceTerm = fvc::grad(vsf1+vsf2); forAll(phi,phii) { label P = mesh.faceOwner()[phii]; label N = mesh.faceNeighbour()[phii]; if(vsf1[P]<(vsf2[N]vsf2[P])) { if(vsf1[N]<vsf3[N]){//manipulating source term before giving it to the right sight of equation... sourceTerm...... = ..... } } } Andy 

January 26, 2012, 20:31 

#2 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
Try something like this in your code.
Code:
volScalarField sumVf = vf1+vf2; surfaceScalarField sumVff = fvc::interpolate(sumVf); volVectorField gradVsf = fvc::grad(sumVff); const unallocLabelList& owner = mesh.owner(); const unallocLabelList& neighbour = mesh.neighbour(); const vectorField& Sf = mesh.Sf(); const scalarField& V = mesh.V(); scalarField& sumVfi = sumVf.internalField(); scalarField& issf = sumVff.internalField(); vectorField& igGrad = gradVsf.internalField(); forAll(owner, facei) { label P = owner[facei]; label N = neighbour[facei]; if (vf1[P] < (vf2[N]vf2[P])) & (vf1[N] < vf3[N]) { igGrad[P] += Sf[facei]*(sumVfi[P]  issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N]  issf[facei]) / V[N]; } } 

February 3, 2012, 05:43 
what about the correctBoundaryConditions

#3 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
hi cliffoi,
many thanks for your hint. Yes, I use the gaussgrad scheme and in the meantime I had a deeper look in that source code. As you propose, I should calculate the gradient over all faces and subtracting the values of the specific faces, where neighour and owner fullfill some condition. Of course dividing the value with the volume of that cell. With volVectorField gradVsf = fvc::grad(sumVff); I use the first templateclass in gaussgradscheme, putting in the surfaceScalarField of the sumVf into the fvc::grad() and I get back the gradient (volVectorField gradVsf). But what about the correctBoundaryConditions(vsf, gGrad); In the original form, the gaussGrad does first the loop with the surfaceScalarField, and then the correctBoundaryConditions(vsf,gGrad) with looping over the patches (except the coupled ones). Of course, I have in the first templateclass the function: gGrad.correctBoundaryConditions(); But do I need the correctBoundaryConditions(vsf, gGrad); correction too? And if so, how can I perform this calculation/correction? Many thanks in advance! Best regards Andy 

February 3, 2012, 12:46 

#4 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
correctBoundaryConditions(vsf, gGrad) modifies the value of the gradient on all normal (noncoupled) boundaries so that the surfacenormal component of gradient matches the surfacenormal of the boundary condition. This correction has no influence on the internal field values of the gradient so, unless the boundary values are important to you, you can safely ignore this. Conveniently this function is static so you can call it yourself to correct the surfacenormal component afterwards.
Code:
volScalarField sumVf = vf1+vf2; surfaceScalarField sumVff = fvc::interpolate(sumVf); volVectorField gradVsf = fvc::grad(sumVff); const unallocLabelList& owner = mesh.owner(); const unallocLabelList& neighbour = mesh.neighbour(); const vectorField& Sf = mesh.Sf(); const scalarField& V = mesh.V(); // Apply correction to internalField scalarField& sumVfi = sumVf.internalField(); scalarField& issf = sumVff.internalField(); vectorField& igGrad = gradVsf.internalField(); forAll(owner, facei) { label P = owner[facei]; label N = neighbour[facei]; if (vf1[P] < (vf2[N]vf2[P])) & (vf1[N] < vf3[N]) { igGrad[P] += Sf[facei]*(sumVfi[P]  issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N]  issf[facei]) / V[N]; } } // Now apply correction to all coupled boundaries forAll(mesh.boundary(), patchi) { const unallocLabelList& pFaceCells = mesh.boundary()[patchi].faceCells(); const vectorField& pSf = mesh.Sf().boundaryField()[patchi]; const fvsPatchField<Type>& pssf = sumVff.boundaryField()[patchi]; if (sumVf.boundaryField()[patchi].coupled()) { const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField(); const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField(); const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField(); const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField(); // const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField(); const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField(); forAll(mesh.boundary()[patchi], facei) { if (vf1P[facei] < (vf2N[facei]vf2P[facei])) & (vf1N[facei] < vf3N[facei]) { label celli = pFaceCells[facei]; igGrad[celli]] += pSf[facei]*(sumVfi[celli]  pssf[facei]) / V[celli]; } } } } fv::gaussGrad<scalar>::correctBoundaryConditions(sumVf, gradVsf); 

February 7, 2012, 17:00 
summarising your last post

#5 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
Hello cliffoi,
thank you for your very fast reply. Please correct me, if I summaries something wrong: 1) First we go through the internal field, calculating the gradients with: volVectorField gradVsf = fvc::grad(sumVff); Then we manipulate the internal field with the first loop. Ok. In file gaussGrad.C, after dividing the cell gradient with the volume, then comes: Code:
gGrad.correctBoundaryConditions(); This, we have to manipulate again with your suggestion: Code:
// Now apply correction to all coupled boundaries forAll(mesh.boundary(), patchi) { const unallocLabelList& pFaceCells = mesh.boundary()[patchi].faceCells(); const vectorField& pSf = mesh.Sf().boundaryField()[patchi]; const fvsPatchField<Type>& pssf = sumVff.boundaryField()[patchi]; if (sumVf.boundaryField()[patchi].coupled()) { const scalarField vf1P = vf1.boundaryField()[patchi].patchInternalField(); const scalarField vf1N = vf1.boundaryField()[patchi].patchNeighbourField(); const scalarField vf2P = vf2.boundaryField()[patchi].patchInternalField(); const scalarField vf2N = vf2.boundaryField()[patchi].patchNeighbourField(); // const scalarField vf3P = vf3.boundaryField()[patchi].patchInternalField(); const scalarField vf3N = vf3.boundaryField()[patchi].patchNeighbourField(); forAll(mesh.boundary()[patchi], facei) { if (vf1P[facei] < (vf2N[facei]vf2P[facei])) & (vf1N[facei] < vf3N[facei]) { label celli = pFaceCells[facei]; igGrad[celli]] += pSf[facei]*(sumVfi[celli]  pssf[facei]) / V[celli]; } } } } Which is similar to your call: Code:
fv::gaussGrad<scalar>::correctBoundaryConditions(sumVf, gradVsf); Thank you very much. You are a great help! Best regards Andy Last edited by andyru; February 7, 2012 at 17:24. Reason: Better explanation 

February 7, 2012, 17:30 

#6 
Member
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 10 
That's pretty much it. Let me know if it actually works...


February 9, 2012, 15:14 
I made progess

#7 
Member
Andreas Ruopp
Join Date: Aug 2009
Location: Stuttgart / Germany
Posts: 31
Rep Power: 10 
Hello Cliffoi,
the implementation of the source worked and I could manipulate the source term of the momentum equation. Here I did one correction to your suggestion: Code:
... igGrad[P] += Sf[facei]*(sumVfi[P]  issf[facei]) / V[P]; igGrad[N] += Sf[facei]*(sumVfi[N] + issf[facei]) / V[N]; ... Then I got almost the expected result in my testcase. However, I manipulated the source term of the momentum equation before the pisoloop. So, I have a fluxfield, which is not divergence free, but approximately satisfies momentum (because this grad term has the same function as "fvc::grad(p)" in my incompressible flow). Now I have to do the same corrections in the pisoloop [having a correction (?) of the grad(p)/grad(vsf1+vsf2) values as I see in Ferziger & Peric] because I have fluctuations in the fluxvalues after the timestep ... So first, I must understand the pisoloop. I think, some new interesting questions will arise from my side soon. Sorry about that in advance... Best regards, Andy 

Tags 
fvc, gauss, grad, manipulate, source 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
swak4foam building problem  GGerber  OpenFOAM Installation  54  April 24, 2015 16:02 
pisoFoam compiling error with OF 1.7.1 on MAC OSX  Greg Givogue  OpenFOAM Programming & Development  3  March 4, 2011 18:18 
OpenFOAM on MinGW crosscompiler hosted on Linux  allenzhao  OpenFOAM Installation  127  January 30, 2009 20:08 
DxFoam reader update  hjasak  OpenFOAM PostProcessing  69  April 24, 2008 01:24 
DecomposePar links against liblamso0 with OpenMPI  jens_klostermann  OpenFOAM Bugs  11  June 28, 2007 17:51 