Skip to content

Commit

Permalink
Added FieldEmissionBC. Not converging. (Issue #3).
Browse files Browse the repository at this point in the history
  • Loading branch information
Alex Lindsay committed Sep 3, 2016
1 parent 09c2955 commit 7c4fce6
Show file tree
Hide file tree
Showing 9 changed files with 15,682 additions and 0 deletions.
62 changes: 62 additions & 0 deletions include/bcs/FieldEmissionBC.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
#ifndef FIELDEMISSIONBC_H
#define FIELDEMISSIONBC_H

#include "IntegratedBC.h"

class FieldEmissionBC;

template<>
InputParameters validParams<FieldEmissionBC>();

class FieldEmissionBC : public IntegratedBC
{
public:

FieldEmissionBC(const InputParameters & parameters);

protected:
virtual Real computeQpResidual();
virtual Real computeQpJacobian();
virtual Real computeQpOffDiagJacobian(unsigned int jvar);

Real _r_units;
Real _r;

// Coupled variables

const VariableGradient & _grad_potential;
unsigned int _potential_id;
const VariableValue & _mean_en;
unsigned int _mean_en_id;
MooseVariable & _ip_var;
const VariableValue & _ip;
const VariableGradient & _grad_ip;
unsigned int _ip_id;

const MaterialProperty<Real> & _muem;
const MaterialProperty<Real> & _d_muem_d_actual_mean_en;
const MaterialProperty<Real> & _massem;
const MaterialProperty<Real> & _e;
const MaterialProperty<Real> & _sgnip;
const MaterialProperty<Real> & _muip;
const MaterialProperty<Real> & _Dip;
const MaterialProperty<Real> & _work_function;
const MaterialProperty<Real> & _field_enhancement;

Real _a;
Real _v_thermal;
RealVectorValue _ion_flux;
Real _n_gamma;
Real _d_v_thermal_d_u;
Real _d_v_thermal_d_mean_en;
RealVectorValue _d_ion_flux_d_potential;
RealVectorValue _d_ion_flux_d_ip;
Real _d_n_gamma_d_potential;
Real _d_n_gamma_d_ip;
Real _d_n_gamma_d_u;
Real _d_n_gamma_d_mean_en;
Real _actual_mean_en;

};

#endif //FIELDEMISSIONBC_H
4 changes: 4 additions & 0 deletions include/materials/Gas.h
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ class Gas : public Material
std::string _potential_units;
Real _voltage_scaling;
Real _user_se_coeff;
Real _user_work_function;
Real _user_field_enhancement;
Real _user_T_gas;
Real _user_p_gas;
bool _use_moles;
Expand All @@ -68,6 +70,8 @@ class Gas : public Material
MaterialProperty<Real> & _massGas;
MaterialProperty<Real> & _massArp;
MaterialProperty<Real> & _se_coeff;
MaterialProperty<Real> & _work_function;
MaterialProperty<Real> & _field_enhancement;
MaterialProperty<Real> & _se_energy;
MaterialProperty<Real> & _ElectronTotalFluxMag;
MaterialProperty<Real> & _ElectronTotalFluxMagSizeForm;
Expand Down
2 changes: 2 additions & 0 deletions src/base/ZapdosApp.C
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@
#include "PenaltyCircuitPotential.h"
#include "CircuitDirichletPotential.h"
#include "SecondaryElectronBC.h"
#include "FieldEmissionBC.h"
#include "HagelaarEnergyBC.h"
#include "HagelaarIonAdvectionBC.h"
#include "HagelaarIonDiffusionBC.h"
Expand Down Expand Up @@ -197,6 +198,7 @@ ZapdosApp::registerObjects(Factory & factory)
registerBoundaryCondition(PenaltyCircuitPotential);
registerBoundaryCondition(CircuitDirichletPotential);
registerBoundaryCondition(SecondaryElectronBC);
registerBoundaryCondition(FieldEmissionBC);
registerBoundaryCondition(HagelaarIonAdvectionBC);
registerBoundaryCondition(HagelaarIonDiffusionBC);
registerBoundaryCondition(HagelaarElectronBC);
Expand Down
189 changes: 189 additions & 0 deletions src/bcs/FieldEmissionBC.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,189 @@
#include "FieldEmissionBC.h"

template<>
InputParameters validParams<FieldEmissionBC>()
{
InputParameters params = validParams<IntegratedBC>();
params.addRequiredParam<Real>("r", "The reflection coefficient");
params.addRequiredCoupledVar("potential","The electric potential");
params.addRequiredCoupledVar("mean_en", "The mean energy.");
params.addRequiredCoupledVar("ip", "The ion density.");
params.addRequiredParam<Real>("position_units", "Units of position.");
return params;
}

FieldEmissionBC::FieldEmissionBC(const InputParameters & parameters) :
IntegratedBC(parameters),

_r_units(1. / getParam<Real>("position_units")),
_r(getParam<Real>("r")),

// Coupled Variables
_grad_potential(coupledGradient("potential")),
_potential_id(coupled("potential")),
_mean_en(coupledValue("mean_en")),
_mean_en_id(coupled("mean_en")),
_ip_var(*getVar("ip",0)),
_ip(coupledValue("ip")),
_grad_ip(coupledGradient("ip")),
_ip_id(coupled("ip")),

_muem(getMaterialProperty<Real>("muem")),
_d_muem_d_actual_mean_en(getMaterialProperty<Real>("d_muem_d_actual_mean_en")),
_massem(getMaterialProperty<Real>("massem")),
_e(getMaterialProperty<Real>("e")),
_sgnip(getMaterialProperty<Real>("sgn" + _ip_var.name())),
_muip(getMaterialProperty<Real>("mu" + _ip_var.name())),
_Dip(getMaterialProperty<Real>("diff" + _ip_var.name())),
_work_function(getMaterialProperty<Real>("work_function")),
_field_enhancement(getMaterialProperty<Real>("field_enhancement")),
_a(0.5),
_v_thermal(0),
_ion_flux(0, 0, 0),
_n_gamma(0),
_d_v_thermal_d_u(0),
_d_v_thermal_d_mean_en(0),
_d_ion_flux_d_potential(0, 0, 0),
_d_ion_flux_d_ip(0, 0, 0),
_d_n_gamma_d_potential(0),
_d_n_gamma_d_ip(0),
_d_n_gamma_d_u(0),
_d_n_gamma_d_mean_en(0),
_actual_mean_en(0)
{}

Real
FieldEmissionBC::computeQpResidual()
{
Real a;
Real b;
Real v;
Real f;
Real je;
Real F;

if ( _normals[_qp] * -1.0 * -_grad_potential[_qp] > 0.0) {
_a = 1.0;
}
else {
_a = 0.0;
}

_v_thermal = std::sqrt(8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]));

// Fowler-Nordheim
// je = (a / wf) * F^2 * exp(-v(f) * b * wf^1.5 / F)
// a = 1.541434E-6 A eV/V^2
// b = 6.830890E9 V/m-eV^1.5
// v(f) = 1 - f + (f/6)*ln(f)
// f = (1.439964E9 eV^2 m/V)*(F/wf^2)

F = _field_enhancement[_qp] * _normals[_qp] * -_grad_potential[_qp] * _r_units;

a = 1.541434E-6; // A eV/V^2
b = 6.830890E9; // V/m-eV^1.5

f = (1.439964E9)*(F / pow(_work_function[_qp], 2) );
v = 1 - f + (f/6)*std::log(f);

je = (a / (_work_function[_qp])) * pow( F , 2) * std::exp(v * b * pow(_work_function[_qp], 1.5) / F);

return _test[_i][_qp] * (je / _e[_qp]);
}

Real
FieldEmissionBC::computeQpJacobian()
{
if ( _normals[_qp] * -1.0 * -_grad_potential[_qp] > 0.0) {
_a = 1.0;
}
else {
_a = 0.0;
}

_actual_mean_en = std::exp(_mean_en[_qp] - _u[_qp]);

_v_thermal = std::sqrt(8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]));
_d_v_thermal_d_u = 0.5 / _v_thermal * 8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]) * -_phi[_j][_qp];

return 0.;
}

Real
FieldEmissionBC::computeQpOffDiagJacobian(unsigned int jvar)
{
Real a;
Real b;
Real v;
Real f;
Real je;
Real F;
Real _E_Flux;
Real _dv_dF;


// Fowler-Nordheim
// je = (a / wf) * F^2 * exp(-v(f) * b * wf^1.5 / F)
// a = 1.541434E-6 A eV/V^2
// b = 6.830890E9 V/m-eV^1.5
// v(f) = 1 - f + (f/6)*ln(f)
// f = (1.439964E9 eV^2 m/V)*(F/wf^2)

F = _field_enhancement[_qp] * _normals[_qp] * -_grad_potential[_qp] * _r_units;

a = 1.541434E-6; // A eV/V^2
b = 6.830890E9; // V/m-eV^1.5

f = (1.439964E9)*(F / pow(_work_function[_qp], 2) );
v = 1 - f + (f/6)*std::log(f);

je = (a / (_work_function[_qp])) * pow( F , 2) * std::exp(v * b * pow(_work_function[_qp], 1.5) / F);
_E_Flux = je / _e[_qp];

if (jvar == _potential_id)
{
if ( _normals[_qp] * -1.0 * -_grad_potential[_qp] > 0.0)
_a = 1.0;
else
_a = 0.0;

_v_thermal = std::sqrt(8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]));

_dv_dF = (-(5.0/6.0) + (1.0/6.0) * std::log(f)) * (1.439964E9 / pow(_work_function[_qp], 2) ) ;

return - _test[_i][_qp] * ( (2.0 / F) - ( b * pow(_work_function[_qp], 1.5) * v / pow(F,2) ) + ( b * pow(_work_function[_qp], 1.5) * _dv_dF / F ) ) * ( _E_Flux ) * (-_grad_phi[_j][_qp] * _normals[_qp]) * _r_units ;
}

else if (jvar == _mean_en_id)
{
if ( _normals[_qp] * -1.0 * -_grad_potential[_qp] > 0.0) {
_a = 1.0;
}
else {
_a = 0.0;
}
_v_thermal = std::sqrt(8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]));
_d_v_thermal_d_mean_en = 0.5 / _v_thermal * 8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]) * _phi[_j][_qp];
_actual_mean_en = std::exp(_mean_en[_qp] - _u[_qp]);

return 0.;
}

else if (jvar == _ip_id)
{
if ( _normals[_qp] * -1.0 * -_grad_potential[_qp] > 0.0) {
_a = 1.0;
}
else {
_a = 0.0;
}

_v_thermal = std::sqrt(8 * _e[_qp] * 2.0 / 3 * std::exp(_mean_en[_qp] - _u[_qp]) / (M_PI * _massem[_qp]));
_d_ion_flux_d_ip = _sgnip[_qp] * _muip[_qp] * -_grad_potential[_qp] * _r_units * std::exp(_ip[_qp]) * _phi[_j][_qp] - _Dip[_qp] * std::exp(_ip[_qp]) * _grad_phi[_j][_qp] * _r_units - _Dip[_qp] * std::exp(_ip[_qp]) * _phi[_j][_qp] * _grad_ip[_qp] * _r_units;

return 0;
}

else
return 0.0;
}
8 changes: 8 additions & 0 deletions src/materials/Gas.C
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ InputParameters validParams<Gas>()
params.addRequiredParam<bool>("use_moles", "Whether to use units of moles as opposed to # of molecules.");

params.addParam<Real>("user_se_coeff", 0.15, "The secondary electron emission coefficient.");
params.addParam<Real>("user_work_function", 4.55, "The work function.");
params.addParam<Real>("user_field_enhancement", 55, "The field enhancement factor.");
params.addParam<Real>("user_T_gas", 300, "The gas temperature in Kelvin.");
params.addParam<Real>("user_p_gas", 1.01e5, "The gas pressure in Pascals.");

Expand All @@ -33,6 +35,8 @@ Gas::Gas(const InputParameters & parameters) :
_ramp_trans_coeffs(getParam<bool>("ramp_trans_coeffs")),
_potential_units(getParam<std::string>("potential_units")),
_user_se_coeff(getParam<Real>("user_se_coeff")),
_user_work_function(getParam<Real>("user_work_function")),
_user_field_enhancement(getParam<Real>("user_field_enhancement")),
_user_T_gas(getParam<Real>("user_T_gas")),
_user_p_gas(getParam<Real>("user_p_gas")),
_use_moles(getParam<bool>("use_moles")),
Expand All @@ -56,6 +60,8 @@ Gas::Gas(const InputParameters & parameters) :
_massGas(declareProperty<Real>("massGas")),
_massArp(declareProperty<Real>("massArp")),
_se_coeff(declareProperty<Real>("se_coeff")),
_work_function(declareProperty<Real>("work_function")),
_field_enhancement(declareProperty<Real>("field_enhancement")),
_se_energy(declareProperty<Real>("se_energy")),
_ElectronTotalFluxMag(declareProperty<Real>("ElectronTotalFluxMag")),
_ElectronTotalFluxMagSizeForm(declareProperty<Real>("ElectronTotalFluxMagSizeForm")),
Expand Down Expand Up @@ -193,6 +199,8 @@ Gas::computeQpProperties()
_n_gas[_qp] = _p_gas[_qp] / (_k_boltz[_qp] * _T_gas[_qp]);

_se_coeff[_qp] = _user_se_coeff;
_work_function[_qp] = _user_work_function;
_field_enhancement[_qp] = _user_field_enhancement;
_se_energy[_qp] = 2. * 3. / 2.; // Emi uses 2 Volts coming off the wall (presumably for Te). Multiply by 3/2 to get mean_en
_e[_qp] = 1.6e-19;
_eps[_qp] = 8.85e-12;
Expand Down
Loading

0 comments on commit 7c4fce6

Please sign in to comment.