Re: Drift Diffusion Equation Setup
Hi Jon, I've been continuing to work on my problem. I've implemented the divergence of the current density mathematically into the positive charge and negative charge density equations. Now, I'm encountering issues with the boundary conditions.. the charges should be flowing towards the center of the system, and ultimately recombine, but the issues I'm seeing is that the charges are flowing out of the region. I would like to keep all of the charges confined to the system and prevent them from flowing out of the left and right boundaries. I've tried implementing neumann boundaries using Pion.faceGrad.constrain(0., where=mesh.exteriorFaces) and the negative charge version, but I'm still seeing the charges flow mostly out of the system. I would appreciate any help! Thank you, Justin import numpy as np import matplotlib.pyplot as plt from scipy import special from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, ImplicitSourceTerm, ConvectionTerm from fipy.tools import numerix ## #''' SET-UP PARAMETERS ''' ## a1,b1,c1,a2,b2,c2 = [1.07114255e+00, 6.50014631e+05, -4.51527221e+00, 1.04633414e+00, 1.99708312e+05, -1.52479293e+00] # Parameters for sum of sines fit (Potential fit) #a = -3930.03590805 #b, c = 3049.38274411, -4.01434474 # Parameters for exponential fit (Charge Density) Not used yet q = 1.602e-19#Elementary Charge pini = 154.1581560721245/q #m^-3 nini = -134.95618729/q #m^-3 k1 = 1.8 p1 = 17 k2 = 17 p2 = 1.8 # Parameters for charge density fit (Susi's fit) l = 0.134901960784314 #Length of system in m nx = 134 #Number of cells in system dx = l/nx #Length of each cell in m x = np.linspace(0,l,nx) #Array to calculate initial values in functions below epsilon_r = 25#Relative permittivity of system epsilon = epsilon_r*8.854e-12 #Permittivity of system C/V*m kb = 1.38e-23 #J/K T = 298 #K f = kb*T/q#Volts mu_n = 1.1e-09/1 #m^2/V*s mu_p = 1.1e-09/1 #m^2/V*s Dn = f * mu_n #m^2/s Dp = f * mu_p #m^2/s k_rec = q*(mu_n+mu_p)/(2*epsilon)*10 #k_rec = 0 #pini*np.exp(a*x) #nini*np.exp(b*x+c) #Equations for exponential charge density fits (Not Used Yet) ## ##''' INITIAL CONDITION FUNCTIONS '''# ## def y01(x): """Initial positive ion charge density""" return pini*((special.gamma(k1+p1))/(special.gamma(k1)*special.gamma(p1))*((x/l)**(k1-1))*(1-(x/l))**(p1-1))/7.3572 def y02(x): Initial negative ion charge density""" return nini*((special.gamma(k2+p2))/(special.gamma(k2)*special.gamma(p2))*((x/l)**(k2-1))*(1-(x/l))**(p2-1))/7.3572 def y03(x): """Initial potential""" return a1*np.sin(b1*x+c1) + a2*np.sin(b2*x+c2) mesh = Grid1D(dx=dx, nx=nx) #Establish mesh in how many dimensions necessary ## #''' SETUP CELLVARIABLES AND EQUATIONS ''' ## #CellVariable - defines the variables that you want to solve for: '''Initial value can be established when defining the variable, or later using 'var.value =' Value defaults to zero if not defined''' Pion = CellVariable(mesh=mesh, name='Positive ion Charge Density', value=y01(x)) Nion = CellVariable(mesh=mesh, name='Negative ion Charge Density', value=y02(x)) potential = CellVariable(mesh=mesh, name='Potential', value=y03(x)) #EQUATION SETUP BASIC DESCRIPTION '''Equations to solve for each varible must be defined: -TransientTerm = dvar/dt -ConvectionTerm = dvar/dx -DiffusionTerm = d^2var/dx^2 -Source terms can be described as they would appear mathematically Notes: coeff = terms that are multiplied by the Term.. must be rank-1 FaceVariable for ConvectionTerm "var" must be defined for each Term if they are not all the variable being solved for, otherwise will see "fipy.terms.ExplicitVariableError: Terms with explicit Variables cannot mix with Terms with implicit Variables." ''' #In English: dPion/dt = -1/q * divergence.Jp(x,t) - k_rec * Nion(x,t) * Pion(x,t) where # Jp = q * mu_p * E(x,t) * Pion(x,t) - q * Dp * grad.Pion(x,t) and E(x,t) = -grad.potential(x,t) # Continuity Equation Pion.equation = TransientTerm(coeff=1, var=Pion) == mu_p * (ConvectionTerm(coeff=potential.faceGrad,var=Pion) + Pion * potential.faceGrad.divergence) + DiffusionTerm(coeff=Dp,var=Pion) - k_rec*Pion*Nion #In English:
Re: Drift Diffusion Equation Setup
When taking the divergence of current density: divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) would the divergence not also be applied to the Pion term, since it is also a function of x. I essentially applied the product rule to obtain the `Pion * potential.faceGrad.divergence` part, which in hindsight should have just been a DiffusionTerm. This is how I would derive the divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) as = DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p * potential.faceGrad, var = Pion) Which would just be the drift term for the current density. So, the entire divergence.(Jp) would be: DiffusionTerm(coeff=-mu_p * Pion, var=potential) + ConvectionTerm(-mu_p * potential.faceGrad, var = Pion) *- DiffusionTerm(coeff=Dp, var=Pion) * to include the diffusion part. Please correct me if this is an incorrect understanding, and thank you so much! Justin Pothoof The University of Washington - Department of Chemistry Pre-Candidacy PhD Student Ginger Group On Fri, Aug 9, 2019 at 1:36 PM Guyer, Jonathan E. Dr. (Fed) via fipy < fipy@nist.gov> wrote: > Justin - > > A couple of things: > - Charge Density is not Pion + Nion, it's Pion - Nion > - What are the terms `Pion * potential.faceGrad.divergence` in > Pion.equation and `Nion * potential.faceGrad.divergence` in Nion.equation? > I don't believe they should be there. > - Your equations are really not coupled. Pion.equation is an implicit > function of only Pion, Nion.equation is an implicit function of only Nion, > and potential.equation is an implicit function of only potential. Binding > these equations together with `&` does not gain you anything. Coupling > comes from having implicit parts more than one variable in your equations, > e.g., k_rec * Nion(x,t) * Pion(x,t) could be > ImplicitSourceTerm(coeff=k_rec*Nion) in Pion.equation and > ImplicitSourceTerm(coeff=k_rec*Pion) in Nion.equation. Likewise, > divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) can either be > ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) or it can be > DiffusionTerm(coeff=-mu_p * Pion, var=potential). > - Boundary conditions in FiPy are no-flux by default, so there's no need > to constrain Pion and Nion. > > - Jon > > > From: Justin Pothoof > Sent: Friday, August 9, 2019 2:51 PM > To: Guyer, Jonathan E. Dr. (Fed); FIPY > Subject: Re: Drift Diffusion Equation Setup > > Hi Jon, > > I've been continuing to work on my problem. I've implemented the > divergence of the current density mathematically into the positive charge > and negative charge density equations. Now, I'm encountering issues with > the boundary conditions.. the charges should be flowing towards the center > of the system, and ultimately recombine, but the issues I'm seeing is that > the charges are flowing out of the region. I would like to keep all of the > charges confined to the system and prevent them from flowing out of the > left and right boundaries. I've tried implementing neumann boundaries using > Pion.faceGrad.constrain(0., where=mesh.exteriorFaces) and the negative > charge version, but I'm still seeing the charges flow mostly out of the > system. I would appreciate any help! > > Thank you, > Justin > > > import numpy as np > import matplotlib.pyplot as plt > from scipy import special > from fipy import Variable, FaceVariable, CellVariable, Grid1D, > ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, > ImplicitSourceTerm, ConvectionTerm > from fipy.tools import numerix > > ## > #''' SET-UP PARAMETERS ''' > ## > > a1,b1,c1,a2,b2,c2 = [1.07114255e+00, 6.50014631e+05, -4.51527221e+00, > 1.04633414e+00, > 1.99708312e+05, -1.52479293e+00] > # Parameters for sum of sines fit (Potential fit) > > #a = -3930.03590805 > #b, c = 3049.38274411, -4.01434474 > # Parameters for exponential fit (Charge Density) Not used yet > > q = 1.602e-19#Elementary Charge > > pini = 154.1581560721245/q #m^-3 > > nini = -134.95618729/q #m^-3 > > > k1 = 1.8 > > p1 = 17 > > k2 = 17 > > p2 = 1.8 > # Parameters for charge density fit (Susi's fit) > > l = 0.134901960784314 #Length of system in m > > nx = 134 #Number of cells in system > > dx = l/nx #Length of each cell in m > > x = np.linspace(0,l,nx) #Array to calculate initial values in functions > below > > > epsilon_r = 25#Relative permittivity of system > > epsilon = epsilon_r*8.854e-12 #Permittivity of system C/V*m > > kb = 1.38e-23 #J/K > > T = 298 #K > > f = kb*T/q#Volts > > mu_n = 1.1e-09/1 #m^2/V*s > > mu_p = 1.1e-09/1 #m^2/V*s > > Dn = f * mu_n #m^2/s > > Dp = f * mu_p #m^2/s > > k_rec = q*(mu_n+mu_p)/(2*epsilon)*10 >
Re: Drift Diffusion Equation Setup
Justin - A couple of things: - Charge Density is not Pion + Nion, it's Pion - Nion - What are the terms `Pion * potential.faceGrad.divergence` in Pion.equation and `Nion * potential.faceGrad.divergence` in Nion.equation? I don't believe they should be there. - Your equations are really not coupled. Pion.equation is an implicit function of only Pion, Nion.equation is an implicit function of only Nion, and potential.equation is an implicit function of only potential. Binding these equations together with `&` does not gain you anything. Coupling comes from having implicit parts more than one variable in your equations, e.g., k_rec * Nion(x,t) * Pion(x,t) could be ImplicitSourceTerm(coeff=k_rec*Nion) in Pion.equation and ImplicitSourceTerm(coeff=k_rec*Pion) in Nion.equation. Likewise, divergence.(-mu_p * grad.potential(x,t) * Pion(x,t)) can either be ConvectionTerm(coeff=-mu_p * potential.faceGrad, var=Pion) or it can be DiffusionTerm(coeff=-mu_p * Pion, var=potential). - Boundary conditions in FiPy are no-flux by default, so there's no need to constrain Pion and Nion. - Jon From: Justin Pothoof Sent: Friday, August 9, 2019 2:51 PM To: Guyer, Jonathan E. Dr. (Fed); FIPY Subject: Re: Drift Diffusion Equation Setup Hi Jon, I've been continuing to work on my problem. I've implemented the divergence of the current density mathematically into the positive charge and negative charge density equations. Now, I'm encountering issues with the boundary conditions.. the charges should be flowing towards the center of the system, and ultimately recombine, but the issues I'm seeing is that the charges are flowing out of the region. I would like to keep all of the charges confined to the system and prevent them from flowing out of the left and right boundaries. I've tried implementing neumann boundaries using Pion.faceGrad.constrain(0., where=mesh.exteriorFaces) and the negative charge version, but I'm still seeing the charges flow mostly out of the system. I would appreciate any help! Thank you, Justin import numpy as np import matplotlib.pyplot as plt from scipy import special from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer, ImplicitSourceTerm, ConvectionTerm from fipy.tools import numerix ## #''' SET-UP PARAMETERS ''' ## a1,b1,c1,a2,b2,c2 = [1.07114255e+00, 6.50014631e+05, -4.51527221e+00, 1.04633414e+00, 1.99708312e+05, -1.52479293e+00] # Parameters for sum of sines fit (Potential fit) #a = -3930.03590805 #b, c = 3049.38274411, -4.01434474 # Parameters for exponential fit (Charge Density) Not used yet q = 1.602e-19#Elementary Charge pini = 154.1581560721245/q #m^-3 nini = -134.95618729/q #m^-3 k1 = 1.8 p1 = 17 k2 = 17 p2 = 1.8 # Parameters for charge density fit (Susi's fit) l = 0.134901960784314 #Length of system in m nx = 134 #Number of cells in system dx = l/nx #Length of each cell in m x = np.linspace(0,l,nx) #Array to calculate initial values in functions below epsilon_r = 25#Relative permittivity of system epsilon = epsilon_r*8.854e-12 #Permittivity of system C/V*m kb = 1.38e-23 #J/K T = 298 #K f = kb*T/q#Volts mu_n = 1.1e-09/1 #m^2/V*s mu_p = 1.1e-09/1 #m^2/V*s Dn = f * mu_n #m^2/s Dp = f * mu_p #m^2/s k_rec = q*(mu_n+mu_p)/(2*epsilon)*10 #k_rec = 0 #pini*np.exp(a*x) #nini*np.exp(b*x+c) #Equations for exponential charge density fits (Not Used Yet) ## ##''' INITIAL CONDITION FUNCTIONS '''# ## def y01(x): """Initial positive ion charge density""" return pini*((special.gamma(k1+p1))/(special.gamma(k1)*special.gamma(p1))*((x/l)**(k1-1))*(1-(x/l))**(p1-1))/7.3572 def y02(x): Initial negative ion charge density""" return nini*((special.gamma(k2+p2))/(special.gamma(k2)*special.gamma(p2))*((x/l)**(k2-1))*(1-(x/l))**(p2-1))/7.3572 def y03(x): """Initial potential""" return a1*np.sin(b1*x+c1) + a2*np.sin(b2*x+c2) mesh = Grid1D(dx=dx, nx=nx) #Establish mesh in how many dimensions necessary ## #''' SETUP CELLVARIABLES AND EQUATIONS ''' ## #CellVariable - defines the variables that you want to solve for: '''Initial value can be established when defining the variable, or later using 'var.value =' Value defaults to
