How can I modify the equation for electric potential?

[rodion]

Hello everyone,

I work with simulation of electric arc plasma torches and currently I am trying to improve the theory of the electric arc model of Code_Saturne,
For that I developed a two-temperature model (one temperature fore electrons and one temperature for heavy species) of electric arc in Code_Saturne. Now in my two-temperature electric arc model I need to introduce the effective field instead of the common electric field in the equation of electric current continuity. Right now Code_Saturne solves the following equation in terms of electric potential:
div(σ∇φ)=0
Basically this equation is the diffusive term being equal zero. It is given in the section "cs elec model routine" of Code_Saturne theory guide. The further improvement of the model that I am trying to introduce consists in correcting the previous equation to the following one:
div(A(x,y,z)+σ∇φ)=0.
Does anyone have any idea how can I implement this correction? How can I add an extra term inside the divergence of the diffusive term?

Best regards,
Rodion

[Erwan Le Coupanec]

Hello,

To be able to give more adapted advices:
can you give some description of how this vector A is computed ?
what does this term stand for, or what phenomena is it modeling ?

Thanks,
Erwan

[rodion]

Hello Erwan,

Thank you very much for your response.

Basically I am trying to add the effect of electron pressure on the electric field. The current form of the equation for electric potential φ is:
div(σ∇φ)=0 (1)
Or if we write it in terms of electric field div(σE)=0

I am trying to replace the electric field E in this equation by the effective field Ep=E+∇Pe/(e*ne) where Pe is electron pressure and ne is the electron number density. The electron number density is a function of temperature.
The expression for effective field can be rewritten as ∇φp=∇φ+∇(ne*k*Te)/(e*ne) where Te is the electron temperature which I compute in my model.
If we substitute the gradient of electric potential in the equation (1) by the expression for ∇φp, we get:
div(σ(∇φ+∇(ne*k*Te)/(e*ne)))=0
Expanding the gradient of the product we get:
div(σ(∇φ+k(Te*∇ne+ne*∇Te)/(e*ne)))=0 (2)

Since in the model for electric arc the solved thermal variable is enthalpy, we can express ∇ne and ∇Te through the gradient of electron enthalpy he as follows:
∇ne = (dne/dhe) *∇he
∇Te = ∇he /Cpe
Where Cpe is the specific heat of the electron gas and defined as Cpe=dhe/dTe. In general Cpe is defined for constant pressure and the plasma flow is weakly compressible, since the gradient of electron pressure is low enough.
Replacing the gradients ∇ne and ∇Te in the equation (2) by the expressions above we get:
div(σ(∇φ+k(Te*(dne/dhe) *∇he+ne*∇he /Cpe)/(e*ne)))=0
Then we need to reorganize this equation as follows:
div(σ(∇φ+(k/e)*((Te/ne)*(dne/dhe) + 1/Cpe)*∇he))=0

The expression (Te/ne)*(dne/dhe) + 1/Cpe) is basically a function of temperature (or enthalpy, does not matter in this case). Let's replace it as Q(he):
div(σ(∇φ+(k/e)*Q(he)*∇he))=0
The left hand side of this equation can be split into two part:
div(σ∇φ)+div((σ*k/e)*Q(he)*∇he)=0
The first term of the resulting equation correspond to the standard equation for electric potential in Code_Saturne: just the diffusive term. The second term is new, it is a divergence of the field that could be computed separately, since the product (σ*k/e)*Q(he) could be precomputed in a lookup table and the gradient ∇he is already computed in my model in Code_Saturne.

For more details and perhaps better explanation I would suggest also to take a look at the article "Non-Equilibrium Modeling of Arc Plasma Torches" by J. P. Trelles, E. Pfender, J. V.R. Heberlein. It could be found by the following link:
https://arxiv.org/ftp/arxiv/papers/1301/1301.0651.pdf
The page 9 contain a good explanation.

I found files the driflu.f90 and EXAMPLES/cs_user_physical_properties-scalar_drift.f90 that deal with some drift, but I am not really sure that I understand what is going on in those files. Perhaps drift is what I need, but I just don't understand how to use it for my purposes.

Please let me know if my explanation is unclear or you need to know some more details.

Best regards,
Rodion

[Erwan Le Coupanec]

Hello Rodion,

Thank you for the description.

Looking at the model, my first remark is that unlike the model described in this old thread viewtopic.php?f=4&t=1917 that you mentioned, the added term does not write as:
div( phi C u_r ) with u_r a vector, and C a generic factor.

In that former thread, using the drift flux generic model allowed to make such a term implicit in the discretized model. In your case, you first would have to exhibit a dependency on phi in the factor (σ*k/e)*Q(he) and if possible make it a linear dependency.

Have you already tried to take the additional term into account explicitly ? I.e. just adding it to the right hand side, by implementing it in cs_user_source_terms (there is a C version of that user source file now in v6.0, not in v5.0 though).
You could make a call to cs_convection_diffusion_scalar, without forgetting to cancel the convection part, and giving as input the he field values and the diffusion coefficient (σ*k/e)*Q(he).

Of course, having an additional explicit source term can bring additional constraint on the time step. You could perform some sensitivity tests to the time step value and see how the solver behaves.

Regards,
Erwan

[rodion]

Hello Erwan,

Thank you very much for your response.

The thing is that the factor (σ*k/e)*Q(he) is known for each cell and it depends on the enthalpy (or temeprature, does not matter). But this factor does not depend on the electric potential φ, since it is a parameter of plasma and not of the electric field.

Also I already computed the field ∇he. Therefore the term (σ*k/e)*Q(he)*∇he is easy to compute for each cell of the model at each time step. It would be great to creat an additional source term, but the main problem is to compute the divergence of this new custom field. I can compute the divergence of (σ*k/e)*Q(he)*∇he in, let's say, ParaView once or twice, but I cannot do it for every time step. Using the approximate values from ParaView in some cells could be a first step in the development for a steady model but for an unsteady model it will not work.

As far as I can see, the subroutine cs_convection_diffusion_scalar computes the diffusive term for some solved variable, meanwhile I need to compute the divergence of a custom user defined field.

Also if I recall correctly, Code_Saturne can calculate gradients only for solved variables. Please correct if I am wrong. Is it the same for cs_convection_diffusion_scalar?

If so, does it mean that I have to make a copy of the enthalpy he as another solved variable, allot it the diffusivity equal to (σ*k/e)*Q(he) and only like this I could compute the term div((σ*k/e)*Q(he)*∇he) for the source term in the equation for electric potential?

Best regards,
Rodion

[Yvan Fournier]

Hello Rodion,

Regarding gradients, code_saturne can also (at least startting with version 6.0) compute gradients of non-resolved fields such as properties, with some limitations:

- for variables, gradients take boundary conditions into account
- for other fields, such as properties, homogeneous Neuman BC's are assumed everywhere; if you need more specific boundary conditions, you can vompute them, but instead of using the "high level" cs_field_gradient functions, you need to use the lower level cs_gradient functions and provide the coefa/coefb arrays for BC's.

For the rest, I'll let Erwan complete this since he is much more familiar with the relevant numerical aspects.

Best regards,

Yvan

[rodion]

Hello Yvan,

Thanks a lot for your comment.
Unfortunately I am still using Code_Saturne 5.0.9 on Windows. The transition to the version 6.0 or 6.1 and to Linux will take some more time, since I have to wait for the end of the lock-down/confinement in France to create a WSL machine on the multiprocessor computer.
Therefore for now I need to somehow compute the divergence of a custom user defined field in Code_Saturne 5.0.9 and add it as a source term.

Best regards,
Rodion

[Yvan Fournier]

Hello Rodion,

In version 5.0, you can still use the low-level gradient functions with any array.

Depending on how long the lock-down lasts though, I can't predict whether spending more time on low-level stuff with v5.0 or waiting to install 6.0 is the best option... unless you boundary conditions are not all "homogeneous Neumann", in which case for non-resolved variable fields, you'll need to use the lower-level functions anyways...

Best regards,

Yvan

[rodion]

Hello Yvan,

Thanks a lot for your answer.
The low-level gradient function that you are talking about, is it cs_field_gradient_scalar or cs_gradient_scalar or cs_gradient_scalar?

I think in my case the homogeneous Neumann boundary condition should be alright. It is the zero flux boundary condition, right? Anyway, I have no idea what values should enthalpy gradient have at the boundaries. Thus I may assume only zero flux.

Best regards,
Rodion

[Yvan Fournier]

Hello Rodion,

cs_field_gradient_scalar is the high level functions, while cs_gradient_scalar is the low-level one.

Best regards,

Yvan