Dielectric Charge Buildup Boundary Condition

[Shane Keniley]

I've been putting this off for a long time, but I think adding a surface charge boundary condition is something I can no longer avoid. The Economou dielectric BC Corey added for his APS test case is good, but as far as I can tell it is only valid for thin dielectric layers. A problem I'm comparing to has a DBD with a 1 mm discharge region and a 1 mm thick dielectric on both electrodes, so I don't think that BC is valid for this situation.  If possible, it would be good to have the same kind of surface charge BC that other plasma fluid solvers have. I've attached the relevant equations here. Basically, I want a BC/InterfaceKernel that accounts for the charge buildup on a dielectric surface. This takes the form of a distributed ODE on the surface, with Je + Ji being the total charged particle flux to that boundary. In the input file I think this would take the form of a BC on the plasma-side and an InterfaceKernel on the dielectric side, so we would have a multi-block problem. 

I have no idea how to implement this beyond that though. As far as I can see, here are the options we have: 

• Stateful material property. This would have the ODE solved explicitly using a simple forward difference scheme. Alex added the capability to access normals to the SigmaMat material last time I tried working on this so the framework for this one is in place. I previously tried this approach and it failed miserably, though I don't know if that was because the approach was flawed or because my implementation was wrong. (Honestly, the latter is more likely. This is worth a second look.)
• Nodal kernels. I know these are built for his kind of thing, but unfortunately I don't think this works for us. This BC requires gradients of nonlinear variables (potential) and material properties, some of which are nonlinear (electron transport coefficients). I don't think either is accessible through a nodal kernel. Maybe the total flux could be computed through a postprocessor and passed to the NodalKernel? That would require running with PJFNK since the jacobian contributions wouldn't be included, but that's not dealbreaking. 
• IntegratedBC. I've seen in this MOOSE user group thread that someone suggested just lumping this kind of BC into an integrated BC, but I don't really understand how to do that. How can you lump an ODE into an IntegratedBC? 

That's all I've got. If anyone else has any insight into this, I would be happy to hear your suggestions. 

[Shane Keniley]

All,

I'm still really struggling with this problem. I tried adding a SideCurrent postprocessor to calculate the total flux to each wall, and then I tried including a simple modification to the InterfaceKernel example shown here to add the surface charge value. Whatever I do simply doesn't work; the problem usually just runs into floating point errors (which I'm still having in opt mode, somehow). The problem is that I cannot tell if my InterfaceKernel is wrong, if there is some problem with the SideCurrent postprocessor, if I am applying my BCs/InterfaceKernels incorrectly in the input file, or if my problem is simply ill-posed. I don't know of any simple model that I can compare to here. Here is my procedure: 

1. Created 1D gmsh file with three regions: two dielectrics with the discharge gap in the middle, grounded wall on the right and driven electrode on the left. 
2. Added postprocessor to compute right hand side of surface charge ODE on each side. 
3. Applied an InterfaceKernel (InterfaceDiffusion, modified to include surface charge) on the plasma side of each dielectric to enforce surface charge BC (see figure attached to previous message)
4. Applied a MatchedValueBC on the dielectric-side of each interface to enforce potential continuity

The InterfaceDiffusion kernel looks like this: 

switch (type)
{
    case Moose::Element:
        r = -_test[_i][_qp] * (_D_neighbor[_qp] * _grad_neighbor_value[_qp] * _normals[_qp]);
    case Moose:Neighbor:
        r = -_test_neighbor[_i][_qp] * (_D[_qp] * _grad_u[_qp] * _normals[_qp] + _sigma[_qp]);
}

Note that I have no idea if this is correct. I've tried every possible combination of signs on each term, along with every possible combination of adding/subtracting _sigma to both the Element side and Neighbor side. The problem either just explodes (e.g. the potential will skyrocket/plummet at one or both of the interfaces, or it runs but eventually runs into floating point errors and breaks). I don't understand how to derive the flux boundary condition when surface charge is involved.

For _sigma, I created a ScalarVariable with an ODETimeDerivative kernel and a second kernel that just takes the value of the SideCurrent postprocessor. That looks like this: 

[Variables]
    [./sigma_left]
        order = FIRST
        family = SCALAR
        initial_condition = 0
    [../]
    [./sigma_right]
        order = FIRST
        family = SCALAR
        initial_condition = 0
    [../]

    ...(other variables)
[]

[ScalarKernels]
    [./sigma_left_dt]
        type = ODETimeDerivative
        variable = sigma_left
    [../]
    [./sigma_left_rhs]
        type = BoundaryFlux
        current = left_current
        variable = sigma_left
    [../]

    [./sigma_right_dt]
        type = ODETimeDerivative
        variable = sigma_right
    [../]
    [./sigma_right_rhs]
        type = BoundaryFlux
        current = right_current
        variable = sigma_right
    [../]
[]


[PostProcessors]
    [./left_current]
        type = SideCurrent
        variable = em
        mobility = muem
        potential = potential
        r = 0
        mean_en = mean_en
        ions = 'Arp Ar2p'
        Arp = Arp  # deprecated; doesn't do anything but hasn't been removed yet
        position_units = ${dom0Scale}
        boundary = 'master10_interface'
    [../]
    [./right_current]
        type = SideCurrent
        variable = em
        mobility = muem
        potential = potential
        r = 0
        mean_en = mean_en
        ions = 'Arp Ar2p'
        Arp = Arp  # deprecated; doesn't do anything but hasn't been removed yet
        position_units = ${dom0Scale}
        boundary = 'master12_interface'
    [../]
[]

where BoundaryFlux is a scalarkernel I added that just applies the value from the postprocessors as the value of the scalar residual. This approach won't work in 2D-3D, but I think it makes sense in 1D since our interface is a single node and as such the current has a single (scalar) value.  And obviously this is all run with PJFNK since I can't add the contributions of the postprocessor to the jacobian. 

An example of this was added to my surface_charge branch. The relevant input file is argon_csv_interface.i. As it is now, the file simply runs into an MPI_Abort error at timestep 215. (Debug mode tells me nothing; it just runs into a floating point error earlier in the simulation.)

Sorry, I know this is a long message but any guidance here would be much appreciated. I've been struggling with this problem on and off for about two years, and I'm just completely stuck at this point. I can't figure out how to make this work. 

[Alexander Lindsay]

If you've been struggling with it for a couple years, then it's probably a pretty challenging problem which will take a nontrivial amount of time to work through. I'm curious whether we could set up some contract work for an INL employee (me) to look at it because I don't think I can justify working on this task with one of my existing charge numbers. I haven't looked into such a thing. I'll talk to  my other MOOSE team members. I'm curious if Davide or Steve has any thought on it.

--
You received this message because you are subscribed to the Google Groups "zapdos-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/06bcde48-2bf4-4c5b-a2ac-e36877356a81%40googlegroups.com.

[Shane Keniley]

If you've been struggling with it for a couple years, then it's probably a pretty challenging problem which will take a nontrivial amount of time to work through.

It's certainly challenging for me! I haven't been working on it consistently over 2 years, just on and off whenever I get a chance. 

I'm curious whether we could set up some contract work for an INL employee (me) to look at it because I don't think I can justify working on this task with one of my existing charge numbers. I haven't looked into such a thing. I'll talk to  my other MOOSE team members. I'm curious if Davide or Steve has any thought on it.

That would certainly be helpful. I have a meeting with Davide tomorrow afternoon so I'll mention it to him. 

In the meantime I'll try to simplify this problem a bit...maybe start with a constant surface charge so I can make sure the interfacekernel is correct and remove all the postprocessor nonsense. 

On a side note, I tried using MeshSideSetGenerator to create a sideset block that the surface charge would live on. It either doesn't work on 1D or I was using it incorrectly (I get the error "Cannot reinit 0D LAGRANGE elements!"), but do you think that would be a viable path forward in 2D? I'm not sure how InterfaceKernels would work in that case, but I was going to try to set up a surface charge variable on that sideset. 

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos...@googlegroups.com.

[Alexander Lindsay]

On Mon, Apr 13, 2020 at 3:14 PM Shane Keniley <keni...@illinois.edu> wrote:

If you've been struggling with it for a couple years, then it's probably a pretty challenging problem which will take a nontrivial amount of time to work through.

It's certainly challenging for me! I haven't been working on it consistently over 2 years, just on and off whenever I get a chance. 

I'm curious whether we could set up some contract work for an INL employee (me) to look at it because I don't think I can justify working on this task with one of my existing charge numbers. I haven't looked into such a thing. I'll talk to  my other MOOSE team members. I'm curious if Davide or Steve has any thought on it.

That would certainly be helpful. I have a meeting with Davide tomorrow afternoon so I'll mention it to him. 

In the meantime I'll try to simplify this problem a bit...maybe start with a constant surface charge so I can make sure the interfacekernel is correct and remove all the postprocessor nonsense. 

On a side note, I tried using MeshSideSetGenerator to create a sideset block that the surface charge would live on. It either doesn't work on 1D or I was using it incorrectly (I get the error "Cannot reinit 0D LAGRANGE elements!"), but do you think that would be a viable path forward in 2D? I'm not sure how InterfaceKernels would work in that case, but I was going to try to set up a surface charge variable on that sideset. 

Variables cannot be boundary restricted; they must be block restricted. It sounds like you want to create a lower dimensional block. We have a LowerDBlockFromSidesetGenerator that can do that.

The problem then is that we don't currently have a means for pulling that lower dimensional solution into an interface kernel or into a boundary condition. Which of these would you want? Both? Let me know and I will create a MOOSE issue. This would be another good thing for me to contract on...

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/21b2c44a-f809-4008-9571-fe6a04876163%40googlegroups.com.

[Shane Keniley]

Variables cannot be boundary restricted; they must be block restricted. It sounds like you want to create a lower dimensional block. We have a LowerDBlockFromSidesetGenerator that can do that.

Yeah, that one sounds more like what I wanted. I was under the impression that MeshSideSetGenerator did that, but I guess I misinterpreted it! 

The problem then is that we don't currently have a means for pulling that lower dimensional solution into an interface kernel or into a boundary condition. Which of these would you want? Both? Let me know and I will create a MOOSE issue. This would be another good thing for me to contract on...

I have no idea if this is possible, but this is what I was trying to set up: 

1. Create lower-dimensional block on dielectric surface, with a variable corresponding to surface charge, "sigma" (e.g. in a 2D domain, sigma would be a 1D variable living on the interface between two blocks)
2.  Accept as input the fluxes from charged particles on the plasma side as sigma's source term (directly calculating these with variable gradients and material properties would be ideal, but I have no idea if that's possible...if not, a postprocessor or userobject of some sort should do the trick with PJFNK)
3. Solve sigma ODE implicitly each timestep
4. Access the variable sigma in the InterfaceKernel to enforce dielectric surface charge boundary condition

So I would like to be able to (a) accept variable, gradient, and material values on lower dimensional mesh from an adjacent block, and (b) access lower dimensional block's variable value(s) on the boundary's InterfaceKernel. (Again, this is assuming either one is possible. I'm not sure about (a); I know material properties don't make sense for nodal values, and I'm assuming that's what a lower dimensional mesh becomes?)

[Alexander Lindsay]

In terms of MOOSE objects what you want boils down to a conforming mesh mortar object. In 2D you can do this right now. Maybe look at the tests in moose/test/tests/mortar/continuity-2d-conforming for guidance, and then you can ask more questions.

Now this will  be a lot slower than it could be because this will be doing geometric searches to determine neighbors and things when in fact we will always know the neighbor. So I've created a ticket here: https://github.com/idaholab/moose/issues/15061 that would add new capability. This new thing would also work right out of the box in 1D.

--
You received this message because you are subscribed to the Google Groups "zapdos-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/deb6b173-8e7c-4f4a-92f2-7a4aaf4c6615%40googlegroups.com.

[Steven Shannon]

It could be an interesting problem in the fusion/mhd effort at INL. We could always put in for a small continuation on the NSF grant if we can justify that this would be a critical problem identified during the initial two year effort that if solved would greatly increase the user base and utility of Zapdos...

--
You received this message because you are subscribed to the Google Groups "zapdos-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/deb6b173-8e7c-4f4a-92f2-7a4aaf4c6615%40googlegroups.com.

[Shane Keniley]

Okay, I have an exciting update! (Exciting to me, anyway.) My approach that I detailed in my post here was sound. Using the SideCurrent postprocessor as the right-hand side of the surface charge ODE and applying the _sigma values as scalar variables works, but there were several sign convention errors that stacked up and caused all of my previous attempts to fail. I think I finally have it right now! I've attached plots of my results here. Note that this is a 1D simulation; the x axis is the gap distance and the y axis is time in these plots. Compare those to figures 2, 3, and 6 in Comsol's DBD example here (there is a pdf available on the "download the application files" tab): 

https://www.comsol.com/model/dielectric-barrier-discharge-8637

It looks very comparable. I think this is a good verification of Zapdos!

While this works, I think the ticket Alex opened will still be very helpful. In principle surface charge is a spatial quantity, so in 2D (and god forbid, 3D) simulations it will need to be some kind of lower-dimensional variable that accepts the spatially varying charged particle fluxes to the surface as a source term. I believe doing it that way would also allow me to implement the proper jacobians which would help convergence -- this simulation behaved alright PJFNK, but I wouldn't be surprised if it starts struggling when there are multiple ionic species and the current calculation becomes increasingly nonlinear. 

[Alexander Lindsay]

Very, very cool Shane! Thank you for sharing! Are you gonna get a PR in with some tests? Love what you guys are doing!

--

You received this message because you are subscribed to the Google Groups "zapdos-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/bcb44a11-379a-4153-af29-4e1357789bb6%40googlegroups.com.

[Shane Keniley]

Yeah, I'll roll this into a PR with its own test by the end of the week. This didn't actually require much - just a scalarkernel, a slight modification to a postprocessor, and a modified InterfaceDiffusion kernel. I just want to clean everything up, make sure it's absolutely correct, and abstract it a bit to allow for an arbitrary number of ionic species before I submit. 

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos...@googlegroups.com.

[Shane Keniley]

Now that I have the units problem figured out, I went back and tried calculating the ODE explicitly via forward euler in a boundary-restricted material again. It works fine! It also works with an AuxKernel, but using a Material also allows me to compute the jacobian terms.

This surface charge calculation only works if I use the same identical charged particle fluxes as I use in their BCs. For example, I now have a material called "HagelaarSurfaceCharge", wherein I essentially just copy-and-pasted the ion and electron flux calculations from HagelaarElectronBC, HagelaarIonAdvectionBC, HagelaarIonDiffusionBC, and HagelaarSecondaryElectronBC. This works, but it means that I need a different material (or several if statements) to handle the three possible charged particle BCs included in Zapdos right now.

So I guess the question is, should I create multiple different materials, one for each boundary condition type (e.g. SakiyamaSurfaceCharge, HagelaarSurfaceCharge, LymberopoulosSurfaceCharge) so that surface charge can be computed with any BC? Or is there maybe a cleaner way of doing this? I don't like having so much duplicate code between the BCs and this material, but I'm not sure if there's an alternative. I don't think it's necessarily bad to have three different materials for this either. I'm just looking for the most easy-to-use solution before I move forward with a PR. 

[Corey DeChant]

Good morning,

I made those BCs to better match the BCs used in the papers I was comparing Zapdos with, but ideally they are just limiting version of the Hagelaar BCs. If you want to include the Sakiyama and Lymberopoulos BCs into the new surface charge calculation, I think the most user friendly option (though not the cleanest for the .C file) is to include all of them into one Object with user options based on the assumptions of the BCs (the default option being the Hagelaar BCs). I have a descriptive section of the difference/assumptions for each BC for a paper I am working, that I can send you, or a can post a short list on this forum.

Thank you,

Corey DeChant

[Shane Keniley]

Corey,

Thanks for the feedback! No need to send that info; I'm familiar with the papers. I know that Hagelaar's BC is the more complete version, but since those other BCs exist I want this to be able to fully support them even if they're not used often. I've also been thinking of adding the BCs COMSOL uses (they are similar to the Lymberopoulos BCs, but include an ion sticking coefficient formulation as well). 

 

 If you want to include the Sakiyama and Lymberopoulos BCs into the new surface charge calculation, I think the most user friendly option (though not the cleanest for the .C file) is to include all of them into one Object with user options based on the assumptions of the BCs (the default option being the Hagelaar BCs).

Yeah, that's probably the most straightforward way of doing this. I agree that it's not going to make for a nice, clean .C file, but it probably is easier from an end user's point of view to have one "SurfaceCharge" material rather than multiple materials to choose from. Unfortunately I don't think there's any perfect solution here.  If I could somehow have a material inherit from a BC and access that BC's residual function that would be a bit nicer from a coding standpoint (no duplicate code!), but I don't think materials work like that. Oh well. 

I'll go with a single material and put together a PR soon! It's always possible to change it later if necessary. 

[Alexander Lindsay]

I think you can keep the code relatively clean. Just make no one function is very long, e.g.:

switch (bc_type)

{

  case Hagelaar:

    computeHagelaarStuffMethod();

    break;

  case Sakiyama:

    computeSakiyamaStuffMethod();

    break;

blah blah blah

}

--
You received this message because you are subscribed to the Google Groups "zapdos-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/409e160f-be8d-45db-ae79-da9e825e6647%40googlegroups.com.

[Shane Keniley]

Tangentially related topic: how difficult would it be to add ADInterfaceKernels? That's the only thing stopping me from being able to use a Newton solver since the jacobian for surface charge is quite complicated. I can't seem to get it 100% correct. I wouldn't mind tackling ADInterfaceKernels myself, but I've never handled anything in the base code myself so I'm not sure where to start. I assume the formulation is similar to ADDGKernels?

On that note, is AD compatible with stateful material properties...?

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos...@googlegroups.com.

[Alexander Lindsay]

Yea quite similar. You can open an issue for MOOSE, and either assign yourself and/or me :-)

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/00a09961-e2b2-4086-aa6b-8427f993c4ee%40googlegroups.com.

[Alexander Lindsay]

ADInterfaceKernel should be a template like ADKernel so that it can handle vector variables as well as standard variables.

[Casey Icenhour]

Hi Shane,

This is something I actually put on my todo list last week, since I'd like to have it for some sintering modeling work I'm doing. I should have time this week to work on it if you don't think you'll be able to, so feel free to assign me if that ends up being the case.

Casey

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/CANFcJrEyD6wn%3Dc4Ge%2BhGGrmN0bs3fpEEyiBLtBkyXSeOk3gtHw%40mail.gmail.com.

[Shane Keniley]

Casey,

It would honestly be great if you could take a look at it! It's something I want to do, but it's not a simulation-breaking problem so I won't be able to justify doing it for a few weeks or so; I'm busy writing at the moment. I'll open an issue now and assign both of us to it just in case though. 

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/00a09961-e2b2-4086-aa6b-8427f993c4ee%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "zapdos-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/CANFcJrEyD6wn%3Dc4Ge%2BhGGrmN0bs3fpEEyiBLtBkyXSeOk3gtHw%40mail.gmail.com.

[Shane Keniley]

Well, the issue is up but I'm not actually sure how to assign people...

https://github.com/idaholab/moose/issues/15307

[Casey Icenhour]

Similar situation here -- it's not breaking my simulations, but soon we'll be moving to temperature and pressure dependent material definitions, so coding up those Jacobians are going to be painful. Would rather get it out of the way now. :-)

Good luck on the writing! Is this the dissertation, or are you putting some journal articles together?

I have assigned both you and me to the issue. Thanks for getting that in!

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/9b027fd7-711d-46fa-af3c-5f7ea06aa6e8%40googlegroups.com.

[Shane Keniley]

Mostly just a chapter of my dissertation, which I will turn into a journal article (or two). I'm working on air plasmas with Zapdos, which means I have a couple dozen species (at least) and somewhere around 100-700 reactions, depending on which species you want to look at. It's really cool that Zapdos and Crane work with something that complicated! 

To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/9b027fd7-711d-46fa-af3c-5f7ea06aa6e8%40googlegroups.com.

[Casey Icenhour]

That's awesome! I look forward to reading it! :-)

To unsubscribe from this group and stop receiving emails from it, send an email to zapdos-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/zapdos-users/e8ff505f-f3ac-4755-9ff5-326d989c9c9c%40googlegroups.com.