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

si, j−1

= si, j+1 , ki, j−1 = ki, j

(на

границе y = 0 ).

Тогда граничные

условия будут иметь вид:

m+1 1

2

m

m

m+1

1

m

m

s1, j

+

(k1, j

+ k1, j−1 )

+ s2, j −

(k1, j

+ k1, j−1 ) +

∆τ

∆2

∆2

m+1

1

m

m+1

1

m

s1,mj

+ s1, j+1

−

k1, j

+ s1, j−1

−

k1, j−1

=1

+

,i =

1, j = 2...N −1;

∆2

∆2

∆τ

m+1 1

2

m

m

m+1

1

m

m

si,1

+

(ki−1,1 + ki,1 ) + si,2 −

(ki−1,1

+ ki,1 )

+

∆τ

∆2

∆2

1

1

sim,1

(37)

m+1

m

m+1

m

+ si+1,1

−

ki,1

+ si−1,1

−

ki−1,1

=1 +

, j =1,i = 2...N −1;

∆2

∆2

∆τ

m+1

1

4

m

m+1

2

m

s1,1

+

k1,1

+ s2,1

−

k1,1

+

∆2

∆2

∆τ

m+1

2

m

s1,1m

+ s1,2

−

k1,1

=

1 +

, i

=1, j =1.

∆2

∆τ

Численное решение задачи (35) реализовано с помощью математического пакета Maple.

На графиках в программе (в конце параграфа) последовательно представлены: расположение ненулевых элементов матрицы системы, уровни поверхности ледникового покрова в различные моменты времени в плоскостях ледораздела ( x = 0 и y = 0 ), изменение толщины ледникового покрова на вершине со

временем, форма поверхности ледника. Время накопления ледниковой массы (время выхода решения задачи (25) на стационарное решение (22)) для выбранных значений входных

данных составляет ≈ 6 103 лет – значение представлено в конце программы.

Результаты получены как для постоянной во времени

аккумуляции a(t) =1 , так	и для	случая периодически
осциллирующей аккумуляции	a(t) =1+sin	6π t	.
		T

41

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

периода изменений климата Земли: 100, 41, 23, и 19 ×103 лет (циклы Миланковича). В программе период изменения аккумуляции

составляет величину ≈ 20 103 лет, что соответствует одному из циклов Миланковича. Период и фаза колебаний толщины льда соответствуют входным данным для аккумуляции.

Аналогично одномерным случаям разностная схема (35), устойчива при условии γ = ∆τ / ∆2 <α , где α ≈1.2 (см. графики в

конце программы). Числа узлов разбиения по координате, начиная с которых для заданного шага по времени проявляется нарушение сходимости численного решения, представлены в табл.3.

Таблица 3. Параметры шагов программы
∆τ	2 10−2	10−2	5 10−3	2 10−3	10−3	5 10−4
Nmax	11	15	18	26	36	51
∆τ / ∆2	2	1.96	1.45	1.25	1.23	1.25

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

>#== 2D-Ice Cap Model, Shallow Ice Approximation =====

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

>restart:

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

>

> N_Glen:=3;

N_Glen := 3

> a:=0.3;

a := 0.3

> dens:=910;

dens := 910

> g:=9.8;

g := 9.8

42

> Ao:=1e-16;

Ao := 0.1 10-15

> Lo:=10^5;

Lo := 100000

>

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

2*(N_Glen+1)));

Zo := 1007.006694

> To:=Zo/a;

To := 3356.688980

>

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

>Nt:=18000;

Nt := 18000

> No:=11;

No := 11

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

dl := 101

> dt:=0.001;

dt := 0.001

>

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

>K_eff:=array(1..No-1,1..No-1):

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

>S1:=array(1..No,1..No): S2:=array(1..No,1..No):

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

>H_Ice_Summit1:=array(0..Nt): H_Ice_Summit2:=array(0..Nt):

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

>

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

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

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

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

43

>H_Ice_Summit1[0]:=S1[1,1]:

>H_Ice_Summit2[0]:=S2[1,1]:

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

>M:=Matrix(No^2);

	121 x 121 Matrix
	Data Type: anything
M :=	Storage: rectangular
	Order: Fortran_order
>

> V:=Vector(No^2);	121 Element Column Vector
	Data Type: anything
V :=	Storage: rectangular
	Order: Fortran_order
>
> #--------------------------
>

> for m from 1 to Nt do

...

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

K_eff[i,j]:=(((S2[i,j]+S2[i+1,j]+S2[i+1,j+1]+S2[i,j+1])/4)^(N_Glen+2)) *((1/(4*dl^2))*((S2[i+1,j]-S2[i,j]+S2[i+1,j+1]-S2[i,j+1])^2+(S2[i,j+1]- S2[i,j]+S2[i+1,j+1]-S2[i+1,j])^2))^((N_Glen-1)/2):

end do: end do:

M[1,1]:=1/dt+(4/dl^2)*K_eff[1,1]:

M[1,No+1]:=(-2/(dl^2))*K_eff[1,1]:

M[1,2]:=(-2/(dl^2))*K_eff[1,1]:

for i from 2 to No-1 do M[(i-1)*No+1,(i-1)*No+1]:=1/dt+(2/dl^2)*(K_eff[i-1,1]+K_eff[i,1]): M[(i-1)*No+1,i*No+1]:=(-1/(dl^2))*K_eff[i,1]: M[(i-1)*No+1,(i-2)*No+1]:=(-1/(dl^2))*K_eff[i-1,1]: M[(i-1)*No+1,(i-1)*No+2]:=(-1/dl^2)*(K_eff[i-1,1]+K_eff[i,1]):

44

end do:

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

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

for j from 1 to No-1 do 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

M[(i-1)*No+j,(i-1)*No+j]:=(1/dt+(1/dl^2)*(K_eff[i,j]+K_eff[i- 1,j]+K_eff[i,j-1]+K_eff[i-1,j-1])): M[(i-1)*No+j,(i-1)*No+j+1]:=(-1/(2*dl^2))*(K_eff[i-1,j]+K_eff[i,j]): M[(i-1)*No+j,(i-1)*No+j-1]:=(-1/(2*dl^2))*(K_eff[i-1,j-1]+K_eff[i,j- 1]):

M[(i-1)*No+j,i*No+j]:=(-1/(2*dl^2))*(K_eff[i,j]+K_eff[i,j-1]): M[(i-1)*No+j,(i-2)*No+j]:=(-1/(2*dl^2))*(K_eff[i-1,j]+K_eff[i-1,j-1]): end do:

end do:

V[1]:=S2[1,1]/dt+1+sin(6*evalf(Pi)*(m-1)/(Nt-1)):

for i from 2 to No-1 do V[(i-1)*No+1]:=S2[i,1]/dt+1+sin(6*evalf(Pi)*(m-1)/(Nt-1)): end do:

for j from 2 to No-1 do V[j]:=S2[1,j]/dt+1+sin(6*evalf(Pi)*(m-1)/(Nt-1)): end do:

45

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

V[(i-1)*No+j]:=S2[i,j]/dt+1+sin(6*evalf(Pi)*(m-1)/(Nt-1)): end do:

end do:

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

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

S2[i,j]:=Surf[(i-1)*No+j]: end do:

end do:

H_Ice_Summit1[m]:=S1[1,1]:

H_Ice_Summit2[m]:=S2[1,1]:

if m=N1 then for i from 1 to No do S1t1[i]:=S1[i,1]: S2t1[i]:=S2[i,1]: end do: end if:

if m=N2 then for i from 1 to No do S1t2[i]:=S1[i,1]: S2t2[i]:=S2[i,1]: end do: end if:

if m=N3 then for i from 1 to No do S1t3[i]:=S1[i,1]: S2t3[i]:=S2[i,1]: end do: end if:

if m=N4 then for i from 1 to No do S1t4[i]:=S1[i,1]: S2t4[i]:=S2[i,1]: end do: end if:

if m=N5 then for i from 1 to No do S1t5[i]:=S1[i,1]: S2t5[i]:=S2[i,1]: end do: end if:

end do:

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

>plots[sparsematrixplot](M,symbol=circle,symbolsize=5,title="Nonzero matrix elements locations");

46

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

> L_x:=[seq((i-1)*dl,i=1..No)]:

> L_y:=[seq((j-1)*dl,j=1..No)]:

>

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

>L_S:=[seq(S1t1[i],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr1:=plot(S_sp,x=0..1,color=red,thickness=2,linestyle=1,legend="Nt= 100"):

>L_S:=[seq(S1t2[i],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr2:=plot(S_sp,x=0..1,color=red,thickness=2,linestyle=2,legend="Nt= 200"):

>L_S:=[seq(S1t3[i],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr3:=plot(S_sp,x=0..1,color=red,thickness=2,linestyle=3,legend="Nt= 300"):

47

>L_S:=[seq(S1t4[i],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr4:=plot(S_sp,x=0..1,color=red,thickness=2,linestyle=4,legend="Nt= 500"):

>L_S:=[seq(S1t5[i],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr5:=plot(S_sp,x=0..1,color=red,thickness=2,linestyle=5,legend="Nt= 800"):

>L_S:=[seq(S1[i,round(No/2)],i=1..No)]: S_sp:=spline(L_x,L_S,x,1): gr6:=plot(S_sp,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, m"],labeldirections=[HORIZONTAL,VERTICAL]);

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

>L_S:=[seq(S1[i,1],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=0 "):

48

>L_S:=[seq(S1[1,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= 0"):

>plots[display]([gr1,gr2],axes=boxed,title="Ice Cap Surface Elevation at the flowlines",labels=["x,y *200, km","Elevation, nondimensional value"],labeldirections=[HORIZONTAL,VERTICAL]);

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

>

> Accum:=1+sin(6*evalf(Pi)*t/(18*To)):

>

gr1:=plot(Accum,t=0..18*To,axes=boxed,color=red,thickness=3,linest yle=4,legend="accumulation"):

>

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

> L_S:=[seq(H_Ice_Summit1[m],m=0..Nt)]: > S_sp:=spline(L_t,L_S,t,1):

>

gr2:=plot(S_sp,t=0..Nt*dt*To,axes=boxed,color=black,thickness=3,li

49

nestyle=1,legend="Ice Thickness at the Summit, steady accumulation"):

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

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

gr3:=plot(S_sp,t=0..Nt*dt*To,axes=boxed,color=blue,thickness=3,lin estyle=3,legend="Ice Thickness at the Summit, periodic accumulation"):

>plots[display]([gr1,gr2,gr3],axes=boxed,labels=["Time, years","Non-dimensional values"],labeldirections=[HORIZONTAL,VERTICAL],labeldirectio ns=[HORIZONTAL,VERTICAL]);

>
> #	-----------------------------------------------------------------------
> #------------------	The Time of Ice Accumulation ------------

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

>T_accumulation:=To*1.8;

T_accumulation := 6042.040164

> #	------------------------------------------------------------------------
> #-------------- --------------	Ice Cap Thickness at the Summit

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

>Thickness:=S1[1,1]*Zo;

Thickness := 1254.965269

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

50