Хөндийн урсгалыг ЛБА-аар кодлох
Хөндийн урсгал буюу квадрат саванд орших шингэн квадратын дээд хүрээний
шилжилтээр хөдөлгөөнд орж саван дотроо эргэлдэн урсах бодлого нь тооцон бодох
шингэний динамик дахь жишиг (benchmark) бодлогуудын нэг юм. Латтис больцманы аргаар энэхүү
хөндийн урсгалыг хэрхэн кодлох талаар авч үзье. Квадрат хөндий нь 0.2м талтай
бөгөөд 15 градусын тосоор дүүргэгдсэн. Уг тосны кинемтик зунгааралт нь alpha=1.2*10-3 м2/с байх ба квадратын дээд
хүрээ ulid=6м/с хурдтайгаар хөдлөнө. Изотерм шингэний урсгалд Латтис больцманы аргыг
амжилттай ашиглахын тулд Рейнольдсын тоо ба геометр нөлөөллийг сайтар
тохируулах шаардлагатай. Бодит Рейнольдсын тоог тодорхойлбол:
Хөндийн энэ бодлогийн хувьд геометр харьцаа нь 1 байна. ЛБА-р тооцоо хийхэд
хэмжээсгүй ямарч хурдыг /сүлжээний хурд гэе/ квадратын дээд талын шилжих
хурдаар авч болох ба сүлжээний зунгааралт нь бодит Рейнольдсын тооны утгыг
хангаж байх ёстой. ЛБА-ын тооцооны шахагдлын нөлөөг бууруулахын тулд сүлжээний
хурдыг 0.1 – ээр, сүлжээний зунгааралтыг 0.01 гэж авъя. Үүнд тохирох геометр
хэмжигдэхүүн нь:
Ингээд
дурын гэхдээ хазаартай /хязгаар/ сонгож авсан хурд ба зунгааралтын коэф-той үед
бодит Рейнольдсын тоог хангахуйц сүлжээний тоо нь х ба у тэнхлэгийн турш 100,
хөндийн бүхэлд 1000 ширхэг байх боллоо. Ингээд кодыг нэг бүрчлэн шинжилж үзье.
parameter (n=100,m=100)
real f(0:8,0:n,0:m)
real feq(0:8,0:n,0:m),rho(0:n,0:m)
real w(0:8), cx(0:8),cy(0:8)
real u(0:n,0:m), v(0:n,0:m)
integer i
Сүлжээний тоо 100х100 байх ба уг сүлжээ нь D2Q9 буюу 9 ширхэг түгэлтийн функцтай байх учир f(0:8,0:n,0:m) гэсэн гурван хэмжээст бүрдэл /array
гэдгийг бүрдэл гэнэ./ хэмжигдэхүүн байна. Мөн үүнтэй адил
хэмжээстэйгээр тэнцвэрт түгэлтийн функц feq(0:8,0:n,0:m)
гэж зарлагдана. Сүлжээний зангилаа бүрд бодох гэж буй шингэний
параметрүүд байх ёстой учир нягт, босоо хурд, хэвтээ хурдуудыг харгалзан rho(0:n,0:m), v(0:n,0:m), u(0:n,0:m) гэж зарласан
байна. Энд байх w(0:8), cx(0:8), cy(0:8) эдгээр
нь жинлэх хүчин зүйл, хэвтээ сүлжээний хурд, босоо сүлжээний хурд гэх мэт
тэнцвэрт түгэлтийн функцийг тооцоход хэрэглэгдэх ЛБА-ын параметрүүд ба эдгээр
нь 0-ээс 8 хүртэл дугаарлагдах нэг хэмжээст бүрдэл байна.
uo=0.10 ! lattice lid velocity
rhoo=5.00
dx=1.0
dy=dx
dt=1.0
alpha=0.01 ! lattice viscosity
Re=uo*m/alpha ! Lattice Reynolds
omega=1.0/(3.*alpha+0.5)
mstep=40000
Ингээд квадрат хөндийн дээд талын шилжих хурдыг uo=0.10,
шингэний нягтыг rhoo=5.00, сүлжээ хоорондох
зайн алхамыг dy=dx=1.0, сүлжээний
зунгааралтыг alpha=0.01, тус тус өгч
Рейнольдсын тоог олж, тайвшралтын давтамжийг omega=1.0/(3.*alpha+0.5)
гэж тодорхойлж байна.
w(0)=4./9.
do i=1,4
w(i)=1./9.
end do
do i=5,8
w(i)=1./36.
end do
cx(0)=0
cx(1)=1
cx(2)=0
cx(3)=-1
cx(4)=0
cx(5)=1
cx(6)=-1
cx(7)=-1
cx(8)=1
cy(0)=0
cy(1)=0
cy(2)=1
cy(3)=0
cy(4)=-1
cy(5)=1
cy(6)=1
cy(7)=-1
cy(8)=-1
Энд жинлэх хүчин зүйл болон сүлжээний нэгж хурднуудын утгыг өгч байна. Дараах
зурагтай тус бүрийг харьцуулж харна уу.
do j=0,m
do i=0,n
rho(i,j)=rhoo
u(i,j)=0.0
v(i,j)=0.0
end do
end do
do i=1, n-1
u(i,m)=uo
v(i,m)=0.0
end do
Дээрх хэсэгт сүлжээ бүр дээр шингэний нягт ба анхны хурдуудыг өгч байна.
Харин хөндийн дээд хүрээний сүлжээнүүд нь uo гэсэн
сүлжээний хурдтай баруун тийш шилжиж байх учир n-1 байрлал
дахь хэвтээ хурдыг u(i,m)=uo гэж өгнө.
1 do kk=1,mstep
call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
call streaming(f,n,m)
call sfbound(f,n,m,uo)
call rhouv(f,rho,u,v,cx,cy,n,m)
print *,
u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
write(8,*) kk,u(n/2,m/2),v(n/2,m/2)
END DO
Энэ нь бодлогийн үндсэн гогцоо бөгөөд гол бодолт энэ дотор явагдана. Тус
тус дуудагдаж байгаа дэд програмуудыг тайлбарлавал: call
collesion нь түгэлтийн функцын мөргөлдөөнийг тодорхойлно, call streaming нь мөргөлдөөний дараах шинэчлэгдсэн
түгэлтийн функцуудыг дараачийн сүлжээрүү сүлжээний нэгж хурдын утгаар
шилжүүлнэ, call sfbound нь урсгалын дараах
хязгаар дээрх сүлжээний хязгаарын нөхцлийг шийдвэрлэнэ, call rhouv гэдэг нь нягт, босоо ба хэвтээ хурдыг түгэлтийн
функцын тусламжтайгаар бодно. Ингээд тус бүр дэд кодуудыг авч үзье.
Мөргөлдөөний дэд код
DO i=0,n
DO j=0,m
t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
DO k=0,8
t2=u(i,j)*cx(k)+v(i,j)*cy(k)
feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
END DO
END DO
END DO
Анхны хугацааны алхамд хөндий дотор нягт ба хурнууд өгөгдсөн тайван байдалд
байх ба энэ үеийн түгэлтийн функцууд нь 0 гэсэн утгатай байна. Харин тэнцвэрт
түгэлтийн функцыг дараах байдлаар тодорхойлно.
Дээрх кодонд бичигдсэн t1=u(i,j)*u(i,j)+v(i,j)*v(i,j) нь u2 бөгөөд t2=u(i,j)*cx(k)+v(i,j)*cy(k)
нь c*u
юм. Тэнцвэрт түгэлтийн функцыг feq(k,i,j)
гэж тодорхойлсоны дараа мөргөлдөөний үе дахь түгэлтийн функцыг f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j) гэж
бодно. Энд омега нь дээр бодогдсон тайвшралтын давтамж юм.
Урсгалын дэд код
DO j=0,m
DO
i=n,1,-1 !RIGHT TO LEFT
f(1,i,j)=f(1,i-1,j)
END DO
DO i=0,n-1 !LEFT TO RIGHT
f(3,i,j)=f(3,i+1,j)
END DO
END DO
DO
j=m,1,-1 !TOP TO BOTTOM
DO
i=0,n
f(2,i,j)=f(2,i,j-1)
END DO
DO
i=n,1,-1
f(5,i,j)=f(5,i-1,j-1)
END DO
DO
i=0,n-1
f(6,i,j)=f(6,i+1,j-1)
END DO
END DO
DO j=0,m-1 !BOTTOM TO TOP
DO i=0,n
f(4,i,j)=f(4,i,j+1)
END DO
DO i=0,n-1
f(7,i,j)=f(7,i+1,j+1)
END DO
DO i=n,1,-1
f(8,i,j)=f(8,i-1,j+1)
END DO
END DO
Дээрхийг тайлбарлавал баруунаас зүүн тийш урсах түгэлтийн функц нь зөвхөн
баруун тийш харсан f1 /дээрх
схем зургийг харах/ байна. Иймд f(1,i-1,j) гэдэг мөргөлдөөний үед тодорхойлогдсон функцын
утга сүлжээгээ орхин дараачийн сүлжээний 1-р шилбэн дээр буюу f(1,i,j) –т ирэх
ёстой. Зүүнээс баруун мөн ялгаагүй. Харин дээшээ харж байрласан f2,
f5, f6-ын хувьд дээрээс доошоо нүүх
ёстой. Гэхдээ f2 нь
зэргэлдээ дээд сүлжээрүү, f5 нь даамны нүүдлээр өрөөлдөж чиглэлийн дагуу, f6 нь мөн өрөөлдөж нөгөө чиглэлийн дагуу дээд мөрний сүлжээрүү
хөдлөнө. Үүнтэй адил доош харж байрласан түгэлтийн функцууд доошоо нүүх ёстой.
Програмын утгын шилжүүлэг бидний бодож байгаагаас эсрэгээр байгааг сайн
анхаарах хэрэгтэй. f(2,i,j)=f(2,i,j-1) гэдэг нь доод сүлжээнээс дээшээ нүүж
байгааг илтгэнэ. Доорхи зураг нь урсгалын процессыг товч тодорхой харуулна.
Хязгаарын нөхцлийн дэд код
! bounce back on west boundary
f(1,0,j)=f(3,0,j)
f(5,0,j)=f(7,0,j)
f(8,0,j)=f(6,0,j)
! bounce back on east boundary
f(3,n,j)=f(1,n,j)
f(7,n,j)=f(5,n,j)
f(6,n,j)=f(8,n,j)
end do
! bounce back on south boundary
do
i=0,n
f(2,i,0)=f(4,i,0)
f(5,i,0)=f(7,i,0)
f(6,i,0)=f(8,i,0)
end do
! moving lid, north boundary
do
i=1,n-1
rhon=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
f(4,i,m)=f(2,i,m)
f(8,i,m)=f(6,i,m)+rhon*uo/6.0
f(7,i,m)=f(5,i,m)-rhon*uo/6.0
end do
Дөрвөн тал тус бүр дээр буцаж ойх схем бүхий хязгаарын нөхцлийг хэрэглэж
байна. Буцаж ойх схемээс энд хамгийн энгийн схем буюу хязгаар дээр сүлжээ
байрласан үеийн схемийг хэрэглэж байна. Жишээ болгон баруун тал дээрх хязгаарын
нөхцлийг тайлбарлая. Саяны урсгалын дэд програм дээр f5,
f1, f8 гэсэн түгэлтийн функцууд
уг баруун хязгаар дээр утгагүй байгаа ба f7, f3, f6 гэсэн функцууд хязгаараас гадагш чиглээд байрлаж байгаа. Тэгэхээр
уг гадагш чиглэсэн функцуудыг өнчин хоцорсон f5, f1, f8-д харгалзуулан буцааж ойлгосоноор энэхүү хязгаарын
нөхцлийн хэрэг гарч асуудал шийдэгдэж байгаа юм. Учирыг нь ойлгохгүй бол асууж
болно. Бусад хязгаарууд яг адилхан байх ба энд онцлог нэг хязгаар нь квадратын
дээд тал дээрх сүлжээнүүд юм. Энэ нь бодлогийн хөдөлгөгч болох хэвтээ хурдыг
шийдэж байх хязгаар юм. Оролтын хурдтай хязгаарын нөхцлийг хязгаартай нормаль
тэнцвэрийн тэгшитгэл ашиглан шийдэврлэх ба уг аргачлалаар дараах тэгшитгэлүүд
гарч ирнэ.
Эндээс босоо
хурд нь тэг учир түүнийг агуулсан гишүүд үгүй болж кодлогдсон байна.
Хамаарах хэмжигдэхүүнийг үнэлэх дэд код
Энэ
бодлогийн хувьд хамаарах хэмжигдэхүүн нь нягт, босоо ба хэвтээ хурд орох ба
эдгээрийг дараах байдлаар кодолж өгсөн байна.
do j=0,m
do i=0,n
ssum=0.0
do k=0,8
ssum=ssum+f(k,i,j)
end do
rho(i,j)=ssum
end do
end do
do i=1,n
rho(i,m)=2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
rho(i,m)=rho(i,m)+f(0,i,m)+f(1,i,m)+f(3,i,m)
end do
DO
i=1,n
DO j=1,m-1
usum=0.0
vsum=0.0
DO k=0,8
usum=usum+f(k,i,j)*cx(k)
vsum=vsum+f(k,i,j)*cy(k)
END DO
u(i,j)=usum/rho(i,j)
v(i,j)=vsum/rho(i,j)
END DO
END DO
Эхний
хэсэгт ssum=ssum+f(k,i,j) гэдэгээр бүх
сүлжээний /1000ширхэг сүлжээний/ 9 ширхэг тзгэлтийн функцын нийлбэрийг олж
байна. энэ нь тухайн зангилаан дээрх нягттай тэнцүү байна гэдгийг дараагийн rho(i,j)=ssum гэдэг томъёоллоор харуулж байна.
Харин квадрат хөндийн дээд хязгаар дээрх нягтыг олохдоо rho(i,m)= 2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))+f(0,i,m)+f(1,i,m)+f(3,i,m) гэсэн
урт томъёоллыг ашиглах ба үүнийг кодонд задалж бичсэн учир хоёр удаа rho(i,m) гарч байна. Харин хурд нь түгэлтийн
функцуудыг түүний нэгж хурдаар үржсэний нийлбэрийг нягтад харьцуулсантай тэнцүү
буюу момент хадгалагдах илэрхийллээс u(i,j)=usum/rho(i,j)
гэж олдож байна. Босоо хурдын хувьд мөн ялгаагүй юм. Хөндийн урсгалд мөн
урсгалын шугамын функцыг шийдвэрлэж өгөх шаардлагатай ба товчоор урсгалын функц
нь дараах байдлаар илэрхийлэгдэнэ.
do i=0,n
rhoav=0.5*(rho(i-1,0)+rho(i,0))
if(i.ne.0)
strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
do j=1,m
rhom=0.5*(rho(i,j)+rho(i,j-1))
strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
end do
end do
Ингээд
үр дүнг сонирхвол:
Урсгалын
вектор
Урсгалын
хэвтээ хурдны контур
Урсгалын
босоо хурдны контур
Урсгалын
шугамын функын утга буюу урсгалын шугам.
Дээрх үр
дүнгүүд нь ЛБА-д тохируулсан хэмжээсгүй үр дүнгүүд ба 100 нэгж нь 0.2м-ийг
илэрхийлж, хурдны 0.1 хэмжээ нь бодит хурдны 6м/с – ийг илтгэх ёстой. Бодлогийг
анх яаж масштбалсан тэр журмаараа хөрвүүлж үр дүнг авна.
Нэг хэмжээст ЛБА-ын кодыг эндээс
Гурван хэмжээст ЛБА-ын кодыг эндээс
Ингээд бүрэн кодыг хавсаргавал
1: ! computer code for lid-driven cavity
2: parameter (n=100,m=100)
3: real f(0:8,0:n,0:m)
4: real feq(0:8,0:n,0:m),rho(0:n,0:m)
5: real w(0:8), cx(0:8),cy(0:8)
6: real u(0:n,0:m), v(0:n,0:m)
7: integer i
8: open(2,file='uvfield.dat')
9: open(3,file='uvely.dat')
10: open(4,file='vvelx.dat') ! main result
11: open(8,file='timeu.dat') ! time with u and v
12: ! Initialazation
13: uo=0.10 ! lattice lid velocity
14: rhoo=5.00
15: dx=1.0
16: dy=dx
17: dt=1.0
18: alpha=0.01 ! lattice viscosity
19: Re=uo*m/alpha ! Lattice Reynolds
20: print *, 'Re=', Re
21: omega=1.0/(3.*alpha+0.5)
22: mstep=40000
23: w(0)=4./9.
24: do i=1,4
25: w(i)=1./9.
26: end do
27: do i=5,8
28: w(i)=1./36.
29: end do
30: cx(0)=0
31: cx(1)=1
32: cx(2)=0
33: cx(3)=-1
34: cx(4)=0
35: cx(5)=1
36: cx(6)=-1
37: cx(7)=-1
38: cx(8)=1
39: cy(0)=0
40: cy(1)=0
41: cy(2)=1
42: cy(3)=0
43: cy(4)=-1
44: cy(5)=1
45: cy(6)=1
46: cy(7)=-1
47: cy(8)=-1
48: do j=0,m
49: do i=0,n
50: rho(i,j)=rhoo
51: u(i,j)=0.0
52: v(i,j)=0.0
53: end do
54: end do
55: do i=1,n-1
56: u(i,m)=uo
57: v(i,m)=0.0
58: end do
59: ! main loop
60: 1 do kk=1,mstep
61: call collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
62: call streaming(f,n,m)
63: ! *********************************
64: call sfbound(f,n,m,uo)
65: call rhouv(f,rho,u,v,cx,cy,n,m)
66: print *, u(0,m/2),v(0,m/2),rho(0,m/2),u(n,m/2),v(n,m/2),rho(n,m/2)
67: write(8,*) kk,u(n/2,m/2),v(n/2,m/2)
68: END DO
69: ! end of the main loop
70: call result(u,v,rho,uo,n,m)
71: stop
72: end
73: ! end of the main program
74: ! ********************************************************!
75: subroutine collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m)
76: real f(0:8,0:n,0:m)
77: real feq(0:8,0:n,0:m),rho(0:n,0:m)
78: real w(0:8), cx(0:8),cy(0:8)
79: real u(0:n,0:m), v(0:n,0:m)
80: DO i=0,n
81: DO j=0,m
82: t1=u(i,j)*u(i,j)+v(i,j)*v(i,j)
83: DO k=0,8
84: t2=u(i,j)*cx(k)+v(i,j)*cy(k)
85: feq(k,i,j)=rho(i,j)*w(k)*(1.0+3.0*t2+4.50*t2*t2-1.50*t1)
86: f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j)
87: END DO
88: END DO
89: END DO
90: return
91: end
92: ! ********************************************************!
93: subroutine streaming(f,n,m)
94: real f(0:8,0:n,0:m)
95: ! streaming
96: DO j=0,m
97: DO i=n,1,-1 !RIGHT TO LEFT
98: f(1,i,j)=f(1,i-1,j)
99: END DO
100: DO i=0,n-1 !LEFT TO RIGHT
101: f(3,i,j)=f(3,i+1,j)
102: END DO
103: END DO
104: DO j=m,1,-1 !TOP TO BOTTOM
105: DO i=0,n
106: f(2,i,j)=f(2,i,j-1)
107: END DO
108: DO i=n,1,-1
109: f(5,i,j)=f(5,i-1,j-1)
110: END DO
111: DO i=0,n-1
112: f(6,i,j)=f(6,i+1,j-1)
113: END DO
114: END DO
115: DO j=0,m-1 !BOTTOM TO TOP
116: DO i=0,n
117: f(4,i,j)=f(4,i,j+1)
118: END DO
119: DO i=0,n-1
120: f(7,i,j)=f(7,i+1,j+1)
121: END DO
122: DO i=n,1,-1
123: f(8,i,j)=f(8,i-1,j+1)
124: END DO
125: END DO
126: return
127: end
128: ! ********************************************************!
129: subroutine sfbound(f,n,m,uo)
130: real f(0:8,0:n,0:m)
131: do j=0,m
132: ! bounce back on west boundary
133: f(1,0,j)=f(3,0,j)
134: f(5,0,j)=f(7,0,j)
135: f(8,0,j)=f(6,0,j)
136: ! bounce back on east boundary
137: f(3,n,j)=f(1,n,j)
138: f(7,n,j)=f(5,n,j)
139: f(6,n,j)=f(8,n,j)
140: end do
141: ! bounce back on south boundary
142: do i=0,n
143: f(2,i,0)=f(4,i,0)
144: f(5,i,0)=f(7,i,0)
145: f(6,i,0)=f(8,i,0)
146: end do
147: ! moving lid, north boundary
148: do i=1,n-1
149: rhon=f(0,i,m)+f(1,i,m)+f(3,i,m)+2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
150: f(4,i,m)=f(2,i,m)
151: f(8,i,m)=f(6,i,m)+rhon*uo/6.0
152: f(7,i,m)=f(5,i,m)-rhon*uo/6.0
153: end do
154: return
155: end
156: ! ********************************************************!
157: subroutine rhouv(f,rho,u,v,cx,cy,n,m)
158: real f(0:8,0:n,0:m),rho(0:n,0:m),u(0:n,0:m),v(0:n,0:m)
159: real cx(0:8),cy(0:8)
160: do j=0,m
161: do i=0,n
162: ssum=0.0
163: do k=0,8
164: ssum=ssum+f(k,i,j)
165: end do
166: rho(i,j)=ssum
167: end do
168: end do
169: do i=1,n
170: rho(i,m)=2.*(f(2,i,m)+f(6,i,m)+f(5,i,m))
171: rho(i,m)=rho(i,m)+f(0,i,m)+f(1,i,m)+f(3,i,m)
172: end do
173: DO i=1,n
174: DO j=1,m-1
175: usum=0.0
176: vsum=0.0
177: DO k=0,8
178: usum=usum+f(k,i,j)*cx(k)
179: vsum=vsum+f(k,i,j)*cy(k)
180: END DO
181: u(i,j)=usum/rho(i,j)
182: v(i,j)=vsum/rho(i,j)
183: END DO
184: END DO
185: return
186: end
187: ! ********************************************************!
188: subroutine result(u,v,rho,uo,n,m)
189: real u(0:n,0:m),v(0:n,0:m)
190: real rho(0:n,0:m),strf(0:n,0:m)
191: ! open(5, file='streamf.dat')
192: ! stream function calculations
193: strf(0,0)=0.
194: do i=0,n
195: rhoav=0.5*(rho(i-1,0)+rho(i,0))
196: if(i.ne.0) strf(i,0)=strf(i-1,0)-rhoav*0.5*(v(i-1,0)+v(i,0))
197: do j=1,m
198: rhom=0.5*(rho(i,j)+rho(i,j-1))
199: strf(i,j)=strf(i,j-1)+rhom*0.5*(u(i,j-1)+u(i,j))
200: end do
201: end do
202: ! ———————————–
203: write(2,*)' VARIABLES =X, Y, U, V, S'
204: write(2,*)'ZONE ','I=',n+1,'J=',m+1,',','F=BLOCK'
205: do j=0,m
206: write(2,*)(i,i=0,n)
207: end do
208: do j=0,m
209: write(2,*)(j,i=0,n)
210: end do
211: do j=0,m
212: write(2,*)(u(i,j),i=0,n)
213: end do
214: do j=0,m
215: write(2,*)(v(i,j),i=0,n)
216: end do
217: do j=0,m
218: write(2,*)(strf(i,j),i=0,n)
219: end do
220: do j=0,m
221: write(3,*)j/float(m),u(n/2,j)/uo,u(n/4,j)/uo,u(3*n/4,j)/uo
222: end do
223: do i=0,n
224: write(4,*) i/float(n),v(i,m/2)/uo
225: end do
226: return
227: end
228: !============end of the program
Төгсөв.
No comments:
Post a Comment