Wassersimulation (SPH)

[Ikaron]

Hallo Leute, ich versuche gerade mit Hilfe dieses Artikels eine Wassersimulation zum Laufen zu bringen. Der Einfachheit halber lasse ich Viskosität und Overflächenspannung erstmal weg, es geht einzig und allein um die Kräfte, die aus Druckunterschieden resultieren.
Hierfür richte ich mich nach Figur 2.3:

und möchte diese Formel implementieren. Dafür braucht man allerdings eine Kernelfunktion (Figur 2.8):

Diese muss abgeleitet werden nach x und y (siehe Nabla-Operator), hier wären meine Ergebnisse. Sind diese korrekt?

Leider sind Fehler bei so langen Ableitungstermen häufig und die Schulmathematik reicht dort auch nicht unbedingt aus.

Hier ist noch mein Code bis jetzt:

private Vector2f calculateForce(int i) {

	Vector2f force = new Vector2f();

	for (int j = 0; j < particles.length; j++) {

		if (i == j)
			continue;

		Particle pi = particles[i];
		Particle pj = particles[j];

		Vector2f pforce = kernel(pi.getPosition().sub(pj.getPosition())
				.length(), h);
		pforce = pforce.mul(pj.getMass() * (pi.getPressure() + pj.getPressure())
				/ pj.getDensity());
		force = force.add(pforce);
	}

	return force.mul(-1);
}

private static final float pi = 3.1415926535f;

private Vector2f kernel(float x, float y) {

	float ypow5 = y * y * y * y * y;

	return new Vector2f((45 * x * x - 90 * x * y + 45 * y * y) / (pi * y * y * y * y * y * y),
 		15 * (39 * ypow5 * y * y * y - 72 * x * ypow5 * y * y + 27 * x * x * ypow5 * y - 6 * x * x * x * ypow5) / (pi * ypow5 * ypow5 * y * y));
	}

Wie zu erwarten war, stimmt das ganze so noch nicht. Ich mach mal 2 Videos, eins mit den Druckkräften und eins ohne, damit ihr das mal seht.

Video

Vielen Dank für eure Zeit.

 

[rme]

Ich kann deine Berechnung nicht so ganz nachvollziehen. Woher hast du plötzlich x und y? Die Funktion ist doch zweistellig mit r und h. Weiterhin wendest du den Differentialoperator auf die Funktion an und nicht den partiellen Differentialoperator, der hier erforderlich wäre.

 

[Ikaron]

Ich habe das nur der Einfachheit halber x und y statt r und h genannt. Desweiteren dachte ich, dass d/dx bereits der partielle Differenzialoperator ist - Man leitet das Feld in x-Richtung ab für einen gegebenen y-Wert.
Ich muss allerdings zugeben, dass ich den Nabla-Operator nur ansatzweise verstanden habe - Daran wird's wohl liegen.