Коновалов Учебно-методическое пособие по курсу 2007

+a1k1 (α21α11 + β21β11 ).

Вматематической постановке (36), разностную схему 2-го порядка точности по координате можно получить, используя основную разностную схему (43) в граничных узлах при x = 0 , y = 0

математической постановки задачи (31) (схемы (43)). На графиках в программе последовательно представлены уровни поверхности ледникового покрова в различные моменты времени в плоскости ледораздела y =1/ 2 , стационарные решения в плоскостях x =1/ 2 и

y =1/ 2 , изменение толщины ледникового покрова на вершине со

61

временем, поверхность ледника. Время накопления ледниковой массы и толщина покрова в центральной части соответствуют результатам, полученным ранее конечно-разностным методом,

наблюдается немонотонность изменения численного решения в отличие от решения, полученного конечно-разностным методом. Впоследствии решение выходит на монотонно изменяющиеся стационарные значения, если выполнено условие устойчивости

γ = ∆τ / ∆2 <α , где α ≈ 0.4 . Устойчивость схемы исследовалась в

математической постановке задачи (36) для ¼ области ледника. Результаты исследования представлены в таб. 4.

>#===========================================

>#=== 2D Ice Sheet Model (Finite-Element Method) =====

>#===========================================

>restart:

>#--------------------

>N_Glen:=3;

N_Glen := 3

> a:=0.3;

a := 0.3

> dens:=910;

dens := 910

62

> g:=9.8;

g := 9.8

> Ao:=1e-16;

Ao := 0.1 10-15

> Lo:=2*10^5;

Lo := 200000

>

Zo:=(a*(N_Glen+2)*Lo^(N_Glen+1)/(2*Ao*(dens*g)^N_Glen))^(1/(

2*(N_Glen+1)));

Zo := 1424.122524

> To:=Zo/a;

To := 4747.075080

>

>#--------------------------------------------------------------

>Nt:=1500;

Nt := 1500

> No:=21;

No := 21

> dl:=1/(No-1);

dl := 201

> dt:=0.001;

dt := 0.001

>

> a_e:=dl^2/2;

a_e := 8001

>

>#--------------------------------------------------------------

>N1:=100: N2:=200: N3:=300: N4:=500: N5:=800:

>#--------------------------------------------------------------

>S:=array(1..No,1..No):

>#--------------------------

>H_Ice_Summit:=array(0..Nt):

>#--------------------------------------------------------------

63

>for i from 1 to No do for j from 1 to No do S[i,j]:=0: end do: end do:

>#--------------------------------------------------------------

>H_Ice_Summit[0]:=S[round(No/2),round(No/2)]:

>#--------------------------------------------------------------

>M:=Matrix(No^2);

>

>#---------------------------------------------------------------

>M_abc:=Matrix(3);

M_abc[2,1]:=i*dl: M_abc[2,2]:=j*dl: M_abc[2,3]:=1: M_abc[3,1]:=(i-1)*dl: M_abc[3,2]:=j*dl: M_abc[3,3]:=1: V_abc[1]:=1: V_abc[2]:=0: V_abc[3]:=0: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[1,1]:=L_abc[1]: b[1,1]:=L_abc[2]: c[1,1]:=L_abc[3]:

64

V_abc[1]:=0: V_abc[2]:=1: V_abc[3]:=0: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[1,2]:=L_abc[1]: b[1,2]:=L_abc[2]: c[1,2]:=L_abc[3]: V_abc[1]:=0: V_abc[2]:=0: V_abc[3]:=1: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[1,3]:=L_abc[1]: b[1,3]:=L_abc[2]: c[1,3]:=L_abc[3]:

>M_abc[1,1]:=(i-1)*dl: M_abc[1,2]:=(j-1)*dl: M_abc[1,3]:=1: M_abc[2,1]:=i*dl: M_abc[2,2]:=j*dl: M_abc[2,3]:=1: M_abc[3,1]:=i*dl: M_abc[3,2]:=(j-1)*dl: M_abc[3,3]:=1: V_abc[1]:=1: V_abc[2]:=0: V_abc[3]:=0: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[6,1]:=L_abc[1]: b[6,1]:=L_abc[2]: c[6,1]:=L_abc[3]: V_abc[1]:=0: V_abc[2]:=1: V_abc[3]:=0: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[6,2]:=L_abc[1]: b[6,2]:=L_abc[2]: c[6,2]:=L_abc[3]: V_abc[1]:=0: V_abc[2]:=0: V_abc[3]:=1: L_abc:=LinearAlgebra[LinearSolve](M_abc,V_abc): a[6,3]:=L_abc[1]: b[6,3]:=L_abc[2]: c[6,3]:=L_abc[3]:

>#-------------------------------------------------------------

>Dlk:=array(1..6,1..3,1..3):

>#-------------------------------------------------------------

>for p from 1 to 6 do

for l from 1 to 3 do for k from 1 to 3 do

Dlk[p,l,k]:=a[p,l]*a[p,k]+b[p,l]*b[p,k]: end do: end do: end do:

>

>#-------------------------------------------------------------

>for m from 1 to Nt do

for i from 1 to No do M[(i-1)*No+1,(i-1)*No+1]:=1: M[(i-1)*No+No,(i-1)*No+No]:=1:

65

end do:

for j from 2 to No-1 do M[j,j]:=1: M[(No-1)*No+j,(No-1)*No+j]:=1: end do:

for i from 2 to No-1 do for j from 2 to No-1 do

K_eff_1:=((S[i,j]+S[i+1,j+1]+S[i,j+1])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[1 ,1,1]+S[i+1,j+1]^2*Dlk[1,2,2]+S[i,j+1]^2*Dlk[1,3,3]+2*S[i,j]*S[i+1,j+ 1]*Dlk[1,1,2]+2*S[i,j]*S[i,j+1]*Dlk[1,1,3]+2*S[i+1,j+1]*S[i,j+1]*Dlk[ 1,2,3])^((N_Glen-1)/2):

K_eff_2:=((S[i,j]+S[i- 1,j]+S[i,j+1])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[2,1,1]+S[i- 1,j]^2*Dlk[2,2,2]+S[i,j+1]^2*Dlk[2,3,3]+2*S[i,j]*S[i- 1,j]*Dlk[2,1,2]+2*S[i,j]*S[i,j+1]*Dlk[2,1,3]+2*S[i- 1,j]*S[i,j+1]*Dlk[2,2,3])^((N_Glen-1)/2):

K_eff_3:=((S[i,j]+S[i-1,j]+S[i-1,j- 1])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[3,1,1]+S[i-1,j]^2*Dlk[3,2,2]+S[i-1,j- 1]^2*Dlk[3,3,3]+2*S[i,j]*S[i-1,j]*Dlk[3,1,2]+2*S[i,j]*S[i-1,j- 1]*Dlk[3,1,3]+2*S[i-1,j]*S[i-1,j-1]*Dlk[3,2,3])^((N_Glen-1)/2):

K_eff_4:=((S[i,j]+S[i,j-1]+S[i-1,j- 1])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[4,1,1]+S[i,j-1]^2*Dlk[4,2,2]+S[i-1,j- 1]^2*Dlk[4,3,3]+2*S[i,j]*S[i,j-1]*Dlk[4,1,2]+2*S[i,j]*S[i-1,j- 1]*Dlk[4,1,3]+2*S[i,j-1]*S[i-1,j-1]*Dlk[4,2,3])^((N_Glen-1)/2):

K_eff_5:=((S[i,j]+S[i,j- 1]+S[i+1,j])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[5,1,1]+S[i,j- 1]^2*Dlk[5,2,2]+S[i+1,j]^2*Dlk[5,3,3]+2*S[i,j]*S[i,j- 1]*Dlk[5,1,2]+2*S[i,j]*S[i+1,j]*Dlk[5,1,3]+2*S[i,j- 1]*S[i+1,j]*Dlk[5,2,3])^((N_Glen-1)/2):

66

K_eff_6:=((S[i,j]+S[i+1,j+1]+S[i+1,j])/3)^(N_Glen+2)*(S[i,j]^2*Dlk[6 ,1,1]+S[i+1,j+1]^2*Dlk[6,2,2]+S[i+1,j]^2*Dlk[6,3,3]+2*S[i,j]*S[i+1,j+ 1]*Dlk[6,1,2]+2*S[i,j]*S[i+1,j]*Dlk[6,1,3]+2*S[i+1,j+1]*S[i+1,j]*Dlk[ 6,2,3])^((N_Glen-1)/2):

M[(i-1)*No+j,(i- 1)*No+j]:=a_e*(1/dt+K_eff_1*Dlk[1,1,1]+K_eff_2*Dlk[2,1,1]+K_eff_ 3*Dlk[3,1,1]+K_eff_4*Dlk[4,1,1]+K_eff_5*Dlk[5,1,1]+K_eff_6*Dlk[6 ,1,1]):

M[(i-1)*No+j,(i- 1)*No+j+1]:=a_e*(1/(6*dt)+K_eff_1*Dlk[1,3,1]+K_eff_2*Dlk[2,3,1]):

M[(i-1)*No+j,(i- 2)*No+j]:=a_e*(1/(6*dt)+K_eff_2*Dlk[2,2,1]+K_eff_3*Dlk[3,2,1]):

M[(i-1)*No+j,(i-2)*No+j- 1]:=a_e*(1/(6*dt)+K_eff_3*Dlk[3,3,1]+K_eff_4*Dlk[4,3,1]):

M[(i-1)*No+j,(i-1)*No+j- 1]:=a_e*(1/(6*dt)+K_eff_4*Dlk[4,2,1]+K_eff_5*Dlk[5,2,1]):

M[(i- 1)*No+j,i*No+j]:=a_e*(1/(6*dt)+K_eff_5*Dlk[5,3,1]+K_eff_6*Dlk[6, 3,1]):

M[(i- 1)*No+j,i*No+j+1]:=a_e*(1/(6*dt)+K_eff_1*Dlk[1,2,1]+K_eff_6*Dlk[ 6,2,1]):

V[(i-1)*No+j]:=a_e*(2+S[i,j]/dt+S[i+1,j+1]/(6*dt)+S[i,j+1]/(6*dt)+S[i- 1,j]/(6*dt)+S[i-1,j-1]/(6*dt)+S[i,j-1]/(6*dt)+S[i+1,j]/(6*dt)):

end do: end do:

H_Ice_Summit[m]:=S[round(No/2),round(No/2)]:

Surf:=LinearAlgebra[LinearSolve](M,V):

67

>L_S:=[seq(S[i,round(No/2)],i=1..No)]: Surf:=spline(L_x,L_S,x,1): gr6:=plot(Surf,x=0..1,color=blue,thickness=3,linestyle=1,legend="Nt =1500"):

>plots[display]([gr1,gr2,gr3,gr4,gr5,gr6],axes=boxed,title="Ice Cap Surface",labels=["x *200, km","Elevation, non-dimensional value"],labeldirections=[HORIZONTAL,VERTICAL]);

>#-------------------------------------------------------------

>L_S:=[seq(S[i,round(No/2)],i=1..No)]: Surf:=spline(L_x,L_S,x,1): gr1:=plot(Surf,x=0..1,color=red,thickness=3,linestyle=3,legend="y=1/ 2"):

>L_S:=[seq(S[round(No/2),j],j=1..No)]: Surf:=spline(L_y,L_S,x,1): gr2:=plot(Surf,x=0..1,color=blue,thickness=3,linestyle=1,legend="x= 1/2"):

69

> plots[display]([gr1,gr2],axes=boxed,title="Ice Cap

Surface",labels=["x,y *200, km","Elevation, non-dimensional value"],labeldirections=[HORIZONTAL,VERTICAL]);

>#-------------------------------------------------------------

>L_t:=[seq(m*dt,m=0..Nt)]:

>L_S:=[seq(H_Ice_Summit[m],m=0..Nt)]:

>Surf:=spline(L_t,L_S,t,1):

>plot(Surf,t=0..Nt*dt,axes=boxed,title="Ice Thickness at the Ice Cap Summit", labels=["time, non-dimensional value","Elevation, nondimensional value"],labeldirections=[HORIZONTAL,VERTICAL],color=black,t hickness=3,linestyle=1);

70