navier_gl.txt Information on numerical transient simulation of water flow This simulation program computes and displays the transient solution of water problems using Navier-Stokes and classical physics equations. Simulation is in three dimensional Cartesian coordinates with time stepping, basically a finite element analysis. Water is assumed incompressible and the top surface open to the atmosphere. The atmosphere is taken at nominal values. Viscosity of water is taken into account but temperature change is not computed. Problems of interest have water from 1 centimeter to 10 meters deep with the surface open to the standard atmosphere. Falling or partially submerged solid objects may be moving in air or water. A vacuum may be created due to separation. The model uses elementary volumes fixed in space. The fluid may move from one element to an adjacent element at any time step. A fraction "full" is maintained in order to maintain conservation of mass. In general the initial state will have potential energy and no kinetic energy. An object dropped into calm water or a dam bursting or a slouch gate opening converts potential energy into kinetic energy and after a long time results in a state of only potential energy. The final potential energy is expected to be less than the initial potential energy because friction will have caused loss of kinetic plus potential energy due to heating. Test cases include modeling and simulating: A tower of water that collapses. A hole in the water that fills in. A drop of water falling into calm water. Water initially "piled up" in the corner of an open container. (dam bursting model) Water in one half of a "U" tube. Rotational flow, vortex, whirlpool. This is not intended to be a general fluid dynamics program. Simulation accuracy considerations The physics equations used in the simulation assume the system is: 1) not moving in inertial space or is in a uniform velocity frame of reference 2) at nominal altitude, pressure, temperature and densities 3) a relatively small mass of water 4) energy can be lost due to friction Notes: 1) not moving in inertial space implies that the rotation of the earth is not taken into account. i.e. no Coriolis effect causing water to rotate as is flows from a bowl into a drain. This could be taken into account by providing the latitude of the system being simulated and adding the force computations. 2) at nominal altitude, pressure, temperature and densities means that the following parameters depend on the location of the system being simulated but are currently constants. Acceleration due to gravity is constant throughout the volume being analyzed. Atmospheric pressure, temperature, density and viscosity are constant throughout the volume being analyzed. Water is incompressible and thus has constant density and viscosity throughout the volume being analyzed. (see choice of constants below) 3) a relatively small mass of water that does not take into account the gravitational effects of the moon, as in tides. Input would need latitude, longitude, date, and time at start of simulation in order to handle gravitational effects. 4) energy can be lost due to friction temperature is not computed and temperature is not used in calculations, thus there can be an apparent loss of kinetic plus potential energy. The numerical approximations used in the simulation introduce quantization errors and numerical errors: 1) quantization of Cartesian coordinates, user input, manual and automatic adjustment. 2) quantization of time step of transient solution, user input of largest time step with automatic reduction as necessary 3) approximation of system geometry, user input 4) finite difference approximation of partial derivatives 5) numerical errors due to roundoff and approximation of functions 6) ignoring chaos theory, Brownian movement and the like. only MKS, metric units are used: Meter, Kilogram, Second for: Length, Mass, Time dimensionality distance in meter L time in second T mass in kilogram M velocity in meter/second L/T acceleration in meter/second^2 L/TT area in meter^2 LL volume in meter^3 LLL force in kilogram*meter/second^2 ML/TT pressure in kilogram/(meter*second^2) M/LTT density in kilogram/meter^3 M/LLL viscosity in kilogram/(meter*second) M/LT momentum in kilogram*meter/second ML/T moment of inertia in kilogram*meter^2 MLL gravity acceleration 9.80665 meter/second^2 acceleration air density 1.223 kilogram/meter^3 15.7 'C air pressure 101.32 kilogram/(meter*second^2) sea level air viscosity 182.7E-6 kilogram/(meter*second) 18 'C water density 998.2 kilogram/meter^3 20 'C water dynamic viscosity 0.001066 kilogram/(meter*second) 20 'C water pressure 9690.93 kilogram/(meter*second^2) at 1m water mass 998.2 kilogram for 1 meter cube wall viscosity 0.002 kilogram/(meter*second) user set sound speed in water 1461. meters/sec water air surface tension 7.28 newton/meter An element of volume in the simulation is dx by dy by dz. The Z axis is positive in the vertical direction. The user can rotate the view of the simulation on the screen. The initial orientation is positive X axis to the right, positive Y axis into the screen. An element or cell, sometimes called a box, has six sides with names xn, xp for x-negative side, x-positive side, yn, yp for y-negative side, y-positive side, zn, zp for z-negative side, z-positive side. These names of sides are used later when describing pressures. The elements are in a Cartesian system that has nx * ny * nz elements. The bottom and four sides are solid walls and the top is open to the atmosphere. The elements are fixed in space and do not move. Some amount of water in an element may move to any of the six adjacent neighboring elements in any time step. At each time step, the simulation computes the pressure, positive outward, for every face of every element. Then the forces normal to the faces are computed. Forces due to viscosity from adjacent cell faces at different velocities are computed. Forces are resolved to the center of mass of the cell. Then the acceleration is computed, then velocities are integrated. Then the positions are integrated and the water volume in each element is updated. The numerical calculations are based on the physics equations and specifically the Navier Stokes partial differential equation for fluids with non-zero viscosity. Currently flow is irrotational. The ability for each cell to have rotation about its three axis can be added. An element of cell may contain a solid, air, vacuum or water. The amount of water in an element is the fraction 'full' that may vary from 0.0 to 1.0 times the volume of the element. The initial condition of every element may be set. The physical geometry may be any shape within the system by placing solid elements to create an arbitrary shape. __ / /| /__/ | Consider a single element of water. | | / Compute the static plus dynamic pressure. |__|/ The top, zp, is open to the atmosphere thus the pressure out of the element at zp is air pressure = 101.32 Kg/m^2. The element contains a mass, m = density * dx * dy * dz * full = 998.2 * dx * dy * dz * full. The force of gravity on the mass is F = m g = mass * 9.80665 The pressure at the bottom, zm, is zp + F/(dx*dy) Simplifying zm = zp + 9690.93 * dz * full. Pressure increases linearly with depth, thus xp, xm, yp, ym all are at the average pressure (zm+zp)/2 The water in the element has three signed components of velocity Vx, Vy, Vz The dynamic pressure is added to the static pressure. Pxp = static xp pressure + sign(Vx) * 1/2 density * Vx * Vx Pxm = static xm pressure - sign(Vx) * 1/2 density * Vx * Vx Pyp = static yp pressure + sign(Vy) * 1/2 density * Vy * Vy Pym = static ym pressure - sign(Vy) * 1/2 density * Vy * Vy Pzp = static zp pressure + sign(Vz) * 1/2 density * Vz * Vz Pzm = static zm pressure - sign(Vz) * 1/2 density * Vz * Vz Pressures for all elements are computer from the top to the bottom elements. Each element has six neighbors that may be other elements or walls or open to the atmosphere. A wall has a pressure equal to the pressure of the element face against the wall. The atmosphere has a pressure of 101.32 Kg/m^2. At this instant in the transient solution there can be a different pressure on the neighbors face than on this elements face. The six differences are computed and summed as pressure differential, Pd. The water is assumed incompressible thus Pd must distribute over all six faces. Each face is summed with Pd/6. e.g. DPxp = Pxp + Pd/6 The three signed components of force on the center of mass Fx, Fy, Fx are computed from the pressure differences divided by area. Fx = (DPxp-DPxm)/(dy*dz) Fy = (DPyp-DPxm)/(dx*dz) Fz = (DPzp-DPzm)/(dx*dy) The force due to viscosity and the neighboring element faces moving at different velocities is viscosity * (difference-in-velocity/length) * area The forces are summed with the force due to viscosity and the neighboring face velocity difference. Vz is the Z component of velocity for this element. Vznxm = Vz of neighbor elements xm. There are four surfaces that can have viscous frictional forces in the Z direction ( and four more in X and four more in Y directions). Using mui for the viscosity that depends on the type of neighbor element, the force Fz must be summed with mu1*((Vz-Vznxm)/dz)*dx*dz + mu2*((Vz-Vznxp)/dz)*dx*dz + mu3*((Vz-Vznym)/dz)*dy*dz + mu4*((Vz-Vznyp)/dz)*dy*dz TFz = Fz + muxm*(Vz-Vznxm)*dx + muxp*(Vz-Vznxp)*dx + muym*(Vz-Vznym)*dy + muyp*(Vz-Vznyp)*dy TFx = Fx + muym*(Vx-Vynym)*dy + muyp*(Vx-Vxnyp)*dy + muzm*(Vx-Vxnzm)*dz + muxp*(Vx-Vxnzp)*dz TFy = Fy + muxm*(Vy-Vynxm)*dx + muxp*(Vy-Vynxp)*dx + muzm*(Vy-Vynzm)*dz + muzp*(Vy-Vynzp)*dz Note that the muij may be water viscosity, air viscosity or wall viscosity depending on the neighbor element contents. The acceleration of the elements mass are now Az = TFz/(density*dx*dy*dz*full) from F = m a, a = F / m Ax = TFx/(density*dx*dy*dz*full) Ay = Tfy/(density*dx*dy*dz*full) The new velocity of the center of mass is the numerical integration of acceleration TVz = Vz + Az * dt TVx = Vx + Ax * dt TVy = Vy + Ay * dt The fraction of volume to try to move in this time step is based on face pressure difference. Basically S=v*dt + 1/2 a dt^2 with S/dx being the fraction of the volume. The time step will be decreased if more than 1/4 of the volume is attempting to move or any face is attempting to move more than 1/6 of the volume Szp = (Vz*dt + 1/2 *dPzp*dAxy*dt*dt/dM)/dz Szm = (Vz*dt + 1/2 *dPzm*dAxy*dt*dt/dM)/dz Sxp = (Vx*dt + 1/2 *dPxp*dAyz*dt*dt/dM)/dx Sxm = (Vx*dt + 1/2 *dPxm*dAyz*dt*dt/dM)/dx Syp = (Vy*dt + 1/2 *dPyp*dAxz*dt*dt/dM)/dy Sym = (Vy*dt + 1/2 *dPym*dAxz*dt*dt/dM)/dy Now, mass must be conserved and the actual volume moved must be adjusted so that the total mass in all elements is constant. Additional computation handles the case of a solid object moving through the water elements and special cases of cavitation and air bubbles. The above is repeated for each delta time step, dt. The user sets the initial dt, the program automatically adjusts dt between the users input and min(dx,dy,dz)/speed-of-sound-in-water. The calculations may dynamically double or halve dt as needed. A graphic plot is drawn each user dt time step. The simulation ends at time max_time. Possible decomposition of an element into 8 elements, each with size dx/2 * dy/2 * dz/2 may occur. Test cases: w_one.dat nx=ny=nz=3 a single element of water on the floor in the center w_tower.dat nx=ny=nz=5 a tower of five elements in the center w_hole.dat nx=ny=nz=3 all 2 deep except center 1 deep W_vortex.dat nx=ny=nz=5 four sloping elements around walls w_drop.dat nx=ny=nz=5 a single element at top center w_drop2.dat nx=ny=nz=10 a single element at top center into 2 elements deep navier_gl.dat a sloping wall of water in one corner