CFD Online Discussion Forums - Simplifying equation help

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Simplifying equation help (https://www.cfd-online.com/Forums/openfoam-programming-development/233938-simplifying-equation-help.html)

mcfdma

February 17, 2021 12:51

Simplifying equation help

3 Attachment(s)

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).
Attachment 82797

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:Attachment 82796
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 16: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 02: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 04: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 04: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 04:22

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

mcfdma

February 19, 2021 04: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 04: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 05: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 07: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.

All times are GMT -4. The time now is 03:14.