Курсовая работа: 2 задание
Описание
Характеристики курсовой работы
Список файлов
- 2 задание
- zad2.f90 8,03 Kb
- zad2.txt 8,01 Kb
Program HDfc
use dflib
integer, parameter :: Nd=500
real(8) r(2,0:Nd),p(0:Nd),Ei(2,0:Nd),x(2,0:Nd+1) ,v(2,0:Nd+1),m(0:Nd),dm,dt,t,r1,r2,e1,e2 ,v1,v2,xmax,m1,m2,rt,gamma1,gamma2,w(0:N d),wg
real(8) pa(0:Nd),cs,cs1,cs2,Dwl1,Dwl2,Dwr1,Dwr2, vc,a1,a2,pi0,Ra1,Ra2,xc,xl1,xl2,xr1,xr2, ra(0:Nd),va(0:Nd),x0,xh,po,em
integer(4) color,paint
integer i,wave1,wave2,im
type(wxycoord) s
call init
call solve
contains
subroutine init
gamma1=7.0/5.0; gamma2=1.4
xmax=5
T1=1500; cv=1.5*0.04155
dt=0.000001
r1=(10**7/cv/T1/(gamma1-1))**0.5; r2=1000;
v1=0; v2=0
e1=r1*cv*T1; e2=10e5/r2/(gamma2-1) !energy
wave1=2; wave2=1
m1=xmax/2.0*r1; m2=m1+xmax/2.0*r2
dm=m2/Nd
em=10**4 !module de Jung
do i=0,Nd
m=dm*(i+0.5)
if (i*dm.le.m1) then
r(1,i)=r1
Ei(1,i)=e1
v(1,i)=v1
im=i
p(i)=pre(r(1,i),Ei(1,i))
else
r(1,i)=r2
Ei(1,i)=e2
v(1,i)=v2
p(i)=pre(r(1,i),Ei(1,i))
end if
x(1,i+1)=x(1,i)+dm/r(1,i)
end do
v(1,Nd+1)=v2
wg=(gamma2+1)/(8*3.14**2)*16*dm**2 !искусственная вязкость
paint=SetWindow(.true.,dble(-0.01),dble( -0.07),dble(xmax+0.05),dble(maxval(r(1,: ))+1))
!call oxyc(dble(0),dble(0),dm*Nd,dble(maxval(r (1,:))+1),dble(0.1),dble(0.1),dble(1),db le(1))
!color=RGBToInteger(0,255,255)
!call res(p,color)
!write(*,*) p(1)-p(Nd)
!call sleepqq(10000)
pi0=(p(0)+p(Nd))/2.0
call iter
end subroutine init
subroutine iter
real(8) :: eps=0.1**6, delta=1
real(8) alpha1,phi1,z1,pi01
do while (delta.gt.eps)
po=0
a1=a(r(1,0),p(0),pi0,wave1,gamma1)
po=em/gamma2
a2=a(r(1,Nd),p(Nd),pi0,wave2,gamma2)
phi1=phi(p(0),p(Nd),v(1,0),v(1,Nd),a1,a2 )
!z1=z(pi0,p(0),p(Nd))
!alpha1=alpha(z1,gamma)
!pi01=(alpha1*pi0+phi1)/(1+alpha1)
pi01=phi1
delta=abs(pi01-pi0)
pi0=pi01
end do
vc=(a1*v(1,0)+a2*v(1,Nd)+p(0)-p(Nd))/(a1 +a2)
po=0
if (wave1.eq.1) then
Ra1=r(1,0)*((gamma1+1)*(pi0+po)+(gamma1- 1)*(p(0))+po)/((gamma1-1)*(pi0+po)+(gamm a1+1)*(p(0)+po))
Dwl1=v(1,0)-a1/r(1,0)
else
cs1=(gamma1*(p(0)+po)/r(1,0))**0.5+0.5*( gamma1-1)*(v(1,0)-vc)
Ra1=gamma1*(pi0+po)/(cs1)**2
Dwl1=v(1,0)-(gamma1*(p(0)+po)/r(1,0))**0 .5
Dwl2=vc-cs1
end if
po=em/gamma2
if (wave2.eq.1) then
Ra2=r(1,Nd)*((gamma2+1)*(pi0+po)+(gamma2 -1)*(p(Nd)+po))/((gamma2-1)*(pi0+po)+(ga mma2+1)*(p(Nd)+po))
Dwr1=v(1,Nd)+a2/r(1,Nd)
else
cs2=(gamma2*(p(Nd)+po)/r(1,Nd))**0.5-0.5 *(gamma2-1)*(v(1,Nd)-vc)
Ra2=gamma2*pi0/(cs2)**2
Dwr1=v(1,Nd)+(gamma2*(p(Nd)+po)/r(1,Nd)) **0.5
Dwr2=vc+cs2
end if
end subroutine iter
subroutine solve
do t=dt,15000*dt,dt
do i=0,Nd-1
v(2,i+1)=v(1,i+1)-dt/dm*(p(i+1)-p(i))
x(2,i+1)=x(1,i+1)+0.5*dt*(v(2,i+1)+v(1,i +1))
rt=1.0/r(1,i)+0.5*dt/dm*(v(2,i+1)-v(2,i) +v(1,i+1)-v(1,i))
r(2,i)=1.0/rt
Ei(2,i)=Ei(1,i)-p(i)*0.5*dt/dm*(v(2,i+1) -v(2,i)+v(1,i+1)-v(1,i))
end do
r(2,0)=r(2,1); Ei(2,0)=Ei(2,1); r(2,Nd)=r(2,Nd-1); Ei(2,Nd)=Ei(2,Nd-1); v(2,0)=v(2,1); v(2,Nd+1)=v(2,Nd)
x(2,0)=x(1,0)+0.5*dt*(v(2,0)+v(1,0)); x(2,Nd+1)=x(1,Nd+1)+0.5*dt*(v(2,Nd+1)+v( 1,Nd+1))
! искусственая вязкость
do i=0,Nd
if (v(2,i+1)-v(2,i).lt.0) then
w(i)=r(2,i)*wg/dm*(v(2,i+1)-v(2,i))**2
else
w(i)=0
end if
if (i*dm.le.m1) then
p(i)=pre(r(2,i),Ei(2,i))+w(i)
else
p(i)=pre1(r(2,i),Ei(2,i))+w(i)
end if
end do
if (mod(int(t/dt),20).eq.0) then
!paint=SetWindow(.true.,dble(-0.01),dble (-0.07),xmax+0.05,dble(maxval(p(:))+100) )
paint=SetBKColor(0)
call ClearScreen(0)
x0=0.5*xmax
xc=0.5*xmax+vc*t
xl1=0.5*xmax+Dwl1*t
xl2=0.5*xmax+Dwl2*t
xr1=0.5*xmax+Dwr1*t
xr2=0.5*xmax+Dwr2*t
do i=0,Nd
xh=x(2,i)
if (wave1.eq.1) then
if (x(2,i).lt.xl1) then
ra(i)=r(1,0)
pa(i)=p(0)
va(i)=v(1,0)
else if ((x(2,i).ge.xl1).and.(x(2,i).lt.xc)) then
ra(i)=Ra1
pa(i)=pi0
va(i)=vc
end if
else
if (x(2,i).lt.xl1) then
ra(i)=r(1,0)
pa(i)=p(0)
va(i)=v(1,0)
else if ((x(2,i).ge.xl1).and.(x(2,i).lt.xl2)) then
va(i)=((gamma1-1)*v(1,0)+2*((gamma1*p(0) /r(1,0))**0.5-(xh-x0)/t))/(gamma1+1)
cs=(2*(gamma1*p(0)/r(1,0))**0.5+(gamma1- 1)*(v(1,0)-(xh-x0)/t))/(gamma1+1)
ra(i)=((cs**2*r(1,0)**gamma1)/(gamma1*p( 0)))**(1.0/(gamma1-1))
pa(i)=cs**2*ra(i)/gamma1
else if ((x(2,i).ge.xl2).and.(x(2,i).lt.xc)) then
ra(i)=Ra1
pa(i)=pi0
va(i)=vc
end if
end if
if (wave2.eq.1) then
if ((x(2,i).ge.xc).and.(x(2,i).lt.xr1)) then
ra(i)=Ra2
pa(i)=pi0
va=vc
else if (x(2,i).ge.xr1) then
ra(i)=r(1,Nd)
pa(i)=p(Nd)
va(i)=v(1,Nd)
end if
else
if ((x(2,i).ge.xc).and.(x(2,i).lt.xr2)) then
ra(i)=Ra2
pa(i)=pi0
va=vc
else if ((x(2,i).ge.xr2).and.(x(2,i).lt.xr1)) then
cs=(2.0*(gamma2*p(Nd)/r(1,Nd))**0.5+((xh -x0)/t-v(1,Nd))*(gamma2-1))/(gamma2+1)
va(i)=(xh-x0)/t-cs
ra(i)=((cs**2*r(1,Nd)**gamma2)/(gamma2*p (Nd)))**(1.0/(gamma2-1))
pa(i)=cs**2*ra(i)/gamma2
else if (x(2,i).ge.xr1) then
ra(i)=r(1,Nd)
pa(i)=p(Nd)
va(i)=v(1,Nd)
end if
end if
end do
paint=SetWindow(.true.,dble(-0.01),dble( -1),xmax+0.05,dble(10e6))
color=RGBToInteger(255,255,0)
call res(x(2,:),p-w,color)
color=RGBToInteger(0,255,255)
call res(x(2,:),r(2,:)*10**3,color)
color=RGBToInteger(255,0,0)
call res(x(2,:),ra(:)*10**3,color)
call res(x(2,:),pa(:),color)
call sleepqq(30)
end if
r(1,:)=r(2,:); Ei(1,:)=Ei(2,:); v(1,:)=v(2,:); x(1,:)=x(2,:)
end do
end subroutine solve
subroutine oxyc(x0,y0,x1,y1,xd,yd,xm,ym)
use dflib
real(8) x0,x1,y0,y1,xd,yd,xm,ym,x,y
integer(4) color
integer(2) paint,font
character(4) l
type(wxycoord) s
paint=SetBKColor(0)
call ClearScreen(0)
font=InitializeFonts()
font=SetFont('t''Arial Cyr''h14w7b')
do x=x0,x1+xd/2,xd
write(l,"(f4.1)") x
if (x/xm-int(x/xm).lt.0.001) then
color=RGBToInteger(0,255,0)
else
color=RGBToInteger(0,100,0)
end if
paint=SetColorRGB(color)
call MoveTo_w(x,y0,s)
paint=LineTo_w(x,y1)
call MoveTo_w(x-0.008,y0-0.02,s)
call outGText(l)
end do
do y=y0,y1+yd/2,yd
write(l,"(f4.1)") y
if (y/ym-int(y/ym).lt.0.001) then
color=RGBToInteger(0,255,0)
else
color=RGBToInteger(0,100,0)
end if
paint=SetColorRGB(color)
call MoveTo_w(x0,y,s)
paint=LineTo_w(x1,y)
call MoveTo_w(x0-0.018,y+0.015,s)
call outGText(l)
end do
end subroutine oxyc
subroutine res(x,y,color)
Integer(4) color
real(8) x(0:Nd),y(0:Nd)
paint=SetColorRGB(color)
Старожил
