|
|
(1) |
NN is the number of node points (elements) in the X and MM is the number of elements in the Y direction. delta is the applied current density at the top (Y =1). tf is the final time for simulation. vel is the velocity constant v in the paper. ki0 is the scaled exchange current density k in the paper. This code can be run for positive values of delta. This simulates plating. At the end of simulation, changing delta to negative values and rerunning the code will automatically used the geometry at the end of plating.
Ydatstore stores the geometry at every point in time. Phiaveadd stores the total liquid phase in the domain at any point in time.
Users can change NN, delta, tf, vel, ki0, MM just in this line and choose Edit execute worksheet to run for different design parameters.
Users can modify the call for y0proc for choosing different models.
Users can modify the call for HF to run first-order upwind, ENO2 or WENO3 methods. NN and MM should be even numbers.
| > |
NN:=100;delta:=0.1;tf:=2.0;vel:=1.0;ki0:=1.0;MM:=NN;
|
|
|
(2) |
| > |
N:=NN+2:h:=1.0/NN:M:=MM+2;Ny:=MM;
|
|
|
(3) |
Initial geometry, Model 1, semicircle in a square
| > |
y0proc1:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
|
| > |
N:=NN+2:h:=1.0/NN:w:=h/2:
|
| > |
for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=xx^2+yy^2;
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((sqrt(rr)-0.3)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
|
Model 2, square in a square
| > |
y0proc2:=proc(NN,MM,Y00)# square inside a circle, Model 2
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
|
| > |
N:=NN+2:h:=1.0/NN:w:=h/2:
|
| > |
for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if xx <=0.3 and yy<=0.3 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
|
Initial geometry, Model 3, electrodeposition problem trenches and via
| > |
y0proc3:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
|
| > |
N:=NN+2:h:=1.0/NN:w:=h/2:
|
| > |
for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
if abs(xx-0.5) >0.2 and yy<=0.5 then Y00[i,j]:=1e-9; else Y00[i,j]:=1.0:end:
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
|
Initial geometry, Model 4, Gaussian Seed at the bottom
| > |
y0proc4:=proc(NN,MM,Y00)
local i,j,xx,yy,N,h,w,M,Ny,rf,f,ff,rr;
|
| > |
N:=NN+2:h:=1.0/NN:w:=h/2:
|
| > |
for i from 2 to N-1 do for j from 2 to M-1 do
xx:=-0+(i-1/2-1)*h:yy:=-0+(j-1/2-1)*h:
rr:=0.1+0.1*exp(-500.*(xx-0.5)^2);
Y00[i,j]:=max(1e-9,0.5+0.5*tanh((yy-rr)/w/sqrt(2.0))):
od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
|
| > |
y0proc:=eval(y0proc1):#choose different models using y0proc2, etc.
|
| > |
Y00:=Matrix(1..N,1..M,datatype=float[8]):
|
| > |
evalhf(y0proc(NN,MM,Y00)):
|
| > |
if delta<0 then read("Y0data.m"):end:
|
| > |
if delta<0 then read("tdata.m"):end:
|
|
|
(4) |
| > |
p0:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):
|
|
|
(5) |
Next, boundary conditiosn at X = 0, X =1, Y = 0, Y = 1 are specified below, but these equations are not used and optimally coded inside the procedure Eqs11.
| > |
for i from 1 to 1 do for j from 1 to M do eq[i,j]:=-Phi[i,j]+Phi[i+1,j]:od:od:
|
| > |
for i from N to N do for j from 1 to M do eq[i,j]:=Phi[i-1,j]-Phi[i,j]:od:od:
|
| > |
for i from 1 to N do for j from 1 to 1 do eq[i,j]:=Phi[i,j+1]-Phi[i,j]:od:od:
|
| > |
for i from 1 to N do for j from M to M do eq[i,j]:=Phi[i,j-1]-Phi[i,j]+delta*h:od:od:
|
|
|
(6) |
Residues at different points in X and Y are coded in the Eqs11 procedure. Y0 is the input phase-field parameter (2D Matrix). The potential y is a vector to expedite the calculation of residue.
| > |
Eqs11:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,ff::Vector(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):
for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
ff[i1]:=
(Y0[i,j]+Y0[i,j+1])*(y[i1+N]-y[i1])
-(Y0[i,j]+Y0[i,j-1])*(y[i1]-y[i1-N])
+(Y0[i+1,j]+Y0[i,j])*(y[i1+1]-y[i1])
-(Y0[i,j]+Y0[i-1,j])*(y[i1]-y[i1-1])
-ki0*y[i1]*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1+N]:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:
|
The Jacobian for the residues is coded int he procedure Jac1.
| > |
Jac1:=proc(Y0::Matrix(datatype=float[8]),y::Vector(datatype=float[8]),delta::float,ki0::float,j00::Matrix(datatype=float[8]))
local i::integer,j::integer,i1::integer,h::float[8];
global N,M;
h:=1.0/(N-2):
for i from 2 to N-1 do for j from 2 to M-1 do
i1:=i+(j-1)*N:
j00[i1,i1]:=-4*Y0[i,j]-Y0[i,j+1]-Y0[i,j-1]-Y0[i+1,j]-Y0[i-1,j]-ki0*h*(1e-24+(Y0[i+1,j]-Y0[i-1,j])^2+(Y0[i,j+1]-Y0[i,j-1])^2)^(1/2):
j00[i1,i1+1]:=Y0[i+1,j]+Y0[i,j]:j00[i1,i1-1]:=Y0[i-1,j]+Y0[i,j]:
j00[i1,i1+N]:=Y0[i,j+1]+Y0[i,j]:j00[i1,i1-N]:=Y0[i,j-1]+Y0[i,j]:
od:od:
for i from 1 to 1 do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1+1]:=1.00:
#ff[i1]:=-y[i1]+y[i1+1]:
od:od:
for i from N to N do for j from 1 to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-1]:=1.00:
#ff[i1]:=-y[i1]+y[i1-1]:
od:od:
for i from 1 to N do for j from 1 to 1 do
i1:=i+(j-1)*N:
#ff[i1]:=-y[i1]+y[i1+N]:
j00[i1,i1]:=-1.0:j00[i1,i1+N]:=1.00:
od:od:
for i from 1 to N do for j from M to M do
i1:=i+(j-1)*N:
j00[i1,i1]:=-1.0:j00[i1,i1-N]:=1.00:
#ff[i1]:=-y[i1]+y[i1-N]+delta*h:
od:od:
#print(1);
#ff;
end proc:
|
| > |
Eqs11:=Compiler:-Compile(Eqs11):
|
|
|
(7) |
|
|
(8) |
First-order Upwind method is coded in the procedure UpW.
| > |
UpW:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8];
#e1:=1e-15:
h:=1.0/(N-2):
vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
vx1:=(Y00[i,j]-Y00[i-1,j])/h:
vx2:=(Y00[i+1,j]-Y00[i,j])/h:
vy1:=(Y00[i,j]-Y00[i,j-1])/h:
vy2:=(Y00[i,j+1]-Y00[i,j])/h:
if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:
if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:
end proc:
|
Second-order ENO2 method is coded in the procedure ENO2.
| > |
ENO2:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],alpha::float[8],beta::float[8];
h:=1.0/(N-2):
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
sdx:=(Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j])/h:sdy:=(Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1])/h:
vxb:=0:vxf:=0:vyb:=0:vyf:=0:
if i = 2 then sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j])/h:else sdxb:=(Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j])/h:end:
if sdx*sdxb>=0 then s1x:=1.0 else s1x:=0.0:end:
vx1:=(Y00[i,j]-Y00[i-1,j])/h+0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxb)):
if i = N-1 then sdxf:=(Y00[i+1,j]-2*Y00[i+1,j]+Y00[i,j])/h:else sdxf:=(Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j])/h:end:
if sdx*sdxf>=0 then s1x:=1.0 else s1x:=0.0:end:
vx2:=(Y00[i+1,j]-Y00[i,j])/h-0.5*signum(sdx)*s1x*min(abs(sdx),abs(sdxf)):
if j = 2 then sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1])/h:else sdyb:=(Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2])/h:end:
if sdy*sdyb>=0 then s1y:=1.0 else s1y:=0.0:end:
vy1:=(Y00[i,j]-Y00[i,j-1])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyb)):
if j = M-1 then sdyf:=(Y00[i,j+1]-2*Y00[i,j+1]+Y00[i,j])/h:else sdyf:=(Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j])/h:end:
if sdy*sdyf>=0 then s1y:=1.0 else s1y:=0.0:end:
vy2:=(Y00[i,j+1]-Y00[i,j])/h+0.5*signum(sdy)*s1y*min(abs(sdy),abs(sdyf)):
if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:
if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
F0[i,j]:=nx:
od:od:
end proc:
|
Third-order WENO3 method is coded in the procedure WENO3.
| > |
WENO3:=proc(Y00::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float,N::integer,M::integer,v0::float)
local i::integer,j::integer,h::float[8],nx::float[8],vel::float[8],vx::float[8],vy::float[8],phix::float[8],phiy::float[8],phiave::float[8],jj::integer,phixb::float[8],phixf::float[8],phixb2::float[8],phixf2::float[8],vxb::float[8],phiyb::float[8],phiyf::float[8],vxf::float[8],phiyb2::float[8],phiyf2::float[8],vyb::float[8],vyf::float[8],uf::float[8],ub::float[8],vf::float[8],vb::float[8],tt::float[8],vv0::float[8],sd::float[8],sdf::float[8],sdb::float[8],sdx::float[8],sdy::float[8],sdxb::float[8],s1x::float[8],sdxf::float[8],sdyb::float[8],sdyf::float[8],s1y::float[8],vx1::float[8],vx2::float[8],vy1::float[8],vy2::float[8],w1::float[8],w2::float[8],r1::float[8],r2::float[8],alpha::float[8],beta::float[8],e1::float[8];
e1:=1e-6:
h:=1.0/(N-2):
vv0:=v0:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
for i from 2 to N-1 do for j from 2 to M-1 do
vx:=0.0:vy:=0.0:phix:=0.0:phiy:=0.0:
phix:=(Y00[i+1,j]-Y00[i-1,j])/2/h:
phiy:=(Y00[i,j+1]-Y00[i,j-1])/2/h:
if i = 2 then sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-1,j]: else sdb:=Y00[i,j]-2*Y00[i-1,j]+Y00[i-2,j]:end:
sd:=Y00[i+1,j]-2*Y00[i,j]+Y00[i-1,j]:
if i = N-1 then sdf:=Y00[i,j]-2*Y00[i+1,j]+Y00[i+1,j]: else sdf:=Y00[i+2,j]-2*Y00[i+1,j]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vx1:=phix-0.5*w1/h*(sd-sdb):
vx2:=phix-0.5*w2/h*(sdf-sd):
if j = 2 then sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-1]:else sdb:=Y00[i,j]-2*Y00[i,j-1]+Y00[i,j-2]:end:
sd:=Y00[i,j+1]-2*Y00[i,j]+Y00[i,j-1]:
if j = M-1 then sdf:=Y00[i,j]-2*Y00[i,j+1]+Y00[i,j+1]: else sdf:=Y00[i,j+2]-2*Y00[i,j+1]+Y00[i,j]:end:
r1:=(e1+sdb^2)/(e1+sd^2):w1:=1/(1+2*r1^2):r2:=(e1+sdf^2)/(e1+sd^2):w2:=1/(1+2*r2^2):
vy1:=phiy-0.5*w1/h*(sd-sdb):
vy2:=phiy-0.5*w2/h*(sdf-sd):
if v0>=0 then
vx1:=max(vx1,0):vx2:=-min(vx2,0): else
vx1:=-min(vx1,0):vx2:=max(vx2,0): end:
if v0>=0 then
vy1:=max(vy1,0):vy2:=-min(vy2,0): else
vy1:=-min(vy1,0):vy2:=max(vy2,0): end:
nx:=sqrt(max(vx1,vx2)^2+max(vy1,vy2)^2):
#nx:=sqrt(vx1^2+vx2^2+vy1^2+vy^2):
F0[i,j]:=nx:
od:od:
end proc:
|
| > |
UpW:=Compiler:-Compile(UpW):
|
| > |
ENO2:=Compiler:-Compile(ENO2):
|
| > |
WENO3:=Compiler:-Compile(WENO3):
|
| > |
F0:=Matrix(1..N,1..M,datatype=float[8]):
|
|
|
(9) |
|
|
(10) |
| > |
PhiAdd:=proc(N::integer,Phi0::Vector(datatype=float[8]),db::Vector(datatype=float[8]))
local i::integer;
for i from 1 to N do Phi0[i]:=Phi0[i]+db[i]:od:
end proc:
|
| > |
PhiAdd:=Compiler:-Compile(PhiAdd):
|
| > |
Phi0:=Vector(1..Ntot,datatype=float[8]):
|
| > |
Phi0:=Vector(1..Ntot,datatype=float[8]):
|
| > |
j00:=Matrix(1..Ntot,1..Ntot,datatype=float[8],storage=sparse):
|
| > |
evalf(Eqs11(Y00,Phi0,delta,ki0,ff));
|
|
|
(11) |
| > |
Jac1(Y00,Phi0,delta,ki0,j00):
|
| > |
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
|
| > |
V[0]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;
|
|
|
(12) |
|
|
(13) |
| > |
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
|
| > |
dt:=min(h/vv0/vel/ki0,tf-tt);
|
|
|
(14) |
| > |
Ymid:=Matrix(1..N,1..M,datatype=float[8]):
|
| > |
phiaveAdd:=proc(Y00::Matrix(datatype=float[8]),N::integer,M::integer,Ny::integer)
local i::integer,j::integer,phiave::float;
phiave:=0.0:
for i from 2 to N-1 do for j from 2 to M-1 do phiave:=phiave+Y00[i,j]:od:od:
phiave/(N-2)/(M-2);
end proc:
|
| > |
phiaveAdd:=Compiler:-Compile(phiaveAdd):
|
|
|
(15) |
|
|
(16) |
The three stages of SSR-RK3 are stored in EFAdd, EFAdd2 and EFAdd3
| > |
EFAdd:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
#for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=Y00[i,j]-dt*F0[i,j]*Phi0[i+(j-1)*N]:od:od:
for i from 2 to N-1 do for j from 2 to M-1 do Ymid[i,j]:=max(1e-9,Y00[i,j]-dt*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:
|
| > |
EFAdd:=Compiler:-Compile(EFAdd):
|
| > |
EFAdd2:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ymid[i,j]:=max(1e-9,Y00[i,j]*3/4.+Ymid[i,j]/4.-dt/4.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Ymid[i,1]:=Ymid[i,2]:
Ymid[i,M]:=Ymid[i,M-1]:
od:
for j from 1 to M do
Ymid[1,j]:=Ymid[2,j]:
Ymid[N,j]:=Ymid[N-1,j]:
od:
end proc:
|
| > |
EFAdd2:=Compiler:-Compile(EFAdd2):
|
| > |
EFAdd3:=proc(Y00::Matrix(datatype=float[8]),Ymid::Matrix(datatype=float[8]),Phi0::Vector(datatype=float[8]),F0::Matrix(datatype=float[8]),dt::float[8],vel::float[8],ki0::float[8],N::integer,M::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Y00[i,j]:=max(1e-9,Y00[i,j]*1/3.+Ymid[i,j]*2/3.-dt*2/3.*vel*ki0*F0[i,j]*Phi0[i+(j-1)*N]):od:od:
for i from 1 to N do
Y00[i,1]:=Y00[i,2]:
Y00[i,M]:=Y00[i,M-1]:
od:
for j from 1 to M do
Y00[1,j]:=Y00[2,j]:
Y00[N,j]:=Y00[N-1,j]:
od:
end proc:
|
| > |
EFAdd3:=Compiler:-Compile(EFAdd3):
|
| > |
YdatStore:=proc(Y00::Matrix(datatype=float[8]),Ydat::Array(datatype=float[8]),N::integer,M::integer,jj::integer)
local i::integer,j::integer;
for i from 2 to N-1 do for j from 2 to M-1 do
Ydat[jj,i,j]:=Y00[i,j]:od:od:
end proc:
|
| > |
YdatStore:=Compiler:-Compile(YdatStore):
|
| > |
Phiave[0]:=phiaveAdd(Y00,N,M,Ny);
|
|
|
(17) |
|
|
(18) |
| > |
Ydat:=Array(1..Nt+1,1..N,1..M,datatype=float[8]):
|
| > |
YdatStore(Y00,Ydat,N,M,1);
|
|
|
(19) |
|
|
(20) |
Different upwind schemes can be called by assign WENO3, UpW or ENO3 scheme.
A while loop is written from t=0 to t= tf.
| > |
while tt<tf do
HF(Y00,Phi0,F0,evalf(dt),N,M,delta):
EFAdd(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd2(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Ymid,Phi0,delta,ki0,ff));
Jac1(Ymid,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
HF(Ymid,Phi0,F0,evalf(dt),N,M,delta):
EFAdd3(Y00,Ymid,Phi0,F0,dt,vel,ki0,N,M);
evalf(Eqs11(Y00,Phi0,delta,ki0,ff));
Jac1(Y00,Phi0,delta,ki0,j00):
db:=LinearAlgebra:-LinearSolve(j00,-ff,method=SparseDirect):
PhiAdd(Ntot,Phi0,db):
ii:=ii+1:#print(i);
V[ii]:=(Phi0[ntot-2*N+N/2]/2+Phi0[ntot-2*N+N/2+1]/2)+h/2*delta;#print(ii,V[ii]);
TT[ii]:=TT[ii-1]+dt:tt:=tt+dt:
YdatStore(Y00,Ydat,N,M,ii+1);
Phiave[ii]:=phiaveAdd(Y00,N,M,Ny);
vv0:=max(abs(Phi0[ntot-2*N+1]),abs(Phi0[ntot-2*N+N/2]),abs(Phi0[ntot-2*N+N/2+1]),abs(Phi0[ntot-N])):
dt:=min(h/vv0/vel/ki0,tf-tt);#print(tt,ii,dt);
#YdatStore(Y00,Ydat,N,M,i+1);
#Phiave[i]:=phiaveAdd(Y00,N,M,Ny);
end:
|
|
|
(21) |
|
|
(22) |
Voltage time curves are plotted below. Voltage is measured at X = 0.5, Y = 1.
| > |
plot([[seq([TT[ii],V[ii]],ii=0..nt)]],style=point);
|
Voltage at the end of plating, cpu time can be documented as
| > |
V[nt];time[real]()-t12;time()-t11;NN;
|
|
|
(23) |
| > |
p1:=plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]):
|
|
|
(24) |
Contour plots at t= 0 and t = 2.0 (at the of plating are given below)
| > |
plots:-display({p0});plots:-display({p1});
|
The liquid phase content as a funciton of time is plotted below
| > |
plot([[seq([TT[ii],Phiave[ii]],ii=0..nt)]],style=point);
|
| > |
plots:-surfdata(Y00,-h/2..NN*h+h/2,-h/2..MM*h+h/2,dimension=2,style=surface,colorscheme = ["Red", "Green", "Blue"]);
|
|
