Коновалов Учебно-методическое пособие по курсу 2007
.pdfsi, 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