Symulacja piłki zbudowanej z sieci punktów materialnych

0

Witam! Mam do napisania program tak jak w temacie dotyczący symulacji piłki zbudowanej z sieci punktów materialnych. W mam symulacji uwzględnić ciśnienie powietrza wewnątrz piłki (zależy jedynie od objętości piłki), które rozpycha punktu na zewnątrz i równoważy siłę ciążenia. Jednak nie mam pomysłu jak ruszyć. Zrobiłem tylko tyle. Prosiłbym o pomoc.

//Wykonać symulację piłki zbudowanej z sieci punktów materialnych.
//   W symulacji uwzględnić ciśnienie powietrza wewnątrz piłki (zależy
//   jedynie od objętości piłki), które rozpycha punktu na zewnątrz
//   i równoważy siłę ciążenia.
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>

#define M_PI 3.14159265358979323846

using namespace std;

double objetosc(double r) //objetosc
{
	return 4/3*M_PI*pow(r,3.0);
}

double gestosc(double m, double V) //gestosc
{
	return m/V;
}

double cisnienie(double V, double n, double R, double T) //cisnienie
{
	return (n*R*T)/V;
}

int main()
{
	double r=0.2386, g=9.81, m=0.65, V, R=8.3144598, T=293, n, p;	
	fstream plik;
	cout << "Symulacji pilki" << endl;
	cout << "" << endl;
	cout << "r - promien kuli (pilki) [m]" << endl;
	cout << "m - masa pilki [kg]" << endl;
	cout << "V - objetosc pilki [m3]" << endl;
	cout << "R - stała gazowa [J/mol*K]" << endl;
	cout << "p - cisnienie pilki [Pa]" << endl;
	cout << "T - temperatura [T]" << endl;
	cout << "n - ilosc gazu [mol]" << endl;
	cout << "" << endl;
	V=objetosc(r);
	while(V<10)
	{
		plik.open("wyniki.txt", ios::out| ios::app);
		p=cisnienie(V, n, R, T);
		cout << " Cisnienie pilki wynosi: " << p << " Pa" << " Objetosc pilki wynosi: " << V << " m3" << endl;	
		plik << p << " " << V << endl;
		V=V+0.1;
		plik.close();
	}
	return 0;
}

0

maras15 napisał(a):

Zrobiłem tylko tyle. Prosiłbym o pomoc.

Czyli nic nie zrobiłeś, bo:

nie masz definicji punktu
nie zdecydowałeś się na algorytm symulacji
twój kod w zasadzie nic nie robi

nawet powiedzenie, że jesteś na początku drogi to już przesada.

Na dodatek po co ci te funkcje objetosc cisnienie? Zdefiniowane w ten sposób mają sens dla idealnej kuli, a nie dla piłki, która się będzie odkształcać.

0

MarekR22 napisał(a):

maras15 napisał(a):

Zrobiłem tylko tyle. Prosiłbym o pomoc.

Czyli nic nie zrobiłeś, bo:

nie masz definicji punktu

nie zdecydowałeś się na algorytm symulacji

twój kod w zasadzie nic nie robi

nawet powiedzenie, że jesteś na początku drogi to już przesada.

Na dodatek po co ci te funkcje objetosc cisnienie? Zdefiniowane w ten sposób mają sens dla idealnej kuli, a nie dla piłki, która się będzie odkształcać.

Czyli jak mam zacząć? Jakieś wskazówki?

0

to na początek:
https://pl.wikipedia.org/wiki/Algorytm_Verleta
https://pl.wikipedia.org/wiki/Algorytm_skokowy

Poza tym masz to zrobić 2d czy 3d? Bo licznie objętości i wpływu ciśnienia dla 3d jest dość skomplikowane.

0

MarekR22 napisał(a):

to na początek:
https://pl.wikipedia.org/wiki/Algorytm_Verleta
https://pl.wikipedia.org/wiki/Algorytm_skokowy

Poza tym masz to zrobić 2d czy 3d? Bo licznie objętości i wpływu ciśnienia dla 3d jest dość skomplikowane.

Oczywiście w 2d wystarczy.

0

Akurat algorytm Verleta ostatnio zaimplementowałem. Tylko nie wiem czy to się przyda? Jeśli tak to co mam dalej zrobić?

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <ctime>
#include <cmath>
#include <iomanip>

#define N 10

using namespace std;

void writePDB(double r[][3], char nazwaPDB[])
{
	int i;
	fstream plik_pdb;
	plik_pdb.open(nazwaPDB, ios::out | ios::app);
	if(plik_pdb == NULL)
	{
		perror ("NIE MOGE OTWORZYC PLIKU PDB! ABY DOPISAC DO NIEGO DANE");
	}

	for(i=0;i<N;i++)
		plik_pdb << "ATOM" << "   " << fixed << setprecision(4) << i+1 << "  " << "AR" << "  " << "AR" << fixed << setprecision(7) << i+1 << "    " 
		<< fixed << setprecision(3) << r[i][0] << "\t" << fixed << setprecision(3) << r[i][1] << "\t" << fixed << setprecision(3) << r[i][2] << "\t" 
		<< "  1.00  0.00      AAA AR" << endl;
	plik_pdb << "END" << endl;
	plik_pdb.close();
}

double Ekin(double v[N][3],float Mas) //funkcja do mierzenia sredniej Energii kin
{
	int i;
	double Ek=0; //energia kinetyczna 
	for(i=0;i<N;i++)
		Ek+=Mas*0.5*(v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
	return Ek;
}

double Eelec(double r[N][3], double qi[N], double qj[N]) //energia elektrostatyczna dla ladunkow losowych 
{
	int i,j;
	double pi=M_PI, epsilon0=8.85*pow(10.0,-19.0), Eelec;
	for(i=1;i<N;i++)
		for(j=1;j<N;j++)
		{
			qi[i]=rand()%21-10;
			qj[j]=rand()%21-10;
			Eelec=qi[i]*qj[j]/4*pi*epsilon0*r[i][j];
		}
	return Eelec;
}

void printVel(double v[N][3])
{
	int i;
	fstream plikV;
	plikV.open("Vel.dat", ios::out | ios::app);
	if(plikV == NULL)
	{
		perror ("NIE MOGE OTWORZYC PLIKU PDB! ABY DOPISAC DO NIEGO DANE");
	}
	
	for(i=0;i<N;i++)
		plikV << v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2] << endl;
	plikV.close();
}


void randVel(double v[N][3], double vMax)
{
	int i, j;	
	for (i = 0; i < N; i++) 
		for (j = 0; j < 3; j++)
			v[i][j] = (vMax * 2.0 * rand() / (RAND_MAX - 1.0))-vMax; //predkooci na razie losowo
}

void Acc(const double r[N][3],double a[N][3], float Mas, float sigma, float eps)
{
	int i,j;
	double F;
	double dx, dy, dz, dr;
	
	for (i = 0; i < N; i++) 
		for (j = 0; j < 3; j++)
			a[i][j] = 0; //zerowanie sil

		for (i = 0; i < N; i++)
		{
			for (j = i +1; j < N; j++) 
			{
				dx=0;
				dy=0;
				dz=0;
				dx = r[i][0] - r[j][0];
				dy = r[i][1] - r[j][1];
				dz = r[i][2] - r[j][2];
				dr=sqrt(dx*dx+dy*dy+dz*dz);
				F= eps*(48 * pow(dr/sigma, -13.0) - 24 * pow(dr/sigma, -7.0));
				a[i][0]+=F*dx/(Mas * dr);
				a[j][0]-=F*dx/(Mas * dr);
				a[i][1]+=F*dy/(Mas * dr);
				a[j][1]-=F*dy/(Mas * dr);
				a[i][2]+=F*dz/(Mas * dr);
				a[j][2]-=F*dz/(Mas * dr);
			
			}
 		}
}


int main()
{
	double r[N][3]={{-7.575,12.030,-3.130},{11.785,-7.800,3.510},{-8.475,-12.030,-5.110},{-9.925,2.290,13.550},{9.465,-0.910,-5.300},{6.865,4.680,-13.550},{5.295,-5.840,-10.000},{-3.615,2.290,-8.100},{2.295,4.160,2.140},{-11.785,4.860,-8.300}};           // polozenia N-indeksuje atomy, x->0,y->1,z->2
	double v[N][3]={{-1.208,3.162,0.881},{-5.775,-2.636,1.659},{2.379,-1.410,0.182},{1.587,-0.773,0.151},{0.040,-2.725,-2.663},{-0.328,-0.097,-5.745},{-3.656,1.071,1.031},{-1.526,-0.832,0.607},{3.428,2.849,2.648},{5.057,1.391,1.249}};           // predkosci N-indeksuje atomy, x->0,y->1,z->2
	double a[N][3]={0};           // przyspieszenia N-indeksuje atomy, x->0,y->1,z->2
	int i, j, k;
	double vMax=5.775;   // A/ps
	double dt=0.001;  // timestep krok czasowy w ps
	int TIME=10000; //-number of timesteps 
	float L=14; //PBC rozmiar skrzynki 2*L x 2*L x 2*L
	float sigma=3.81638, eps=0.23846; //parametryVDW
	float Mas= 39.95; //MASA
	char nPDB[20]="AR_MD.pdb"; //nazwa pliku PDB
	double kb=1.38064852;  //Stala boltzmanna w J/K
	double Ek; //energia kinetyczna
	double Ec; //energia elektrostatyczna
	double qi[N], qj[N];

	fstream plik_pdb, data, Ecsave;

	srand(time(NULL));	
	
	cout << "Dynamika molekularna Argonu na podstawie algorytmu velocity Verlet rozszerzona o obliczenie oddzialywan elektrostatycznych " << endl;
	cout << endl;
	cout << "W pliku Multi Frame PDB zapisywane sa trajektorie wynikowe" << endl;
	cout << endl;
	cout << "W pliku data z rozszerzeniem dat zapisywane sa informacje o kroku czasowym, energii kinetycznej" << endl;
	cout << endl;
	cout << "W pliku Ecsave z rozszerzeniem dat zapisywane sa informacje o kroku czasowym i oddzialywaniach elektrostatycznych" << endl;
	
	plik_pdb.open(nPDB, ios::out | ios::app);
	if(plik_pdb == NULL)
	{
		perror("NIE MOGE OTWORZYC PLIKU PDB!");
		return 1;
	}
	plik_pdb.close();
	
	data.open("data.dat", ios::out | ios::app);
	if(plik_pdb == NULL)
	{
		perror("NIE MOGE OTWORZYC PLIKU PDB!");
		return 1;
	}

	randVel(v,vMax); //generowanie losowych predkosci

	writePDB(r,nPDB); //zapis pierwszych koordynatow
	Ecsave.open("Ec.dat", ios::out | ios::app);
	for(k=1;k<=TIME;k++){
		Acc(r,a, Mas, sigma, eps);
		for (i = 0; i < N; i++)
			for (j = 0; j < 3; j++) 
			{
				r[i][j] += v[i][j] * dt + 0.5 * a[i][j] * dt * dt;
				// PBC
				if (r[i][j] < -L) r[i][j] += 2*L;
				if (r[i][j] > L) r[i][j] -= 2*L;
				v[i][j] += 0.5 * a[i][j] * dt;
			}
		Acc(r,a, Mas, sigma, eps);
		for (i = 0; i < N; i++)
			for (j = 0; j < 3; j++)
				v[i][j] += 0.5 * a[i][j] * dt;
			if(k%10==0)
			{
				writePDB(r,nPDB);
				Ek=Ekin(v,Mas);
				Ec=Eelec(r, qi, qj);
				//cout << Ec << endl;
				data << k << fixed << setprecision(6) << "\t" << Ek  << fixed << setprecision(6) << "\t" << 2.0*1.66054*Ek/(3*N*kb) << "\t"
				<< endl; 
				Ecsave << k << "\t" << Ec << endl; 
				//2.0*1.66054*Ek/(3*N*kb) bierze sie z przeliczenia EK na J (1.66054) oraz podstawienia kb w J/K (1.38064852) 10^34 w liczniku i mianowniku skraca sie
			}
	}
	Ecsave.close();
	data.close();
	cout << "GOTOWE!" << endl;
	system("Pause");
	return 0;
}

0

Nie wiem czy się zgodzisz ale takie rzeczy to chyba najłatwiej pisać używając paradygmatu obiektowego (czy jak to się nazywa).

Stwórz sobie klasę punkt która będzie zawierać własności takie jak położenie, prędkość, przyspieszenie oraz metodę która ustawi siłę na nią działającą (na tej podstawie wyliczysz przyspieszenie potem prędkość i położenie w funkcji czasu, zwykłe prawa newtona + małe przedziały czasu i powinno wyjść w miarę dokładnie, ale jak wolisz dokładniej to stosuj bardziej skomplikowane algorytmy).

class Point
{
    vector2d position;
    vector2d velocity;
    vector2d acceleration;

public:
    void appendForce(vector2d force);
}

Następnie stwórz klasę Piłka która będzie się składała z tablicy takich punktów a w niej będziesz liczył siłę dla każdego punktu.

Symulacja piłki zbudowanej z sieci punktów materialnych

maras15 napisał(a):

MarekR22 napisał(a):

maras15 napisał(a):

MarekR22 napisał(a):

1 użytkowników online, w tym zalogowanych: 0, gości: 1

Praca dla programistów

Forum dyskusyjne

Sprawy administracyjne

O nas

Skontaktuj się z nami