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)
- {
- double delta = 1e-3;
- if (std::abs(tab1[0] - tab2[2]) <= abs(delta * tab1[0])
- and
- std::abs(tab1[n-1]-tab2[n-3] <=abs(delta * tab1[n-1])))
- {
- return true;
- }
- else {
- return false;
- }
- }
- //************************************
- int main(int argc, const char* argv[]) {
- double ro = 1.225;
- 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 a = -(dt / (2 * h));
- double b = dt / (h);
- 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 deltp[n]{};
- double D[n]{};
- double S[n]{};
- long iteracja = 0;
- long iterp = 0;
- double maxPressureIteration = 10000;
- double maxVelocityIteration = 10000;
- //************************************
- // 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 - 2; 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 - 2; 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 - 3; 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
- iterp = 0;
- do
- {
- p[1] = p[3];
- p[n - 2] = p[n - 4];
- /* for(int i=3;i<n-3;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]);
- }*/
- for (int i = 3; i < n - 3; i += 2)
- {
- deltp[i] = (-a * (deltp[i - 2] + deltp[i + 2]) - (u1[i + 1] - u1[i - 1])) / b;
- }
- for (int i = 3; i < n - 3; i += 2)
- {
- p1[i] = p[i] + deltp[i];
- }
- for (int i = 3; i < n - 3; i += 2) {
- if (u1[i + 1] - u1[i - 1] < 0.00000001)
- {
- break;
- }
- else {
- for (int i = 3; i < n - 3; i += 2) {
- p[i] = p1[i];
- }
- }
- }
- ++iterp;
- } while (iterp != maxPressureIteration);//max->2
- //p[1]=p[3];
- //p[n-2]=p[n-4];
- //p[i]=p1[i];
- //************************************
- // korekcja prędkości
- for (int i = 2; i < n - 2; 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++;
- std::cout << "Licze dalej " << iteracja << "/" << maxVelocityIteration << std::endl;
- } while (iteracja != maxVelocityIteration);
- //************************************
- // wyświetlenie wyników
- std::cout << "Po osiągnięciu " << iteracja << " iteracji uzyskano następujące wyniki prędkości: " << std::endl;
- for (int i = 0; i < n; i += 2)
- {
- //std::cout << round(u1[i]) << std::endl;
- std::cout << u1[i] << std::endl;
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement