Математическое моделирование процесса удара твердого тела о водную поверхность.

Два из них равны скорости звука левого и правого состояний.где гамма — адиабатический показатель. Первый из них — положение начала волны разрежения, второй — скорость распространения ударной волны.Обозначим:Состояния после скачка (область 4) связаны условиями скачка Ренкина-Гюгонио.Но чтобы рассчитать плотность в области 4, нужно знать давление в этойобласти. Это связано с разрывом контакта с давлением в области 3 соотношениемP4 = P3.Давление в области 3 можно рассчитать только итеративно, правильное решение находится, когдаu3 = u4.u3 – u4 = 0.Эту функцию можно вычислить с произвольной точностью, что дает давление в области 3.P3= calculate (P3,s,s,,)наконец можем рассчитатьu4 = u3и ρ33.2. Программная реализация в среде МатЛабИнициализация переменных задается функцией function y = sod_func(P)%определяются функции для дальнейшего использования в analytic%начальные условияrho_l = 1;P_l = 1;u_l = 0;rho_r = 0.125;P_r = 0.1;u_r = 0;gamma = 1.4;mu = sqrt( (gamma-1)/(gamma+1) );y = (P - P_r)*(( ((1 - mu^2)^2)*((rho_r*(P + mu*mu*P_r))^-1) )^(0.5))... - 2*(sqrt(gamma)/(gamma - 1))*(1 - power(P, (gamma - 1)/(2*gamma)));endОсновные аналитические соотношения, изложенные выше, реализованы функцией analytic_sod(t), для которой требуется указать время: function [data] = analytic(t)%% | | | | |% | | | | |% | | | | |%___|_______|___|_____|_________|_______________% x1 x2 x0 x3 x4%%требуется ввести время: tifnargin < 1%установить по умолчаниюt = 0.2;end%Начальные условияx0 = 0;rho_l = 1;P_l = 1;u_l = 0;rho_r = 0.125;P_r = 0.1;u_r = 0;gamma = 1.4;mu = sqrt( (gamma-1)/(gamma+1) );%скорость звукаc_l = power( (gamma*P_l/rho_l),0.5);c_r = power( (gamma*P_r/rho_r),0.5);P_post = fzero('sod_func',pi);v_post = 2*(sqrt(gamma)/(gamma - 1))*(1 - power(P_post, (gamma - 1)/(2*gamma)));rho_post = rho_r*(( (P_post/P_r) + mu^2 )/(1 + mu*mu*(P_post/P_r)));v_shock = v_post*((rho_post/rho_r)/( (rho_post/rho_r) - 1));rho_middle = (rho_l)*power((P_post/P_l),1/gamma);%ключевые положения областейx1 = x0 - c_l*t;x3 = x0 + v_post*t;x4 = x0 + v_shock*t;%определяем x2c_2 = c_l - ((gamma - 1)/2)*v_post;x2 = x0 + (v_post - c_2)*t;%начало установки значенийn_points = 1000; %устанавливается пользователем%границы (нужно задать)x_min = -0.5;x_max = 0.5;x = linspace(x_min,x_max,n_points);data.x = x';data.rho = zeros(n_points,1); %плотностьdata.P = zeros(n_points,1); %давлениеdata.u = zeros(n_points,1); %скоростьdata.e = zeros(n_points,1); %внутренняяэнергияfor index = 1:n_pointsifdata.x(index) < x1%решение b4 x1data.rho(index) = rho_l;data.P(index) = P_l;data.u(index) = u_l;elseif (x1 <= data.x(index) && data.x(index) <= x2)%решениеb/wx1 иx2c = mu*mu*((x0 - data.x(index))/t) + (1 - mu*mu)*c_l; data.rho(index) = rho_l*power((c/c_l),2/(gamma - 1));data.P(index) = P_l*power((data.rho(index)/rho_l),gamma);data.u(index) = (1 - mu*mu)*( (-(x0-data.x(index))/t) + c_l);elseif (x2 <= data.x(index) && data.x(index) <= x3)%решениеb/wx2 иx3data.rho(index) = rho_middle;data.P(index) = P_post;data.u(index) = v_post;elseif (x3 <= data.x(index) && data.x(index) <= x4)%решениеb/wx3 иx4data.rho(index) = rho_post;data.P(index) = P_post;data.u(index) = v_post;elseif x4 < data.x(index)%решение для x4data.rho(index) = rho_r;data.P(index) = P_r;data.u(index) = u_r;enddata.e(index) = data.P(index)/((gamma - 1)*data.rho(index));endendДля выполнениязадачи введем в командной строке следующие комнды:time = 0.1;data = analytic_sod(time);figure,subplot(2,2,1),plot(data.x,data.rho,'-b','LineWidth',2);xlabel('x (m)');ylabel('Плотность (кг/м^3)');title('Плотность как функция позиции');gridon;subplot(2,2,2),plot(data.x,data.P,'-g','LineWidth',2);xlabel('x (m)');ylabel('Давление (Пa)');title('Функциядавленияотпозиции');gridon;subplot(2,2,3),plot(data.x,data.u,'-r','LineWidth',2);xlabel('x (м)');ylabel('Скорость (м/с)');title('Скорость как функция позиции');gridon;subplot(2,2,4),plot(data.x,data.e,'-k','LineWidth',2);xlabel('x (м)');ylabel('Внутренняя энергия (Дж/кг)');title('Внутренняяэнергиякакфункцияпозици');gridon;Получаем результат:Точки х1 и х2 представляют собой расположение головы и хвоста волны разрежения (движение влево). Хотя решение является непрерывным в этой области (область 2), некоторые из производных величин жидкости могут не быть непрерывными. Точка x3 — это положение, которого элемент жидкости, первоначально находившийся в точке x0, достиг к моменту времени t.Точка х3 называется контактным разрывом. Видно, что поперек контактного разрыва давление и нормальная составляющая скорости непрерывны. Однако плотность и удельная энергия не являются непрерывными на контактном разрыве.Точка х4 — местоположение ударной волны (движение вправо). Во время скачка все величины (p, m, e и p) вообще будут разрывными.Точка x3 — это положение, которого элемент жидкости, первоначально находившийся в x0, достиг к моменту времени t. Точка х3 называется контактным разрывом. Видно, что поперек контактного разрыва давление и нормальная составляющая скорости непрерывны. Однако плотность и удельная энергия не являются непрерывными на контактном разрыве.Точка х4 — местоположение ударной волны (движение вправо). Во время скачка все величины (p, m, e и p) вообще будут разрывными.Список использованной литературыШелкович В. М. Условия Ренкина—Гюгонио и балансовые законыдля δ-ударных волн, Фундамент. и прикл. матем., 12:6 (2006), 213–229; J. Math. Sci., 151:1 (2008), 2781–2792Данилов В. Г., Шелкович В. М. Распространение и взаимодействие δ-ударных волн гиперболических систем законов сохранения // Докл. РАН.—2004.—Т. 394, № 1.—С. 10—14.Ковеня В.М., Чирков Д.В. Методы конечных разностей и конечныхобъемов для решения задач математической физики Электронное учебное пособие. — Новосибирск: Новосиб. гос. ун-т, 2013. — 87 с .Годунов С. К. Рябенький В. С. Разностные схемы: Введение в теорию. М.: Наука, 1973. 400 с.Самарский А. А. Теория разностных схем. Изд. 3. М.: Наука, 1989. 616 с.Овсянников Л. В. Лекции по основам газовой динамики. М.: Наука, 1981, 368 с.Годунов С. К. Разностный метод численного расчета разрывных решений: Матем.сборник. 1959. Т. 47. Вып. 3. С. 271-306.Баландин М.Ю. Шурина Э.П. Методы решения СЛАУ большой размерности. Учебное пособие / Новосибирск: Изд-во НГТУ, 2000. 70 с.Пинчуков В. И., Шу Ч.-В. Численные методы высоких порядков для задачаэрогидродинамики. – Новосибирск: Издательство СО РАН, 2000. – 232 с.Jameson A., Schmidt W., Turkel E. Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes. Technical Report AIAA-81-1259, 1981.Lynn J. F..Multigrid Solution of the Euler Equations with Local Preconditioning. PhD thesis, University of Michigan, 1995.Roe P. L. Approximate Riemann solvers, parameter vectors and difference schemes // J. of Comp. Physics. – 1981. – Vol. 43. – P. 337-372.. Toro E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics: 2nd Ed., SpringerVerlag, 1999.Weiss J. M., Maruszewski J. P., Smith W. A. Imlicit Solution of Preconditioned NavierStokes Equations Using Algebraic Multigrid // AIAA J. - 1999. - Vol 37, № 1. - P. 29-36.Anderson W. K., Thomas J. L., van Leer B. Comparison of Finite Volume Flux Vector Splittings for the Euler Equations // AIAA J. – 1986. – Vol 24, №9. – P. 1453-1460.Чёрный С. Г., Чирков Д. В., Лапин В. Н. и др. Численное моделирование течений в турбомашинах. – Новосибирск: Наука, 2006. – 206 с.