2014 dxdy logo

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

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




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


11/06/11
66
МИФИ
Здравствуйте.

Я пытаюсь запрограммировать одномерный метод крупных частиц первого порядка без модификаций в точности по книжке Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике.
Решил протестировать на задаче Сода (то есть задача о распаде разрыва с начальными условиями слева $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 
Заслуженный участник
Аватара пользователя


15/10/08
11536
Roy Rogers в сообщении #687151 писал(а):
Можно ли это как-то улучшить?

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

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

 Профиль  
                  
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение16.03.2013, 15:10 
Заслуженный участник
Аватара пользователя


11/04/07
1352
Москва
Цитата:
метод Белоцерковского - незамедлительно в печку...

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

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


15/10/08
11536
Боевая колесница ездить, конечно, могёт, но приходит время и вместе с ним танки.

 Профиль  
                  
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение23.03.2013, 13:05 
Заслуженный участник
Аватара пользователя


11/04/07
1352
Москва
Ваши "TVD, ENO, SPH" - убогая реинкарнация по времени, без повышения порядка аппроксимации на разрывах по пространственным координатам( не говоря о излучении и ряду других важных физических процессов), с заметным увеличением вычислительной вязкости и дикой недокументированной дисперсии. К тому же они сейчас уже бордаты.

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


15/10/08
11536
Zai
Думайте как хотите.

 Профиль  
                  
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение08.05.2013, 20:44 


07/05/13
1
я тоже осваиваю метод крупных частиц
Скажите, а в каком компиляторе fortran вы работаете?, потому что у меня при вводе вашего кода кучу ошибок выкидывает

 Профиль  
                  
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение08.05.2013, 21:23 
Заслуженный участник
Аватара пользователя


15/10/08
11536
kor_viktor
У него стандарт 95 или выше.

 Профиль  
                  
 
 Re: Азбучный вопрос по методу крупных частиц
Сообщение16.05.2016, 11:20 


16/05/16
1
Здравствуйте, на основании чего Вами был выбран шаг?
Код:
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 ] 

Модераторы: photon, profrotter, Парджеттер, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group