Hello All,
I hope you are fine.
I am trying to solve 4 coupled PDE as below
∇.∇V=−q(Np-Ne-Nn)/ϵ
∂Ne/∂t +∇⋅(Ne.(-We)-De.∇Ne) = Se
∂Np/∂t +∇⋅(Np.(Wp)-Dp.∇Np) = Sp
∂Nn/∂t +∇⋅(Nn.(-Wn)-Dn.∇Nn) = Sn
I have tried to model the system but I faced instability I tried artificial diffusion and SUPG but also I face instability and divergence. I need your kind support please
here is the code with SUPG trail.
load "iovtk";
verbosity = 1;
//load "UMFPACK64";
real Uapp=-5500; // Applied voltage
real q=1.6022e-19; //[C] electron charge
real eps0=8.85418782e-12; // permittivity of free sace
real mup = 2.24e-4;//[m^2/(V*s)]
real mun = 2.16e-4;//[m^2/(V*s)]
real De = 0.18 ;//[m^2/s]
real Dp = 0.028e-4;//[m^2/s]
real Dn = 0.043e-4;//[m^2/s]
real beta=2e-13;//[m^3/s]
real gamma= 0.01;// secondary emission coefficient
real w1 = -1.65e7; //[V/m]
real w2 = 3.5e5; //1/m
real w3 = 1500;// [1/m]
real w4 = -2.5e6;// [V/m]
real w5 = -60.6;// [m/s]
real w6 = 2.43e-4;// [m²/(V·s)]
real w7 = -2.7e-4; //[m²/(V·s)]
real w8 = 1.9163;// [m^2/(V*s)]
real dt = 5e-13; // Time step size
real T = 2e-6; // Final Time
int numsteps = T / dt; // Number of time steps
real tt=0.;
real I;
real[int] It(numsteps);
int GND =1, needle =2, ax = 3, out1 = 4;// boundary numbering
real x0=0.03; // x-coordinate for domain
real y0=(35e-6*((29)^2)+6e-3);// y-coordinate for domain
real h=6e-3; // needle position
real yy1=0.86; //needle split point
real yy2=2.3; //needle split point
real xn1=(2*35e-6*yy1); //x-coordinate for needle split point
real yn1= (35e-6*(yy1^2)+6e-3); //y-coordinate for needle split point
real xn2=(2*35e-6*yy2); //x-coordinate for needle split point
real yn2= (35e-6*(yy2^2)+6e-3); //y-coordinate for needle split point
real w = 2*35e-6*29;
border b1 (t=0. , yy1){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};
border b2 (t=yy1 , yy2){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};
border b3 (t=yy2, 29 ){x=2*35e-6*t; y=(35e-6*(t^2)+6e-3); label=needle;};
border b4 (t=5.8e-3, 6e-3 ){x=0; y=t; label=ax;};
border b5 (t=5.6e-3, 5.8e-3){x=0; y=t; label=ax;};
border b6 (t=0 , 5.6e-3){x=0; y=t; label=ax;};
border b7 (t=5.8e-3, yn1){x=xn1; y=t;};
border b8 (t=5.6e-3, yn2){x=xn2; y=t;};
border b9 (t=0.,x0){x=t; y=0.; label=GND;};
border b10(t=0, y0){x=x0; y=t; label=out1;};
border b11(t=w,x0){x=t; y=y0; label=out1;};
mesh Th;
Th = buildmesh (
b1 (-150) //needle
+ b2 (-150) //needle
+ b3 (-50) //needle
+ b4 (-200) //ax
+ b5 (-200) //ax
+ b6 (-2000) //ax
+ b7 (150) //mesh control edge
+ b8 (300) //mesh control edge
+ b9 (50) //GND
+ b10(30) //out1
+ b11(-30) //out2
);
//finite element space
fespace VE(Th, P1); //
fespace VNe(Th, P1); //
// Define trial functions for electric potential and concentrations
VE u, v, E1, E2, Emag, E1L, E2L, EL, f, rohv;
VNe Ne, Neold, Se, S1, S2, S3, S4, We, Wex, Wey, alpha, eta, mue, tau1, dtf, Dstab, Np, Npold, Sp, Wp, Wpy, Wpx, tau2, Nn, Nnold, Sn, Wnx, Wn, Wny, tau3, v1, v2, v3;
macro grad(u) [dx(u), dy(u)] //
// initial condition for e, p, n distribution
func a = (x^2); func b = (y-h)^2; real c=2*((25e-6)^2);
real Nm = 1e16; //[m^-3] maximum density
Neold = (Nm)*exp(-(x^2/(c))-(y-0.006)^2/(c));
Npold = (Nm)*exp(-(x^2/(c))-(y-0.006)^2/(c));
Nnold = 0;
problem Poisson(u,v) = int2d(Th)(grad(u)'*grad(v))
-int2d(Th)(f*v)
+on (1,u = 0)
+on (2,u = Uapp)
;
//matrix AP;
//real [int] bP(VE.ndof);
//solve to get laplacian Electric field
f=0;
// AP = Poisson(VE,VE);
// bP = Poisson(0 ,VE);
// u[] = AP^-1*bP ;
Poisson;
E1L = -dx(u);
E2L = -dy(u);
EL = sqrt(E1L*E1L + E2L*E2L);
// plot(EL, cmm="Laplacian electric field Mag",fill=true);//, hsv=colorhsv);
problem NeL(Ne, v1) = int2d(Th)(Ne*v1/dt
+(De*(grad(Ne)'*grad(v1))
+(-Wex * dx(Ne) + -Wey * dy(Ne))*v1))
-int2d(Th)(Se*v1)
-int2d(Th)(Neold*v1/dt)
+int1d(Th,needle)(gamma*Npold*mup*Emag*v1)
+int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Ne/dt-De*(dxx(Ne)+dyy(Ne))+-Wex*dx(Ne)+-Wey*dy(Ne))) //SUPG
-int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Neold/dt)) //SUPG
-int2d(Th)((tau1*(-Wex*dx(v1)+-Wey*dy(v1)))*(Se)) //SUPG
;
problem NpL(Np, v2) = int2d(Th)(Np*v2/dt
+(Dp*(grad(Np)'*grad(v2))
+(Wpx * dx(Np) + Wpy * dy(Np))*v2))
-int2d(Th)(Sp*v2)
-int2d(Th)(Npold*v2/dt)
+on (GND, Np = 0)
+int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Np/dt-Dp*(dxx(Np)+dyy(Np))+ Wpx*dx(Np)+ Wpy*dy(Np))) //SUPG
-int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Npold/dt)) //SUPG
-int2d(Th)((tau2*( Wpx*dx(v2)+ Wpy*dy(v2)))*(Sp)) //SUPG
;
problem NnL(Nn, v3) = int2d(Th)(Nn*v3/dt
+(Dn*(grad(Nn)'*grad(v3))
+(-Wnx * dx(Nn) + -Wny * dy(Nn))*v3))
-int2d(Th)(Sn*v3)
-int2d(Th)(Nnold*v3/dt)
+on (needle, Nn = 0)
+int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Nn/dt-Dn*(dxx(Nn)+dyy(Nn))+-Wnx*dx(Nn)+-Wny*dy(Nn))) //SUPG
-int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Nnold/dt)) //SUPG
-int2d(Th)((tau3*(-Wnx*dx(v3)+-Wny*dy(v3)))*(Sn)) //SUPG
;
//matrix AE, APo, An;
//createMat(Th, AE, P1);
//createMat(Th, APo, P1);
//createMat(Th, An, P1);
//real [int] bE(VNe.ndof), bPo(VNp.ndof), bn(VNn.ndof);
real alphas=1;
VNe hr1=hTriangle;
cout << hr1[].min << endl;
cout << hr1[].max << endl;
//int step=0;
//while (tt<T){
for (int step = 0; step < numsteps; ++step){
cout<< "time:" << tt << endl;
cout<< "step #:" << step << endl;
rohv = q*(Npold-Neold-Nnold);
f = -rohv/eps0;
// AP=Poisson(VE,VE);
// bP=Poisson(0,VE);
// u[] = AP^-1*bP ;
// cout << "Volt Max " << u[].max << endl ;
// cout << "Volt Min " << u[].min << endl << endl ;
Poisson;
E1 = -dx(u); E2 = -dy(u); Emag = sqrt(E1*E1 + E2*E2);
mue = w8/((Emag)^0.25);
alpha = w2*exp(w1/Emag);
eta = w3*exp(w4/Emag);
We = mue*Emag;
Wp = mup*Emag;
Wn = mun*Emag;
Wex = E1*mue;
Wey = E2*mue;
Wpx = E1*mup;
Wpy = E2*mup;
Wnx = E1*mun;
Wny = E2*mun;
tau1= alphas*hr1/2/sqrt(Wex^2+Wey^2);
tau2= alphas*hr1/2/sqrt(Wpx^2+Wpy^2);
tau3= alphas*hr1/2/sqrt(Wnx^2+Wny^2);
// dtf = hr/We;
// dt = dtf[].min;
// cout<< "Time step #:" << dt << endl;
// cout<< "Calculated Time step #:" << dtf[].min << endl;
S1 = Neold*alpha*abs(We); // [1/m^3.s] Ionizatoin
S2 = Neold*eta*abs(We); // [1/m^3.s] Attachment
S3 = Neold*Npold*beta; // [1/m^3.s] Recombination between electrons and positive ions
S4 = Nnold*Npold*beta; // [1/m^3.s] Recombination between positive ions and negative ions
Se = S1-S2-S3;
Sp = S1-S3-S4;
Sn = S2-S4;
// Dstab = De+(0.25*hr*We);
NpL;
// AE=NeL(VNe,VNe);
// bE=NeR(0,VNe);
// Ne[] = AE^-1*bE ;
// cout << "Electron Max " << Ne[].max << endl ;
// cout << "Electron Min " << Ne[].min << endl << endl ;
NeL;
// APo=NpL(VNp,VNp);
// bPo=NeR(0,VNp);
// Np[] = APo^-1*bPo ;
// cout << "Positive Max " << Np[].max << endl ;
// cout << "Positive Min " << Np[].min << endl << endl ;
NnL;
// An=NnL(VNn,VNn);
// bn=NnR(0,VNn);
// Nn[] = An^-1*bn ;
// cout << "Negative Max " << Nn[].max << endl ;
// cout << "Negative Min " << Nn[].min << endl << endl ;
//Peclet Number
// VNe hK = hTriangle ;
VNe PeK = We*hr1/2./De ;
cout << "Pe max= " << PeK[].max << endl ;
cout << "Pe min= " << PeK[].min << endl ;
if (step % 50 == 0)
{
savevtk("SUPG 3SP dt5e-13-Newmesh2/TNS "+step+".vtu", Th, Emag, Ne, Np, Nn, PeK);
}
for (int i=0; i<VNe.ndof; ++i){if(Ne[][i]<0.)Ne[][i]=0.;}
for (int i=0; i<VNe.ndof; ++i){if(Np[][i]<0.)Np[][i]=0.;}
for (int i=0; i<VNe.ndof; ++i){if(Nn[][i]<0.)Nn[][i]=0.;}
I = int2d(Th)(2*pi*x*q*(mue*Ne + mup*Np + mun*Nn)*((E1L*E1)+(E2L*E2)));
I = I/abs(Uapp);
cout << I << endl;
It(step)=I;
// Tt(step)=tt;
{
ofstream file("SUPG 3SP-dt-5e-13.txt");
file << It << endl;
}
Neold = Ne;
Npold = Np;
Nnold = Nn;
tt += dt;
// ++step;
}
