2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




 
 Азбучный вопрос по методу крупных частиц
Сообщение22.02.2013, 23:00 
Аватара пользователя
Здравствуйте.

Я пытаюсь запрограммировать одномерный метод крупных частиц первого порядка без модификаций в точности по книжке Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике.
Решил протестировать на задаче Сода (то есть задача о распаде разрыва с начальными условиями слева $p = 1$, $\rho = 1$, $u = 0$ и справа $p = 0.1$, $\rho = 0.125$, $u=0$). Получается не то, чтобы очень хорошо, но терпимо

Изображение

Хотя провал давления в центре меня несколько смущает. Явно от аналитического решения отличается. Можно ли это как-то улучшить?

И более важно то, что когда ударная волна доходит до правой границы расчетной области, схему разносит. Я там ставлю условие отражения $v(n, Nx+1) = - v(n, Nx)$.

Как быть? Буду очень благодарен за любую помощь.

Прилагаю свой код на фортране.
код: [ скачать ] [ спрятать ]
Используется синтаксис Fortran
program KRUCHA_1D_Sode

    implicit none

    ! Variables
    integer, parameter:: Nt=2500, Nx=300
    real Lx, dx, dt, gamma, Ku, TIME, pr, pl, ur, ul, flowR, flowL, dMR, dML, Xij, vc, vmax
    real u(0:Nt,0:Nx+1), rho(0:Nt,0:Nx+1), p(0:Nt,0:Nx+1), E(0:Nt,0:Nx+1)
    real ui(0:Nx+1), rhoi(0:Nx+1), pi(0:Nx+1), Ei(0:Nx+1)
    integer n, i, DL, DR

    ! Body of KRUCHA_1D_Sode
    Lx=2.0
    gamma=1.4
   
    ! Initial conditions
    do i=0, Nx+1
        if(i.le.(Nx/2)) then
            u(0,i)=0.0
            rho(0,i)=1.0
            p(0,i) = 1.0
            E(0,i) = p(0,i)/(rho(0,i)*(gamma-1))
            else
            u(0,i)=0.0
            rho(0,i)=0.125
            p(0,i) = 0.1
            E(0,i) = p(0,i)/(rho(0,i)*(gamma-1))
         endif
    end do
   
    ! Boundary conditions at initial step
    u(0,0) = - u(0,1); u(0,Nx+1) = -u(0,Nx)
   
    ! Time & scheme parameters
    TIME = 0.0
    Ku = 0.1
    dt = Ku*dx/(sqrt(gamma*p(0,1)/rho(0,1)))
    dx=Lx/Nx
   
    ! Main cycle
    do n=0, Nt-1
        TIME = TIME + dt
        vmax = 0.0
        ! I. Euler step
        do i=1,Nx
            pr = 0.5*(p(n,i+1)+p(n,i)); pl = 0.5*(p(n,i-1)+p(n,i))
            ur = 0.5*(u(n,i+1)+u(n,i)); ul = 0.5*(u(n,i-1)+u(n,i))
            ui(i) = u(n,i) - (pr - pl)*dt/(dx*rho(n,i))
            Ei(i) = E(n,i) - (pr*ur - pl*ul)*dt/(dx*rho(n,i))
        end do
       
        ! Boundary conditions at Euler step
        ui(0) = -ui(1); Ei(0) = Ei(1)
        ui(Nx+1) = -ui(Nx); Ei(Nx+1) = Ei(Nx)
       
        ! Lagrange & final step cycle
        do i=1,Nx
            ! II. Lagrange step
            flowR = 0.5*(ui(i) + ui(i+1))
            if (flowR.ge.0) then
                dMR = rho(n,i)*flowR*dt; DR = 0
                else
                dMR = rho(n,i+1)*flowR*dt; DR = 1
            end if
           
            flowL = 0.5*(ui(i) + ui(i-1))
            if (flowL.ge.0) then
                dML = rho(n,i-1)*flowL*dt; DL = 1
                else
                dML = rho(n,i)*flowL*dt; DL = 0
            end if
            ! III. Final step
            rho(n + 1,i) = (rho(n,i)*dx + (dML - dMR))/dx
            Xij = rho(n+1,i)*dx - (1-DR)*dMR - (1-DL)*dML
            u(n+1,i) = (DL*dML*ui(i-1) - DR*dMR*ui(i+1) + ui(i)*Xij)/(rho(n+1,i)*dx)
            E(n+1,i) = (DL*dML*Ei(i-1) - DR*dMR*Ei(i+1) + Ei(i)*Xij)/(rho(n+1,i)*dx)
            p(n+1,i) = (gamma - 1)*rho(n+1,i)*(E(n+1,i) - 0.5*u(n+1,i)*u(n+1,i))
           
            ! Maximum velocity
            vc = Abs(u(n+1,i))+Sqrt(gamma*p(n+1,i)/rho(n+1,i))
            if (vc.gt.vmax) then ; vmax = vc; endif
        end do
       
        ! Boundary conditions at final step
        rho(n+1,0) = rho(n+1,1); u(n+1,0) = -u(n+1,1); p(n+1,0) = p(n+1,1); E(n+1,0) = E(n+1,1)
        rho(n+1,Nx+1) = rho(n+1,Nx); u(n+1,Nx+1) = -u(n+1,Nx); p(n+1,Nx+1) = p(n+1,Nx); E(n+1,Nx+1) = E(n+1,Nx);
       
        ! Time step
        dt = Ku*dx/vmax
    end do
   
    end program

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение26.02.2013, 21:32 
Аватара пользователя
Roy Rogers в сообщении #687151 писал(а):
Можно ли это как-то улучшить?

Можно, конечно. У вас грубо говоря три дороги: TVD, ENO, SPH.

P.S. А сам доисторический метод Белоцерковского - незамедлительно в печку, разумеется.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение16.03.2013, 15:10 
Аватара пользователя
Цитата:
метод Белоцерковского - незамедлительно в печку...

Утундрий, то что в методе Харлоу-Белоцерковского использовано, обгорело в сотнях возвращаемых тел на Махе 20 без видимых расхождений с экспериментом не сопряжено с Вашими суждениями.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение17.03.2013, 01:25 
Аватара пользователя
Боевая колесница ездить, конечно, могёт, но приходит время и вместе с ним танки.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение23.03.2013, 13:05 
Аватара пользователя
Ваши "TVD, ENO, SPH" - убогая реинкарнация по времени, без повышения порядка аппроксимации на разрывах по пространственным координатам( не говоря о излучении и ряду других важных физических процессов), с заметным увеличением вычислительной вязкости и дикой недокументированной дисперсии. К тому же они сейчас уже бордаты.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение23.03.2013, 18:05 
Аватара пользователя
Zai
Думайте как хотите.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение08.05.2013, 20:44 
я тоже осваиваю метод крупных частиц
Скажите, а в каком компиляторе fortran вы работаете?, потому что у меня при вводе вашего кода кучу ошибок выкидывает

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение08.05.2013, 21:23 
Аватара пользователя
kor_viktor
У него стандарт 95 или выше.

 
 
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение16.05.2016, 11:20 
Здравствуйте, на основании чего Вами был выбран шаг?
Код:
dt = Ku*dx/(sqrt(gamma*p(0,1)/rho(0,1)))

и что Вы здесь вычисляете:
Код:
! Maximum velocity
            vc = Abs(u(n+1,i))+Sqrt(gamma*p(n+1,i)/rho(n+1,i))
            if (vc.gt.vmax) then ; vmax = vc; endif

 
 
 [ Сообщений: 9 ] 


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group