Advertisement

Fluid Simulation

Started by December 19, 2014 08:23 PM
1 comment, last by Stainless 10 years, 1 month ago

Hi all,

I am trying to implement a particle based fluid simulation in unity3d, implemented in this paper (http://www.csc.kth.se/utbildning/kth/kurser/DD143X/dkand13/Group10Pawel/report/DanielM_report.pdf)

The paper has everything written is pseudocode (although, seems as though there are one or two errors)

I tried today to get it working, and it is very iratic and basically doesnt work at all.

As a note,My implementation uses the Physics2D.OverlapCircleAll to find the neighboring particles instead of the grid system

I have attached my two scripts, one particle manager with everything and then the particle object, also I made this over the last couple hours and it is VERY messy coding

If anyone here could take a look and point towards my errors i would really appreciate it, I also have no real clue on what values to use for the different constants.

Regards

Particle Manager


using UnityEngine;
using System.Collections;
using System.Collections.Generic;
public class ParticleManager : MonoBehaviour {

	public GameObject particlePrefab;
	public float radius =1; // Maximum distance particles affect each other
	public float collisionRadius; // Distance from a wall tat counts collision
	public float p0=1; // rest density
	public float viscosityLinear=1;
	public float viscosityQuadratic=1;
	public float k=0.5f;
	public float k_near=0.5f;
	public Vector2 gravity;

	public List<Particle> particles;
	public List<List<Particle>> neighbors;


	int currentIndex=0;
	// Use this for initialization
	void Start () {
		gravity = Physics2D.gravity;

		particles = new List<Particle>();
		neighbors = new List<List<Particle>>();
	}
	
	// Update is called once per frame
	void Update () {
	
		if (Input.GetButton ("Fire1")) 
		{
			var mousePos = Input.mousePosition;
			mousePos.z = 10;
			
			var objectPos = Camera.main.ScreenToWorldPoint(mousePos);
			GameObject newP = Instantiate(particlePrefab,
			                              objectPos,
			                              Quaternion.identity) as GameObject;
			Particle n = newP.GetComponent<Particle>();
			n.index = currentIndex;
//
//
			particles.Add(n);
			neighbors.Add(new List<Particle>());
			currentIndex++;


		}

		ApplyExternalForce();
		ApplyViscosity();
		AdvanceParticles();
		UpdateNeighbors();
		DoubleDensityRelaxation();
		worldCollision();
		UpdateVelocity();

	}

	void ApplyExternalForce()
	{

		foreach( Particle p in particles)
		{
			p.vel = p.vel + gravity*Time.deltaTime;
			p.vel = p.vel+p.externalForce;
		}
	}
	void ApplyViscosity()
	{
		foreach( Particle p in particles)
		{
			int index = p.index;
			foreach (Particle n in neighbors[index])
			{
				Vector2 v_pn = n.gameObject.transform.position-
					p.gameObject.transform.position;
				float vel_inward = Vector2.Dot(p.vel - n.vel,v_pn);

				if(vel_inward>0)
				{
					float length = v_pn.magnitude;
					vel_inward = vel_inward/length;
					v_pn = v_pn/length;

					float q = length/radius;

					Vector2 I = 0.5f*Time.deltaTime*(1-q)
						*(viscosityLinear*vel_inward+viscosityQuadratic*Mathf.Pow(vel_inward,2))
							*v_pn;
					p.vel = p.vel - I;
				}
			}
		}
	}

	void AdvanceParticles()
	{
		foreach(Particle p in particles)
		{
			p.previousPosition = p.Pos;
			p.Pos =p.Pos + p.vel*Time.deltaTime;
		}
	}

	void UpdateNeighbors()
	{
		foreach(Particle p in particles)
		{
			neighbors[p.index].Clear();
			Collider2D[] particleColliders = Physics2D.OverlapCircleAll(p.gameObject.transform.position,radius,1<<LayerMask.NameToLayer("Particle"));

			foreach(Collider2D c in particleColliders)
			{
				neighbors[p.index].Add(c.gameObject.GetComponent<Particle>());
			}

			neighbors[p.index].Remove(p);
		}
	}

	void DoubleDensityRelaxation()
	{
		foreach(Particle i in particles)
		{
			float p =0;
			float p_near = 0;

			foreach(Particle n in neighbors[i.index])
			{
				float tempn = (i.Pos - n.Pos).magnitude;
				float q = tempn/radius;


					p = p+ Mathf.Pow(1-q,2);
					p_near = p_near+ Mathf.Pow(1-q,3);

			
			}

			float P = k*(p-p0);
			float Pnear = k_near*p_near;
			Vector2 delta=Vector2.zero;

			foreach(Particle n in neighbors[i.index])
			{
				float tempn = (i.Pos - n.Pos).magnitude;

				float q = tempn/radius;

					Vector2 v_pn = (i.Pos-n.Pos)/tempn;

					Vector2 D = 0.5f*Mathf.Pow(Time.deltaTime,2)
						*(P*(1-q) + Pnear*Mathf.Pow(1-q,2))*v_pn;
					n.Pos = n.Pos+D;
					delta = delta - D;
			
			}

			i.Pos = i.Pos +delta;
		}
	}

	void UpdateVelocity()
	{
		foreach(Particle p in particles)
		{
			p.vel = (p.Pos - p.previousPosition)/Time.deltaTime;
		}
	}

	public float collisionSoftness = 0.4f;
	void worldCollision()
	{
		foreach(Particle p in particles)
		{

			float dist = -1*(p.Pos.y + 4);
			if(p.Pos.y <-4)
			{
				p.Pos=p.Pos + collisionSoftness*dist*Vector2.up;
			}
			if(p.Pos.x > 3)
			{
				p.Pos=p.Pos - collisionSoftness*(p.Pos.x-3)*Vector2.right;
			}
			else if(p.Pos.x < -3)
			{
				p.Pos=p.Pos - collisionSoftness*(p.Pos.x+3)*Vector2.right;
			}
		}
	}


}

Particle


using UnityEngine;
using System.Collections;

public class Particle : MonoBehaviour {

	public Vector2 previousPosition;
	public Vector2 vel;
	public int index;

	public Vector2 externalForce;
	void Update()
	{
		rigidbody2D.velocity = Vector2.zero;
	}

	public Vector2 Pos
	{
		get{
			return transform.position;}
		set {

			transform.position = value;
		}
	}

	void OnCollisionEnter2D(Collision2D coll)
	{
		vel = Vector2.zero;
	}
}

Maybe you could narrow in the bug or bugs by setting all variables to 0.0 except one and then seeing what happens. Then you will know which parts (if any)of the code that works and which ones that doesn't. If nothing works, then maybe your collision detection doesn't work.

Advertisement

The first thing that jumps out at me is this LOC

p.vel = p.vel+p.externalForce;

If that is a force then this is way wrong. Check to see if it is a force, an acceleration, or if it is really a velocity.

if it's a force then you need to change this line to


p.vel += (p.externalForce / mass) * Time.deltaTime;

if it's an acceleration then it should be


p.vel += p.extrnalForce * Time.deltaTime;

This topic is closed to new replies.

Advertisement