Simplifying equation help

mcfdma · February 17, 2021, 13:51

Hello.

I need some guidance with writing equations in my header file.

I am working on electrohydrodynamic (EHD) solver that uses interFoam+electrostatics elements (http://openfoamwiki.net/index.php/Co...ectromagnetics)

As you can see from the above link, the third equation that is the charge conservation equation is written as (the two equations are the same where 1st one is used in the link above and the other one is the same but written using different variables).
Screenshot 2021-02-17 174400.png

In my header file, I have written

Code:

        tmp<fvScalarMatrix> ECDEqn
        (
            fvm::ddt(rhoE) 
          - fvc::laplacian(sgm, Ue)
          + fvm::div(phi, rhoE)
        );

        solve(ECDEqn);

However, I have read an article [1] that simplifies this equation into:Screenshot 2021-02-17 173711.png
where K is equivalent to the previous equation's sigma.

I have tried rewriting my header file equation but only this seems to be working. Is it correct?

Code:

        tmp<fvScalarMatrix> ECDEqn
        (
          fvm::laplacian(sgm, Ue)

        );

        solve(ECDEqn);

I tried many ways including solve(ECDEqn == 0) but was unsuccessful.

Can someone guide me on how to write the simplified equation?

Many thanks

[1] https://iopscience.iop.org/article/1...17/23/1/015004

Tobermory · February 18, 2021, 17:43

Umm ... I am struggling to relate your coding to your equations. Just for clarity:

Code:

        tmp<fvScalarMatrix> ECDEqn
        (
            fvm::ddt(rhoE) 
          - fvc::laplacian(sgm, Ue)
          + fvm::div(phi, rhoE)
        );

        solve(ECDEqn);

translates into the following:
$\frac{\partial \rho_E}{\partial t} - \nabla \cdot (\sigma \nabla U_e) + \nabla \cdot (\rho U \rho_E)$

or perhaps the following if the code takes an incompressible flow approach ( $\rho=1$ ):
$\frac{\partial \rho_E}{\partial t} - \nabla \cdot (\sigma \nabla U_e) + \nabla \cdot (U \rho_E)$

which is probably not what you are after ... or is it? If not, then it's probably a good idea to read up on the meaning of laplacian & div. Good luck my friend.

mcfdma · February 19, 2021, 03:55

Many thanks for your reply. Yes, I am aware of laplacian and divergence.

I believe I made a blunder mistake in writing my second equation (Apologies):

$\frac{\partial \rho_E}{\partial t}+\nabla \cdot \vec{J} =0$

and $\vec{J} = \rho_E \vec{u} + \sigma \vec{E}$

$\frac{\partial \rho_E}{\partial t}+\nabla \cdot (\rho_E \vec{u}) + \nabla \cdot (\sigma \vec{E}) = 0$

Using Gauss Law $\nabla \cdot (\varepsilon \vec{E}) = \rho_E$ along with $\vec{E} = -\nabla Ue$

I get,

$\frac{\partial \rho_E}{\partial t}+\nabla \cdot \rho_E u - \nabla \cdot (\sigma \nabla Ue) = 0$

Code:

 tmp<fvScalarMatrix> ECDEqn
        (
            fvm::ddt(rhoE) 
          + fvm::div(phi, rhoE)
          - fvc::laplacian(sgm, Ue)
        );

        solve(ECDEqn);

which is equivalent to the below equation EHDFoam (equation 3 within the link),

$\frac{\partial \rho_E}{\partial t}+\nabla \cdot \rho_E u = -\frac{\sigma}{\varepsilon}\rho_E$

Have I done anything wrong? I have been getting appropriate results but wanted to implement this simplified part,

$\nabla . (\sigma \vec{E})= 0$

to observe the change but wanted to know if what I had written was correct.

Code:

        tmp<fvScalarMatrix> ECDEqn
        (
          fvm::laplacian(sgm, Ue)

        );

        solve(ECDEqn);

This was the only success but not sure changing it from fvc to fvm is right. Any suggestions? Many thanks in advance

Tobermory · February 19, 2021, 05:12

Aaah - that makes sense now, and yes your coding looks good to me.

In case you are interested, the difference between fvm and fvc is subtle: fvm builds the matrix coefficients for an implicit representation of the term (i.e. expresses the equation in terms of the variable at the new time step), whilst fvc uses an explicit representation (using old time values only) ... so it makes sense that at least one of your terms in the fvMatrix must be an implicit (fvm) term.

mcfdma · February 19, 2021, 05:18

Thanks for the quick fvm and fvc explanation.

Could I double check when you say my coding looks okay, do you mean the original equation with 3 terms or the simplified one or all of them?

Tobermory · February 19, 2021, 05:22

Both look fine to me - did they compile and run okay?

mcfdma · February 19, 2021, 05:29

Thanks for the quick response, much appreciated.

Yes, it does compile alright but since the original equation has now been simplified, I have run a case but there has been no charge density $\rho_E$ at that surface of the fluid. I am left confused which is why I posted thinking I made a mistake in my coding.

Would it be possible for you to have a quick look at my solver?

Tobermory · February 19, 2021, 05:54

By all means, I can take a look at your code, although it has been 30yrs or so since I did my course on electromagnetism!

I would suggest checking your logic first, though. Note that your second (simplified) equation just calculates a divergence-free

field, but you stated earlier that:
$\rho_E = \nabla \cdot (\varepsilon \vec{E}) \approx \varepsilon \nabla \cdot \vec{E}=0$
and so does this not constrain the charge density $\rho_E$ to zero? Perhaps I am missing something.

Note also that I am assuming that variables $\varepsilon$ and $\sigma$ are constants, since you are moving them inside and outside of the integrals.

mcfdma · February 19, 2021, 06:19

Thank you very much

30 years? wow, that's amazing!

This part is correct
$\rho_E = \nabla \cdot (\varepsilon \vec{E}) \approx \varepsilon \nabla \cdot \vec{E}$

but your equation ending with equals to 0 should be for the simplified equation including conductivity $\nabla . (\sigma \vec{E})= 0$
This above equation is the leaky-dielectric model.

I have created a dropbox link that contains my actual solver (not the leaky-dielectric model).

You are absolutely correct, I am assuming these 2 parameters $\varepsilon$ and $\sigma$ constant along with the surface tension at the interface.

Tobermory · February 19, 2021, 08:23

I am probably missing something here, but let me try reiterating my above point. We start with the charge conservation equation:

$\frac{\partial \rho_E}{\partial t}+\nabla \cdot (\rho_E \vec{u}) + \nabla \cdot (\sigma \vec{E}) = 0$

$\frac{D \rho_E}{D t}= -\nabla \cdot (\sigma \vec{E})$

where D/Dt is the material derivative since the fluid is incompressible; and then you want to set:

$\nabla \cdot(\sigma E) = 0$

which means that:

$\frac{D \rho_E}{D t}= 0$

i.e. the charge density field will be unchanging, or at least is just convected through the domain if the inlet boundary boundary value is non-zero.

Indeed, looking briefly at your EHDEqn.H code, the ECDEqn is where the rhoE value is calculated/updated and so if you change this to your new, simplified version of the equation then where is rhoE calculated now? It is probably just remaining at its initial field definition. You would need to add another expression to calculate rhoE from the electric field since the EHD equation no longer solves for rhoE.

Does that explain why you are not seeing any charge density when you use this simplified version?

I don't have access to the paper that you referenced at the start, nor to your model BCs, but I would suggest thinking again about those to make sure that the simplified equation makes sense for your application, and that your BCs are correct. Good luck.

February 18, 2021, 17:43		#2
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 746 Rep Power: 14	Umm ... I am struggling to relate your coding to your equations. Just for clarity: Code: tmp<fvScalarMatrix> ECDEqn ( fvm::ddt(rhoE) - fvc::laplacian(sgm, Ue) + fvm::div(phi, rhoE) ); solve(ECDEqn); translates into the following: $\frac{\partial \rho_E}{\partial t} - \nabla \cdot (\sigma \nabla U_e) + \nabla \cdot (\rho U \rho_E)$ or perhaps the following if the code takes an incompressible flow approach ( $\rho=1$ ): $\frac{\partial \rho_E}{\partial t} - \nabla \cdot (\sigma \nabla U_e) + \nabla \cdot (U \rho_E)$ which is probably not what you are after ... or is it? If not, then it's probably a good idea to read up on the meaning of laplacian & div. Good luck my friend.

February 19, 2021, 03:55		#3
mcfdma Member X Join Date: Jan 2019 Posts: 63 Rep Power: 7	Many thanks for your reply. Yes, I am aware of laplacian and divergence. I believe I made a blunder mistake in writing my second equation (Apologies): $\frac{\partial \rho_E}{\partial t}+\nabla \cdot \vec{J} =0$ and $\vec{J} = \rho_E \vec{u} + \sigma \vec{E}$ $\frac{\partial \rho_E}{\partial t}+\nabla \cdot (\rho_E \vec{u}) + \nabla \cdot (\sigma \vec{E}) = 0$ Using Gauss Law $\nabla \cdot (\varepsilon \vec{E}) = \rho_E$ along with $\vec{E} = -\nabla Ue$ I get, $\frac{\partial \rho_E}{\partial t}+\nabla \cdot \rho_E u - \nabla \cdot (\sigma \nabla Ue) = 0$ Code: tmp<fvScalarMatrix> ECDEqn ( fvm::ddt(rhoE) + fvm::div(phi, rhoE) - fvc::laplacian(sgm, Ue) ); solve(ECDEqn); which is equivalent to the below equation EHDFoam (equation 3 within the link), $\frac{\partial \rho_E}{\partial t}+\nabla \cdot \rho_E u = -\frac{\sigma}{\varepsilon}\rho_E$ Have I done anything wrong? I have been getting appropriate results but wanted to implement this simplified part, $\nabla . (\sigma \vec{E})= 0$ to observe the change but wanted to know if what I had written was correct. Code: tmp<fvScalarMatrix> ECDEqn ( fvm::laplacian(sgm, Ue) ); solve(ECDEqn); This was the only success but not sure changing it from fvc to fvm is right. Any suggestions? Many thanks in advance

February 19, 2021, 05:12		#4
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 746 Rep Power: 14	Aaah - that makes sense now, and yes your coding looks good to me. In case you are interested, the difference between fvm and fvc is subtle: fvm builds the matrix coefficients for an implicit representation of the term (i.e. expresses the equation in terms of the variable at the new time step), whilst fvc uses an explicit representation (using old time values only) ... so it makes sense that at least one of your terms in the fvMatrix must be an implicit (fvm) term. tonnykz likes this.

February 19, 2021, 08:23		#10
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 746 Rep Power: 14	I am probably missing something here, but let me try reiterating my above point. We start with the charge conservation equation: $\frac{\partial \rho_E}{\partial t}+\nabla \cdot (\rho_E \vec{u}) + \nabla \cdot (\sigma \vec{E}) = 0$ $\frac{D \rho_E}{D t}= -\nabla \cdot (\sigma \vec{E})$ where D/Dt is the material derivative since the fluid is incompressible; and then you want to set: $\nabla \cdot(\sigma E) = 0$ which means that: $\frac{D \rho_E}{D t}= 0$ i.e. the charge density field will be unchanging, or at least is just convected through the domain if the inlet boundary boundary value is non-zero. Indeed, looking briefly at your EHDEqn.H code, the ECDEqn is where the rhoE value is calculated/updated and so if you change this to your new, simplified version of the equation then where is rhoE calculated now? It is probably just remaining at its initial field definition. You would need to add another expression to calculate rhoE from the electric field since the EHD equation no longer solves for rhoE. Does that explain why you are not seeing any charge density when you use this simplified version? I don't have access to the paper that you referenced at the start, nor to your model BCs, but I would suggest thinking again about those to make sure that the simplified equation makes sense for your application, and that your BCs are correct. Good luck. tonnykz and mcfdma like this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Solving coupled steady and unsteady equation in OpenFOAM	jd87	OpenFOAM Running, Solving & CFD	2	July 22, 2020 12:36
mass flow in is not equal to mass flow out	saii	CFX	12	March 19, 2018 06:21
Some problem of "Qcriterion.mcr&	yuyuxuan	Tecplot	9	February 12, 2016 04:27
Need help:about energy equation in CFX	Stein	CFX	4	July 2, 2009 23:31
Diffusion Equation	izardy amiruddin	Main CFD Forum	2	July 4, 2002 09:14

February 19, 2021, 05:18		#5
mcfdma Member X Join Date: Jan 2019 Posts: 63 Rep Power: 7	Thanks for the quick fvm and fvc explanation. Could I double check when you say my coding looks okay, do you mean the original equation with 3 terms or the simplified one or all of them?

February 19, 2021, 05:22		#6
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 746 Rep Power: 14	Both look fine to me - did they compile and run okay?

February 19, 2021, 05:29		#7
mcfdma Member X Join Date: Jan 2019 Posts: 63 Rep Power: 7	Thanks for the quick response, much appreciated. Yes, it does compile alright but since the original equation has now been simplified, I have run a case but there has been no charge density $\rho_E$ at that surface of the fluid. I am left confused which is why I posted thinking I made a mistake in my coding. Would it be possible for you to have a quick look at my solver?

February 19, 2021, 05:54		#8
Tobermory Senior Member Join Date: Apr 2020 Location: UK Posts: 746 Rep Power: 14	By all means, I can take a look at your code, although it has been 30yrs or so since I did my course on electromagnetism! I would suggest checking your logic first, though. Note that your second (simplified) equation just calculates a divergence-free field, but you stated earlier that: $\rho_E = \nabla \cdot (\varepsilon \vec{E}) \approx \varepsilon \nabla \cdot \vec{E}=0$ and so does this not constrain the charge density $\rho_E$ to zero? Perhaps I am missing something. Note also that I am assuming that variables $\varepsilon$ and $\sigma$ are constants, since you are moving them inside and outside of the integrals.

February 19, 2021, 06:19		#9
mcfdma Member X Join Date: Jan 2019 Posts: 63 Rep Power: 7	Thank you very much 30 years? wow, that's amazing! This part is correct $\rho_E = \nabla \cdot (\varepsilon \vec{E}) \approx \varepsilon \nabla \cdot \vec{E}$ but your equation ending with equals to 0 should be for the simplified equation including conductivity $\nabla . (\sigma \vec{E})= 0$ This above equation is the leaky-dielectric model. I have created a dropbox link that contains my actual solver (not the leaky-dielectric model). You are absolutely correct, I am assuming these 2 parameters $\varepsilon$ and $\sigma$ constant along with the surface tension at the interface.