решение уравнеия растекания нефти

@Жананра · Регистрация: 24.04.2013

Студворк — интернет-сервис помощи студентам

помогите написать программу вычисляющую радиус нефтяного пятна с заданным уравнением растекания.
все исходные данные в фаиле:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
function drdt=spread(t,r);
rhow=1020;
rhoo=869;
Voilin=30000;
gr=9.82;
sigma=0.024;
nu=0.000001;
nuinit=0.0001;
delta=(rhow-rhoo)/rhow;
r0=(Voilin/pi)^1/2;
 
drdt=zeros(2,1);
drdt(1)= r(2);
drdt(2)=3/2*delta*Voilin*gr/(pi*r(1)^3)+3*pi*r(1)*sigma/(rhoo*Voilin)+r(2)^2/r(1)-2.175*rhow*nu^1/2*r(1)^3/2*r(2)^3/2/(rhoo*Voilin)-3*nuinit*r(2)/(2*r(1)^2);

Matlab M
1
2
>> T=[0:1:3600];
>> [t,r]=ode45('spread',T,r0);

ну это мой, вариант программы, который не идет.
покажите ошибки, пожалуйста. и желательно вывести график в виде кругов с полученными радиусами

@Жананра · 14.06.2013, 13:17 **[ТС]**

H у меня задается в файл-функции spread.
А если график у меня выходит прямая со значением в каждый момент времени равным начальному, то получается система не решается?

@tvoretsmira · 14.06.2013, 14:00

Сообщение от Жананра

просто параметры: rhoo, Voilin, sigma, nuinit - должны изменяться со временем.

Не вижу чтоб они от времени хоть как-то зависили.

Сообщение от Жананра

А если график у меня выходит прямая со значением в каждый момент времени равным начальному, то получается система не решается?

В конце функции rastekanie наверно так должно быть, так как x имеет размер [количество точек времени на 9]:
plot(t, x(:,1));

@Жананра · 14.06.2013, 15:07 **[ТС]**

Сообщение от tvoretsmira

Не вижу чтоб они от времени хоть как-то зависили.

я потом приняла их на начальные значения, а изменяющиеся заменила на функции х(t)

Сообщение от tvoretsmira

plot(t, x(:,1));

здесь мы берем только первый х?

Добавлено через 1 минуту
да! график возрастающей функции получился!

Добавлено через 1 минуту
теперь буду проверять подходит ли это решение под реальные условия

Добавлено через 6 минут
Помните, вы мне помогли нарисовать круги по заданным r

Matlab M
1
2
3
4
5
6
7
8
9
10
function krugi
T=[0:300:3600];
r0=[22.58, 0];
[t,r]=ode45('spread1',T,r0);
figure
axes
hold on
for ik=1:size(r,1)
circle(0,0,r(ik,1))
end

Что мне надо изменить в этом коде чтоб нарисовать новые круги для полученного х1?

@tvoretsmira · 14.06.2013, 15:19

Всё так само, только r это теперь x

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function rastekanie
global r0 Voilin sigma nuinit rhoo PM H Uair Tend
k14=40000000;
emscst=100;
r0=(0.17*Voilin)^1/2;
Tend=k14./(emscst.*H);
T=[0:300:3600];
x0=[22.58, 0, 830, 0, 1800, 0, 0, 0, 3000];
[t,x]=ode45('spread', T, x0);
figure
axes
plot(t, x(:,1));
 
figure
axes
hold on
for ik=1:size(x,1)
circle(0,0,x(ik,1))
end

Кружки тоже можно раскрасить:

Matlab M
1
2
3
4
5
6
function circle(x,y,r)
ang=0:0.01:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);
plot(x+xp,y+yp,'Color',rand(1,3));
end

@Жананра · 15.06.2013, 13:49 **[ТС]**

Еще вопрос: почему вычисляя r0 у меня получается 22.58, а программа выдает 225?

Добавлено через 11 минут
255

@tvoretsmira · 15.06.2013, 13:56

f как вы вычисляете r0? не так как в программе?
(0.17*Voilin)^(1/2) не равно (0.17*Voilin)^1/2=1/2*(0.17*Voilin)^1

@Жананра · 17.06.2013, 10:04 **[ТС]**

Как найти значение t, для которого х(9)=0.

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
function dxdt=spread(t,x);%задаем название функции%
global r0 Voilin sigma nuinit rhoo PM H Uair Tend %общие переменные%
rhow=1020;%плотность воды% 
rhoo=830;%начальная плотность нефти%
Voilin=3000;%объем вылитой нефти%
gr=9.82;%гравитационная постоянная%
sigma=0.024;%начальное суммарное поверхностное натяжение%
nu=0.000001;%вязкость воды%
nuinit=0.00001;%начальная вязкость нефти%
k88=0.5e-10;%константы%
k89=5.0e-6;
k14=40000000;
k15=3.0e-6;
k16=1;
k18=0.01;
c3=3.0e-10;
c5=1.0e-8;
c6=1.0e-10;
fractn=0.4;
kev=0.000000012;
emscst=100;
n=0.25; %коэффициент турбулентности%
PM=1800; %начальное значение РМ%
H=0.25; %высота волны%
Uair=4; %скорость ветра%
alpha=(2-n)/(2+n);
beta=n/(2+n);
r0=(0.17*Voilin)^(1/2);%начальный радиус%
Tend=k14/(emscst*H); %конечное время эмульсификации%
 
dxdt=zeros(9,1);
dxdt(1)= x(2);
dxdt(2)=3*gr*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)))*(rhow-x(3))/(2*pi*rhow*x(1)^3)+3*pi*x(1)*sigma*x(9)/(x(3)*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)))*(Voilin-(x(4)+x(6)+x(7)+x(8))))+x(2)^2/x(1)-2.175*rhow*nu^1/2*x(1)^3/2*x(2)^3/2/(x(3)*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8))))-3*nuinit*x(2)/(2*x(1)^2)*(1-4*(Voilin-(x(4)+x(6)+x(7)+x(8)+x(9)))/(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8))))^(-4);
dxdt(3)=k88+k89*dxdt(4)/(Voilin-(x(4)+x(6)+x(7)+x(8)))-(4*(dxdt(4)+dxdt(6)+dxdt(7)+dxdt(8)+dxdt(9))*(rhow-x(3))+(dxdt(4)+dxdt(6)+dxdt(7)+dxdt(8))*(rhoo-x(3)))/(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)));
dxdt(4)=k18*fractn*0.031*kev*2^(2-beta)*Uair^alpha*x(5)*x(1)^(2-beta);
dxdt(5)=-PM*dxdt(4)/(fractn*x(9));
dxdt(6)=c5*H*x(9);
dxdt(7)=k15*k16*H*x(9);
dxdt(8)=(c6+c3)*x(9);
dxdt(9)=-(Voilin-(x(4)+x(6)+x(7)+x(8)))/Tend-dxdt(4)-dxdt(6)-dxdt(7)-dxdt(8);

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
function rastekanie
global r0 Voilin sigma nuinit rhoo PM H Uair Tend
k14=40000000;
emscst=100;
Voilin=3000;
H=0.25;
r0=(0.17.*Voilin).^(1/2)
Tend=k14./(emscst.*H)
T=[0:3600];
x0=[22.58, 0, 830, 0, 1800, 0, 0, 0, 3000];
[t,x]=ode45('spread', T, x0);
figure
axes
plot(t, x(:, 1));

@tvoretsmira · 17.06.2013, 14:15

Добавить в конце в rastekanie:

Matlab M
1
2
3
4
figure
axes
plot(t, x(:, 9));
grid on

Из графика видно, что пересечение нуля имеет место при $https://www.cyberforum.ru/cgi-bin/latex.cgi?t=1.3304*10^6$ .

@Жананра · 18.06.2013, 09:44 **[ТС]**

Спасибо! я теперь поняла как выводить графики каждой функции отдельно!

А как мне можно в программу втиснуть условие:
Если х(9) становится <0, то заменять далее х(9) на ноль

@tvoretsmira · 18.06.2013, 14:50

Matlab M
1
x(x(:,9)<0,9)=0;

@Жананра · 20.06.2013, 14:21 **[ТС]**

это в конце spread надо добавить?

@tvoretsmira · 20.06.2013, 14:30

Нет, это в rastekanie полсе решения, в то место начиная с которого должна произойти замена. Тоесть если построить график х(:,9) до замены, то на нём будут отрицательные числа. Если сделать замену перед построением графика, то и на графике вместо отрицательных будут нули.

Зосима · 20.06.2013, 16:12

Не по теме:

Сообщение от tvoretsmira

x(x(:,9)<0,9)=0;

(. Y .) не сразу сообразил, что 0 и 9 это не 0.9 :D

@Жананра · 28.08.2013, 17:44 **[ТС]**

tvoretsmira, здравствуйте!
скажите пожалуйста: как можно узнать значение t, для которого х(9)=0?
вот конечный вариант программы:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 function dxdt=spread(t,x);%задаем название функции%
global r0 Voilin sigma nuinit rhoo PM H Uair Tend %общие переменные%
rhow=1020;%плотность воды, кг/м3% 
rhoo=830;%начальная плотность нефти, кг/м3%
Voilin=3000;%объем вылитой нефти, м3%
gr=9.82;%гравитационная постоянная, м/с2%
sigma=0.024;%начальное суммарное поверхностное натяжение, Н/м%
nu=0.000001;%вязкость воды, м2/с%
nuinit=0.00001;%начальная вязкость нефти, м2/с%
k88=0.5e-10;%константы%
k89=5.0e-6;
k14=40000000;
k15=3.0e-6;
k16=1;
k18=0.01;
c3=3.0e-10;
c5=1.0e-8;
c6=1.0e-10;
fractn=0.4;
kev=0.000000012;
emscst=100;
n=0.25; %коэффициент турбулентности%
PM=1800; %начальное значение РМ%
H=0.25; %высота волны, м%
Uair=4; %скорость ветра, м/с%
alpha=(2-n)/(2+n);
beta=n/(2+n);
r0=(0.17*Voilin)^(1/2);%начальный радиус, м%
Tend=k14/(emscst*H); %конечное время эмульсификации, с%
 
dxdt=zeros(9,1);
dxdt(1)= x(2);
dxdt(2)=3*gr*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)))*(rhow-x(3))/(2*pi*rhow*x(1)^3)+3*pi*x(1)*sigma*x(9)/(x(3)*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)))*(Voilin-(x(4)+x(6)+x(7)+x(8))))+x(2)^2/x(1)-2.175*rhow*nu^1/2*x(1)^3/2*x(2)^3/2/(x(3)*(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8))))-3*nuinit*x(2)/(2*x(1)^2)*(1-4*(Voilin-(x(4)+x(6)+x(7)+x(8)+x(9)))/(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8))))^(-4);
dxdt(3)=k88+k89*dxdt(4)/(Voilin-(x(4)+x(6)+x(7)+x(8)))-(4*(dxdt(4)+dxdt(6)+dxdt(7)+dxdt(8)+dxdt(9))*(rhow-x(3))+(dxdt(4)+dxdt(6)+dxdt(7)+dxdt(8))*(rhoo-x(3)))/(5*Voilin-4*x(9)-5*(x(4)+x(6)+x(7)+x(8)));
dxdt(4)=k18*fractn*0.031*kev*2^(2-beta)*Uair^alpha*x(5)*x(1)^(2-beta);
dxdt(5)=-PM*dxdt(4)/(fractn*x(9));
dxdt(6)=c5*H*x(9);
dxdt(7)=k15*k16*H*x(9);
dxdt(8)=(c6+c3)*x(9);
dxdt(9)=-(Voilin-(x(4)+x(6)+x(7)+x(8)))/Tend-dxdt(4)-dxdt(6)-dxdt(7)-dxdt(8);

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
function rastekanie
global r0 Voilin sigma nuinit rhoo PM H Uair Tend
k14=40000000; %константы%
emscst=100;
Voilin=3000;%объем вылитой нефти, м3%
H=0.25;%высота волны, м%
r0=(0.17.*Voilin).^(1/2)%начальный радиус, м%
Tend=k14./(emscst.*H)%конечное время эмульсификации, с%
T=[0:Tend];%интервал времени, с%
x0=[r0, 0, 830, 0, 1800, 0, 0, 0, Voilin];%начальные условия%
[t,x]=ode45('spread', T, x0);%решение системы%
x(x(:,9)<0,9)=0;%дополнительное условие%
figure
axes
plot(t, x(:, 1)); %график функции х(1)%
grid on %сетка%

Новые блоги и статьи Все статьи Все блоги /
GitHub Copilot в JetBrains: модель по умолчанию, Fix и контроль MCP and_y87 12.10.2026 10 октября 2026 года GitHub обновил Copilot для JetBrains: добавил настройки стартовой модели, исправления из диагностики и управление автозапуском MCP. Поддержка IDE 2025. 1 прекращена — требуется. . .	GitHub Copilot в JetBrains: модель по умолчанию, Fix и контроль MCP and_y87 12.10.2026 10 октября 2026 года GitHub обновил Copilot для JetBrains: добавил настройки стартовой модели, исправления из диагностики и управление автозапуском MCP. Поддержка IDE 2025. 1 прекращена — требуется. . .	День 5 из 28: Codex подсказывает следующий запрос and_y87 12.10.2026 9 октября 2026 года OpenAI запустила бету Composer predictions — подсказок следующего сообщения в текущем диалоге Codex. Это помощь с формулировкой запроса, а не автоматическое продолжение задачи. . . .	Запустил конкурс "тем и промптов для текстовых квестов созданных почти чисто ИИ" Adler 06.10.2026 Всем привет! За последние три-четыре дня я создал более 16 текстовых квестовых игр используя преимущественно по одному запросу к ИИ на игру. Мне так понравилось смотреть все ветки/ сцены во всех. . .
ИИ не может найти нужный язык в списке Supersumestria 05.10.2026 Я ему даю вот такое изображение и прошу найти и подчеркнуть немецкий язык. Возвращает он вот это: https:/ / i. **********/ vqBWLe2. png Нужную строчку в 3й колонке просто выдумал. . Это. . .	Новая последняя моя музыка в SUNO zorxor 05.10.2026 Здравствуйте, дорогие мои друзья! С большой радостью я хотел бы представить вам свою новую последнею музыку, которую сгенерировала мне по моей просьбе нейросеть SUNO. С уважением, zorxor. Это. . .	Программный домашний кинотеатр russiannick 27.09.2026 Сподобился на программный домашний кинотеатр. В качестве ЯВУ по традиции выбрал js. В помощники взял Яндекс-Алису. Было создано три зала на разные интересы. исторические и ретро сериал Хичкок. . .	Беседа с ИИ о программистах, недопускающих к созданию и правке кода генеративные ИИ и причины этого zorxor 21.09.2026 Раньше я радовался или получал некоторые эмоции, пусть небольшие, но всё же, от самого процесса написания кода, рекомпиляции и запуска, видя постепенное развитие программы и прочее. А теперь лень. . .

@Жананра 0 / 0 / 0 Регистрация: 24.04.2013 Сообщений: 24
	14.06.2013, 13:17 [ТС]
	H у меня задается в файл-функции spread. А если график у меня выходит прямая со значением в каждый момент времени равным начальному, то получается система не решается? 0

@Жананра 0 / 0 / 0 Регистрация: 24.04.2013 Сообщений: 24
	15.06.2013, 13:49 [ТС]
	Еще вопрос: почему вычисляя r0 у меня получается 22.58, а программа выдает 225? Добавлено через 11 минут 255 0

@tvoretsmira 462 / 452 / 56 Регистрация: 28.05.2013 Сообщений: 699
	15.06.2013, 13:56
	f как вы вычисляете r0? не так как в программе? (0.17Voilin)^(1/2) не равно (0.17Voilin)^1/2=1/2(0.17Voilin)^1 1

@Жананра 0 / 0 / 0 Регистрация: 24.04.2013 Сообщений: 24
	18.06.2013, 09:44 [ТС]
	Спасибо! я теперь поняла как выводить графики каждой функции отдельно! А как мне можно в программу втиснуть условие: Если х(9) становится <0, то заменять далее х(9) на ноль 0

@Жананра 0 / 0 / 0 Регистрация: 24.04.2013 Сообщений: 24
	20.06.2013, 14:21 [ТС]
	это в конце spread надо добавить? 0

@tvoretsmira 462 / 452 / 56 Регистрация: 28.05.2013 Сообщений: 699
	20.06.2013, 14:30
	Нет, это в rastekanie полсе решения, в то место начиная с которого должна произойти замена. Тоесть если построить график х(:,9) до замены, то на нём будут отрицательные числа. Если сделать замену перед построением графика, то и на графике вместо отрицательных будут нули. 1