
[Sponsors] 
August 13, 2015, 23:33 
fvc::div( )==0

#1 
New Member
Gergely Acs
Join Date: Jun 2015
Location: Hong Kong
Posts: 4
Rep Power: 9 
Hello Foamers,
I'm developing a solver for electrochemical flow simulations, and I have to solve an equation, which has the general form (similar to the incompressible continuity equation): div(flux)=0, which I'll to solve as fvScalarMatrix UeEqn ( fvc::div(Flux) ); where Flux is a volVectorField. However, when I try to compile my solver, I always get error messages, I've tried also some other ways to code it, but without success so far. I would be glad if you could help me with this problem. Regards, Gergely 

August 14, 2015, 09:38 

#2 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
You will likely need a reformulation such as div(grad(phi))=0 and then calculate flux=grad(phi) or so.
An example for this is magnetostatics with a constant current, phi is electrostatic potential, and what you want to solve is div(j)=0 with j being current density (j=sigma*grad(phi)). What you end up solving is laplacian(sigma, phi)=0, i.e. div(sigma*grad(phi))=0. Also keep in mind that div(flux)=0 doesn't have a unique solution. 

August 17, 2015, 04:20 

#3 
New Member
Gergely Acs
Join Date: Jun 2015
Location: Hong Kong
Posts: 4
Rep Power: 9 
Thank you, chriss85, indeed expanding the equations did help, my solver could be compiled, however, it lead me only to "incompatible fields", thus I expanded my equation even deeper.
My eqaution is div( (D1D2)*grad(C)+(D1+D2)*F/R/T*C*grad(Ue) ). Where D1,D2,F,R & T are dimensioned scalar constants, C and Ue are my variables, both are scalar fields. After expansion I've found: laplacian( (D1D2)*grad(C))+ (D1+D2)*F/R/T* ( grad(C)*grad(Ue)+C*laplacian(Ue) ). Translated into OpenFOAM: fvm::laplacian(D1D2,C)+ (D1+D2)*F/R/T* (fvc::grad(C,Ue)+C*fvc::laplacian(Ue)) In my understading it should be correct, however the "fvc::grad(C,Ue)" terms with the following error message: "ElectricEqn.H:8:16: note: ‘Foam::volScalarField {aka Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>}’ is not derived from ‘const Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >’ fvc::grad(C,Ue)" I'm really not sure what does that mean, and how could it be solved. I would be glad for any help. 

August 17, 2015, 06:55 

#4 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
I'm not quite sure where the problem comes from without seeing the type definitions of the variables.
Some ideas: Are you trying to solve for both C and Ue? If so, you will need a second equation for Ue. Also keep in mind that by explicitly using C in the equation the solver will use the field from the last iteration there. A tmp<> is basically a pointer to the object to save performance by avoiding copies at assignments. You can dereference it with (), which is sometimes needed. 

August 19, 2015, 05:00 

#5 
New Member
Gergely Acs
Join Date: Jun 2015
Location: Hong Kong
Posts: 4
Rep Power: 9 
Thank you for the suggestions.
Yes, I have two equations, one is the previously discussed transport equation, the other is a "classic" NernstPlanck equation, where the timederivative of the concentration equals the gradient of the concentration and the potential field. With that I don't have any problems. It is actually very similar to the other equation, the main difference is that the timederivative is missing due to some analytical reasons. Lately I've made a formulation which will enable my code to run, where I've separated my equation in the following way: FluxUe=(D1+D2)*F/R/T*C*fvc::grad(Ue); fvScalarMatrix UeEqn ( fvm::laplacian(D1D2,C)+fvc::div(FluxUe) ); However, this fails to run (Floating point exception) wherever my potential field is not uniform, thus I fear this formulation is essentially incorrect. I'm attaching my code with a very simple 1D testcase, if that would help. In this state may it doesn't makes sense, why would I use OpenFoam, but the next step will be to couple this system with the incompressible NavierStokes equations. 

August 20, 2015, 06:49 

#6 
Senior Member
Join Date: Oct 2013
Posts: 397
Rep Power: 17 
In which of these two lines is it failing?
Have you checked for division by zero in the first line? 

August 20, 2015, 08:39 

#7 
New Member
Gergely Acs
Join Date: Jun 2015
Location: Hong Kong
Posts: 4
Rep Power: 9 
Actually none of these two lines are failing but when I call the solve command with "UeEqn.solve();", only in such a case when the first line creates a nonzero field.
The divisions in the first line are only constant positive scalars, and the whole FluxUe field had only finite values, in other words there wasn't any NaN or Inf values. 

Tags 
custom solver, divergence 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
fvc::div for surfaceVectorField  ARTem  OpenFOAM Programming & Development  4  November 2, 2018 12:41 
fvc::div() instead of fvc::flux()?  Franko  OpenFOAM Programming & Development  1  December 23, 2014 11:11 
fvc::div() strange behaviour  ivan_cozza  OpenFOAM Running, Solving & CFD  2  February 6, 2010 06:09 
help needed with fvc::div()  johanna  OpenFOAM Programming & Development  2  August 31, 2009 07:12 
What type return fvcdiv  su_junwei  OpenFOAM Running, Solving & CFD  6  October 13, 2008 07:09 