Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- //
- // main.cpp
- // Kolejna_proba
- //
- // Created by Syn Kury on 13/10/2019.
- // Copyright © 2019 Syn Kury. All rights reserved.
- //
- #include <iostream>
- #include<cmath>
- //************************************
- bool czyZgodne(double tab1[], double tab2[], const int n)
- {
- for(int i=2;i<n-1;i+=2)
- {
- if(tab1[i]-tab2[i]<0.000000000000001 and tab1[4]-tab2[4]) return true;
- else return false;
- }
- return false;
- }
- //************************************
- int main(int argc, const char * argv[]) {
- double ro = 1.0;
- double ni = 23.13 * 1e-6;
- double g = 10;
- double dd = 0.1;
- const int n = 41; //ilość węzłów
- double h = dd / n;
- double dt = 0.0001;
- double u[n]{}; //wektor prędkości w kroku f
- double u1[n]{}; //wektor prędkości w kroku f+1
- double p[n]{}; //wektor ciśnienia w kroku f
- double p1[n]{}; //wektor ciśnienia w kroku f+1
- double D[n]{};
- double S[n]{};
- int iteracja=1;
- //************************************
- // zadanie warunków brzegowych
- u[0]=10;
- u[n-1]=10;
- u1[0]=10;
- u1[n-1]=10;
- //************************************
- // zadanie prędkości początkowych
- for(int i=2;i<n-1;i+=2)
- {
- u[i]=0;
- }
- //************************************
- // rozpoczęcie algorytmu iteracyjnego
- do
- {
- //************************************
- // obliczenie prędkości w węzłach k=2,4...,n-2
- for(int i=2;i<n-1;i+=2)
- {
- u1[i]=u[i]+ni*dt*(u[i-2]-2*u[i]+u[i+2])/(4*h*h)-dt*((u[i]+u[i+2])*(u[i]+u[i+2])-(u[i-2]+u[i])*(u[i-2]+u[i]))/(8*h)-dt/ro*(p[i+1]-p[i-1])/(2*h)+dt*g;
- }
- //************************************
- // obliczenie wektora D oraz S w węzłach k=3,5,...n-3
- for(int i=3;i<n-4;i+=2)
- {
- D[i]=u[i+1]-u[i-1];
- S[i]= (-((u[i+1]+u[i+3])*(u[i+1]+u[i+3]))+2*(u[i-1]+u[i+1])*(u[i-1]+u[i+1])-((u[i-3]+u[i-1])*(u[i-3]+u[i-1])))/(8*h);
- }
- D[1]=D[3];
- S[1]=S[3];
- D[n-2]=D[n-4];
- S[n-2]=D[n-4];
- //************************************
- // obliczenie wektora p w węzłach k=3,5,...n-3
- for(int i=3;i<n-4;i+=2)
- {
- p1[i]=0.5*(p[i+2]+p[i-2])-h*ro/dt*D[i]-ni*ro/(4*h)*(D[i-2]-2*D[i]+D[i+2]);
- }
- p1[1]=p[3];
- p[n-2]=p[n-4];
- //************************************
- // korekcja prędkości
- for(int i=2;i<n-1;i+=2)
- {
- u1[i]=u[i]+ni*dt*(u[i-2]-2*u[i]+u[i+2])/(4*h*h)-dt*((u[i]+u[i+2])*(u[i]+u[i+2])-(u[i-2]+u[i])*(u[i-2]+u[i]))/(8*h)-dt/ro*(p1[i+1]-p1[i-1])/(2*h)+dt*g;
- }
- //************************************
- // sprawdzenie warunku zaczymania algorytmu
- if(czyZgodne(u1, u, n))
- {
- break; //wyjście z pętli
- }
- else
- {
- for(int i=0;i<n;i++)
- {
- u[i]=u1[i]; //przypisanie wartości z kolejnego kroku (f+1) do tablicy (f)
- }
- }
- iteracja++;
- } while (iteracja!=1000000);
- //************************************
- // wyświetlenie wyników
- std::cout<<iteracja<<std::endl;
- for(int i=0;i<n;i+=2)
- {
- std::cout<<(u1[i])<<std::endl;
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement