Translate

Google+ Badge

Follow by Email

Search This Blog

Thursday, October 23, 2014

General Circulation Model

General Circulation Model

From Wikipedia, the free encyclopedia

Climate models are systems of differential equations based on the basic laws of physics, fluid motion, and chemistry. To “run” a model, scientists divide the planet into a 3-dimensional grid, apply the basic equations, and evaluate the results. Atmospheric models calculate winds, heat transfer, radiation, relative humidity, and surface hydrology within each grid and evaluate interactions with neighboring points.[1]
This visualization shows early test renderings of a global computational model of Earth's atmosphere based on data from NASA's Goddard Earth Observing System Model, Version 5 (GEOS-5).

A general circulation model (GCM), a type of climate model, is a mathematical model of the general circulation of a planetary atmosphere or ocean and based on the Navier–Stokes equations on a rotating sphere with thermodynamic terms for various energy sources (radiation, latent heat). These equations are the basis for complex computer programs commonly used for simulating the atmosphere or ocean of the Earth.
Atmospheric and oceanic GCMs (AGCM and OGCM) are key components of global climate models along with sea ice and land-surface components. GCMs and global climate models are widely applied for weather forecasting, understanding the climate, and projecting climate change. Versions designed for decade to century time scale climate applications were originally created by Syukuro Manabe and Kirk Bryan at the Geophysical Fluid Dynamics Laboratory in Princeton, New Jersey.[2] These computationally intensive numerical models are based on the integration of a variety of fluid dynamical, chemical, and sometimes biological equations.

Note on nomenclature

The initialism GCM stands originally for general circulation model. Recently, a second meaning has come into use, namely global climate model. While these do not refer to the same thing, General Circulation Models are typically the tools used for modelling climate, and hence the two terms are sometimes used as if they were interchangeable. However, the term "global climate model" is ambiguous, and may refer to an integrated framework incorporating multiple components which may include a general circulation model, or may refer to the general class of climate models that use a variety of means to represent the climate mathematically with differing levels of detail.

History: general circulation models

In 1956, Norman Phillips developed a mathematical model which could realistically depict monthly and seasonal patterns in the troposphere, which became the first successful climate model.[3][4] Following Phillips's work, several groups began working to create general circulation models.[5] The first general circulation climate model that combined both oceanic and atmospheric processes was developed in the late 1960s at the NOAA Geophysical Fluid Dynamics Laboratory.[2] By the early 1980s, the United States' National Center for Atmospheric Research had developed the Community Atmosphere Model; this model has been continuously refined into the 2000s.[6] In 1996, efforts began to initialize and model soil and vegetation types, which led to more realistic forecasts.[7] Coupled ocean-atmosphere climate models such as the Hadley Centre for Climate Prediction and Research's HadCM3 model are currently being used as inputs for climate change studies.[5] The role of gravity waves was neglected within these models until the mid-1980s. Now, gravity waves are required within global climate models to simulate regional and global scale circulations accurately, though their broad spectrum makes their incorporation complicated.[8]

Atmospheric vs oceanic models

There are both atmospheric GCMs (AGCMs) and oceanic GCMs (OGCMs). An AGCM and an OGCM can be coupled together to form an atmosphere-ocean coupled general circulation model (CGCM or AOGCM). With the addition of other components (such as a sea ice model or a model for evapotranspiration over land), the AOGCM becomes the basis for a full climate model. Within this structure, different variations can exist, and their varying response to climate change may be studied (e.g., Sun and Hansen, 2003).

Modeling trends

A recent trend in GCMs is to apply them as components of Earth system models, e.g. by coupling to ice sheet models for the dynamics of the Greenland and Antarctic ice sheets, and one or more chemical transport models (CTMs) for species important to climate. Thus a carbon CTM may allow a GCM to better predict changes in carbon dioxide concentrations resulting from changes in anthropogenic emissions. In addition, this approach allows accounting for inter-system feedback: e.g. chemistry-climate models allow the possible effects of climate change on the recovery of the ozone hole to be studied.[9]

Climate prediction uncertainties depend on uncertainties in chemical, physical, and social models (see IPCC scenarios below).[10] Progress has been made in incorporating more realistic chemistry and physics in the models, but significant uncertainties and unknowns remain, especially regarding the future course of human population, industry, and technology.

Note that many simpler levels of climate model exist; some are of only heuristic interest, while others continue to be scientifically relevant.

Model structure

Three-dimensional (more properly four-dimensional) GCMs discretise the equations for fluid motion and integrate these forward in time. They also contain parameterisations for processes – such as convection – that occur on scales too small to be resolved directly. More sophisticated models may include representations of the carbon and other cycles.

A simple general circulation model (SGCM), a minimal GCM, consists of a dynamical core that relates material properties such as temperature to dynamical properties such as pressure and velocity. Examples are programs that solve the primitive equations, given energy input into the model, and energy dissipation in the form of scale-dependent friction, so that atmospheric waves with the highest wavenumbers are the ones most strongly attenuated. Such models may be used to study atmospheric processes within a simplified framework but are not suitable for future climate projections.

Atmospheric GCMs (AGCMs) model the atmosphere (and typically contain a land-surface model as well) and impose sea surface temperatures (SSTs). A large amount of information including model documentation is available from AMIP.[11] They may include atmospheric chemistry.
  • AGCMs consist of a dynamical core which integrates the equations of fluid motion, typically for:
    • surface pressure
    • horizontal components of velocity in layers
    • temperature and water vapor in layers
  • There is generally a radiation code, split into solar/short wave and terrestrial/infra-red/long wave
  • Parametrizations are used to include the effects of various processes. All modern AGCMs include parameterizations for:
A GCM contains a number of prognostic equations that are stepped forward in time (typically winds, temperature, moisture, and surface pressure) together with a number of diagnostic equations that are evaluated from the simultaneous values of the variables. As an example, pressure at any height can be diagnosed by applying the hydrostatic equation to the predicted surface pressure and the predicted values of temperature between the surface and the height of interest. The pressure diagnosed in this way then is used to compute the pressure gradient force in the time-dependent equation for the winds.

Oceanic GCMs (OGCMs) model the ocean (with fluxes from the atmosphere imposed) and may or may not contain a sea ice model. For example, the standard resolution of HadOM3 is 1.25 degrees in latitude and longitude, with 20 vertical levels, leading to approximately 1,500,000 variables.

Coupled atmosphere–ocean GCMs (AOGCMs) (e.g. HadCM3, GFDL CM2.X) combine the two models. They thus have the advantage of removing the need to specify fluxes across the interface of the ocean surface. These models are the basis for sophisticated model predictions of future climate, such as are discussed by the IPCC.

AOGCMs represent the pinnacle of complexity in climate models and internalise as many processes as possible. They are the only tools that could provide detailed regional predictions of future climate change. However, they are still under development. The simpler models are generally susceptible to simple analysis and their results are generally easy to understand. AOGCMs, by contrast, are often nearly as hard to analyse as the real climate system.

Model grids

The fluid equations for AGCMs are discretised using either the finite difference method or the spectral method. For finite differences, a grid is imposed on the atmosphere. The simplest grid uses constant angular grid spacing (i.e., a latitude / longitude grid), however, more sophisticated non-rectantangular grids (e.g., icosahedral) and grids of variable resolution[12] are more often used.[13] The "LMDz" model can be arranged to give high resolution over any given section of the planet. HadGEM1 (and other ocean models) use an ocean grid with higher resolution in the tropics to help resolve processes believed to be important for ENSO. Spectral models generally use a gaussian grid, because of the mathematics of transformation between spectral and grid-point space. Typical AGCM resolutions are between 1 and 5 degrees in latitude or longitude: the Hadley Centre model HadCM3, for example, uses 3.75 in longitude and 2.5 degrees in latitude, giving a grid of 96 by 73 points (96 x 72 for some variables); and has 19 levels in the vertical. This results in approximately 500,000 "basic" variables, since each grid point has four variables (u,v, T, Q), though a full count would give more (clouds; soil levels). HadGEM1 uses a grid of 1.875 degrees in longitude and 1.25 in latitude in the atmosphere; HiGEM, a high-resolution variant, uses 1.25 x 0.83 degrees respectively.[14] These resolutions are lower than is typically used for weather forecasting.[15] Ocean resolutions tend to be higher, for example HadCM3 has 6 ocean grid points per atmospheric grid point in the horizontal.

For a standard finite difference model, uniform gridlines converge towards the poles. This would lead to computational instabilities (see CFL condition) and so the model variables must be filtered along lines of latitude close to the poles. Ocean models suffer from this problem too, unless a rotated grid is used in which the North Pole is shifted onto a nearby landmass. Spectral models do not suffer from this problem. There are experiments using geodesic grids[16] and icosahedral grids, which (being more uniform) do not have pole-problems. Another approach to solving the grid spacing problem is to deform a Cartesian cube such that it covers the surface of a sphere.[17]

Flux buffering

Some early incarnations of AOGCMs required a somewhat ad hoc process of "flux correction" to achieve a stable climate (not all model groups used this technique). This resulted from separately prepared ocean and atmospheric models each having a different implicit flux from the other component than the other component could actually provide. If uncorrected this could lead to a dramatic drift away from observations in the coupled model. However, if the fluxes were 'corrected', the problems in the model that led to these unrealistic fluxes might be unrecognised and that might affect the model sensitivity. As a result, there has always been a strong disincentive to use flux corrections, and the vast majority of models used in the current round of the Intergovernmental Panel on Climate Change do not use them. The model improvements that now make flux corrections unnecessary are various, but include improved ocean physics, improved resolution in both atmosphere and ocean, and more physically consistent coupling between atmosphere and ocean models. Confidence in model projections is increased by the improved performance of several models that do not use flux adjustment. These models now maintain stable, multi-century simulations of surface climate that are considered to be of sufficient quality to allow their use for climate change projections.[18]

Convection

Moist convection causes the release of latent heat and is important to the Earth's energy budget. Convection occurs on too small a scale to be resolved by climate models, and hence it must be parameterized. This has been done since the earliest days of climate modelling, in the 1950s. Akio Arakawa did much of the early work, and variants of his scheme are still used,[19] although there are a variety of different schemes now in use.[20][21][22] Clouds are typically parametrized, not because their physical processes are poorly understood, but because they occur on a scale smaller than the resolved scale of most GCMs. The causes and effects of their small scale actions on the large scale are represented by large scale parameters, hence "parameterization". The fact that cloud processes are not perfectly parameterized is due in part to a lack of understanding of clouds, but not due to some inherent shortcoming of the method.[23]

Output variables

Most models include software to diagnose a wide range of variables for comparison with observations or study of processes within the atmosphere. An example is the 1.5-metre temperature, which is the standard height for near-surface observations of air temperature. This temperature is not directly predicted from the model but is deduced from the surface and lowest-model-layer temperatures. Other software is used for creating plots and animations.

Projections of future climate change

File:Animation of projected annual mean surface air temperature from 1970-2100, based on SRES emissions scenario A1B (NOAA GFDL CM2.1).webmPlay media

Projected annual mean surface air temperature from 1970-2100, based on SRES emissions scenario A1B, using the NOAA GFDL CM2.1 climate model (credit: NOAA Geophysical Fluid Dynamics Laboratory).[24]

Coupled ocean–atmosphere GCMs use transient climate simulations to project/predict future temperature changes under various scenarios. These can be idealised scenarios (most commonly, CO2 increasing at 1%/yr) or more realistic (usually the "IS92a" or more recently the SRES scenarios). Which scenarios should be considered most realistic is currently uncertain, as the projections of future CO2 (and sulphate) emission are themselves uncertain.

The 2001 IPCC Third Assessment Report figure 9.3 shows the global mean response of 19 different coupled models to an idealised experiment in which CO2 is increased at 1% per year.[25] Figure 9.5 shows the response of a smaller number of models to more realistic forcing. For the 7 climate models shown there, the temperature change to 2100 varies from 2 to 4.5 °C with a median of about 3 °C.

Future scenarios do not include unknowable events – for example, volcanic eruptions or changes in solar forcing. These effects are believed to be small in comparison to GHG forcing in the long term, but large volcanic eruptions, for example, are known to exert a temporary cooling effect.

Human emissions of GHGs are an external input to the models, although it would be possible to couple in an economic model to provide these as well. Atmospheric GHG levels are usually supplied as an input, though it is possible to include a carbon cycle model including land vegetation and oceanic processes to calculate GHG levels.

Emissions scenarios

In the 21st century, changes in global mean temperature are projected to vary across the world
Projected change in annual mean surface air temperature from the late 20th century to the middle 21st century, based on SRES emissions scenario A1B (credit: NOAA Geophysical Fluid Dynamics Laboratory).[24]

For the six SRES marker scenarios, IPCC (2007:7–8) gave a "best estimate" of global mean temperature increase (2090–2099 relative to the period 1980–99) that ranged from 1.8 °C to 4.0 °C.[26] Over the same time period, the "likely" range (greater than 66% probability, based on expert judgement) for these scenarios was for a global mean temperature increase of between 1.1 and 6.4 °C.[26]

Pope (2008) described a study where climate change projections were made using several different emission scenarios.[27] In a scenario where global emissions start to decrease by 2010 and then decline at a sustained rate of 3% per year, the likely global average temperature increase was predicted to be 1.7 °C above pre-industrial levels by 2050, rising to around 2 °C by 2100. In a projection designed to simulate a future where no efforts are made to reduce global emissions, the likely rise in global average temperature was predicted to be 5.5 °C by 2100. A rise as high as 7 °C was thought possible but less likely.

Sokolov et al. (2009) examined a scenario designed to simulate a future where there is no policy to reduce emissions. In their integrated model, this scenario resulted in a median warming over land (2090–99 relative to the period 1980–99) of 5.1 °C. Under the same emissions scenario but with different modeling of the future climate, the predicted median warming was 4.1 °C.[28]

Accuracy of models that predict global warming

SST errors in HadCM3
North American precipitation from various models.
Temperature predictions from some climate models assuming the SRES A2 emissions scenario.

AOGCMs represent the pinnacle of complexity in climate models and internalise as many processes as possible. However, they are still under development and uncertainties remain. They may be coupled to models of other processes, such as the carbon cycle, so as to better model feedback effects. Most recent simulations show "plausible" agreement with the measured temperature anomalies over the past 150 years, when forced by observed changes in greenhouse gases and aerosols, and better agreement is achieved when both natural and man-made forcings are included.[29][30]

No model – whether a wind-tunnel model for designing aircraft, or a climate model for projecting global warming – perfectly reproduces the system being modeled. Such inherently imperfect models may nevertheless produce useful results. In this context, GCMs are capable of reproducing the general features of the observed global temperature over the past century.[29]

A debate over how to reconcile climate model predictions that upper air (tropospheric) warming should be greater than surface warming, with observations some of which appeared to show otherwise[31] now appears to have been resolved in favour of the models, following revisions to the data: see satellite temperature record.

The effects of clouds are a significant area of uncertainty in climate models. Clouds have competing effects on the climate. One of the roles that clouds play in climate is in cooling the surface by reflecting sunlight back into space; another is warming by increasing the amount of infrared radiation emitted from the atmosphere to the surface.[32] In the 2001 IPCC report on climate change, the possible changes in cloud cover were highlighted as one of the dominant uncertainties in predicting future climate change;[33] see also[34]

Thousands of climate researchers around the world use climate models to understand the climate system. There are thousands of papers published about model-based studies in peer-reviewed journals – and a part of this research is work improving the models. Improvement has been difficult but steady (most obviously, state of the art AOGCMs no longer require flux correction), and progress has sometimes led to discovering new uncertainties.

In 2000, a comparison between measurements and dozens of GCM simulations of ENSO-driven tropical precipitation, water vapor, temperature, and outgoing longwave radiation found similarity between measurements and simulation of most factors. However the simulated change in precipitation was about one-fourth less than what was observed. Errors in simulated precipitation imply errors in other processes, such as errors in the evaporation rate that provides moisture to create precipitation. The other possibility is that the satellite-based measurements are in error. Either indicates progress is required in order to monitor and predict such changes.[35]

A more complete discussion of climate models is provided in the IPCC's Third Assessment Report.[36]
  • The model mean exhibits good agreement with observations.
  • The individual models often exhibit worse agreement with observations.
  • Many of the non-flux adjusted models suffered from unrealistic climate drift up to about 1 °C/century in global mean surface temperature.
  • The errors in model-mean surface air temperature rarely exceed 1 °C over the oceans and 5 °C over the continents; precipitation and sea level pressure errors are relatively greater but the magnitudes and patterns of these quantities are recognisably similar to observations.
  • Surface air temperature is particularly well simulated, with nearly all models closely matching the observed magnitude of variance and exhibiting a correlation > 0.95 with the observations.
  • Simulated variance of sea level pressure and precipitation is within ±25% of observed.
  • All models have shortcomings in their simulations of the present day climate of the stratosphere, which might limit the accuracy of predictions of future climate change.
    • There is a tendency for the models to show a global mean cold bias at all levels.
    • There is a large scatter in the tropical temperatures.
    • The polar night jets in most models are inclined poleward with height, in noticeable contrast to an equatorward inclination of the observed jet.
    • There is a differing degree of separation in the models between the winter sub-tropical jet and the polar night jet.
  • For nearly all models the r.m.s. error in zonal- and annual-mean surface air temperature is small compared with its natural variability.
    • There are problems in simulating natural seasonal variability.[citation needed] ( 2000)
      • In flux-adjusted models, seasonal variations are simulated to within 2 K of observed values over the oceans. The corresponding average over non-flux-adjusted models shows errors up to about 6 K in extensive ocean areas.
      • Near-surface land temperature errors are substantial in the average over flux-adjusted models, which systematically underestimates (by about 5 K) temperature in areas of elevated terrain. The corresponding average over non-flux-adjusted models forms a similar error pattern (with somewhat increased amplitude) over land.
      • In Southern Ocean mid-latitudes, the non-flux-adjusted models overestimate the magnitude of January-minus-July temperature differences by ~5 K due to an overestimate of summer (January) near-surface temperature. This error is common to five of the eight non-flux-adjusted models.
      • Over Northern Hemisphere mid-latitude land areas, zonal mean differences between July and January temperatures simulated by the non-flux-adjusted models show a greater spread (positive and negative) about observed values than results from the flux-adjusted models.
      • The ability of coupled GCMs to simulate a reasonable seasonal cycle is a necessary condition for confidence in their prediction of long-term climatic changes (such as global warming), but it is not a sufficient condition unless the seasonal cycle and long-term changes involve similar climatic processes.
  • Coupled climate models do not simulate with reasonable accuracy clouds and some related hydrological processes (in particular those involving upper tropospheric humidity). Problems in the simulation of clouds and upper tropospheric humidity, remain worrisome because the associated processes account for most of the uncertainty in climate model simulations of anthropogenic change.
The precise magnitude of future changes in climate is still uncertain;[37] for the end of the 21st century (2071 to 2100), for SRES scenario A2, the change of global average SAT change from AOGCMs compared with 1961 to 1990 is +3.0 °C (4.8 °F) and the range is +1.3 to +4.5 °C (+2 to +7.2 °F).

In the IPCC's Fifth Assessment Report, it was stated that there was "...very high confidence that models reproduce the general features of the global-scale annual mean surface temperature increase over the historical period." However, the report also observed that the rate of warming over the period 1998-2012 was lower than that predicted by 111 out of 114 Coupled Model Intercomparison Project climate models.[38]

Relation to weather forecasting

The global climate models used for climate projections are very similar in structure to (and often share computer code with) numerical models for weather prediction but are nonetheless logically distinct.

Most weather forecasting is done on the basis of interpreting the output of numerical model results. Since forecasts are short—typically a few days or a week—such models do not usually contain an ocean model but rely on imposed SSTs. They also require accurate initial conditions to begin the forecast—typically these are taken from the output of a previous forecast, with observations blended in. Because the results are needed quickly the predictions must be run in a few hours; but because they only need to cover a week of real time these predictions can be run at higher resolution than in climate mode. Currently the ECMWF runs at 40 km (25 mi) resolution[39] as opposed to the 100-to-200 km (62-to-124 mi) scale used by typical climate models. Often nested models are run forced by the global models for boundary conditions, to achieve higher local resolution: for example, the Met Office runs a mesoscale model with an 11 km (6.8 mi) resolution[40] covering the UK, and various agencies in the U.S. also run nested models such as the NGM and NAM models. Like most global numerical weather prediction models such as the GFS, global climate models are often spectral models[41] instead of grid models. Spectral models are often used for global models because some computations in modeling can be performed faster thus reducing the time needed to run the model simulation.

Computations involved

Climate models use quantitative methods to simulate the interactions of the atmosphere, oceans, land surface, and ice. They are used for a variety of purposes from study of the dynamics of the climate system to projections of future climate.

All climate models take account of incoming energy as short wave electromagnetic radiation, chiefly visible and short-wave (near) infrared, as well as outgoing energy as long wave (far) infrared electromagnetic radiation from the earth. Any imbalance results in a change in temperature.

The most talked-about models of recent years have been those relating temperature to emissions of carbon dioxide (see greenhouse gas). These models project an upward trend in the surface temperature record, as well as a more rapid increase in temperature at higher altitudes.[42]

Three (or more properly, four since time is also considered) dimensional GCM's discretise the equations for fluid motion and energy transfer and integrate these over time. They also contain parametrisations for processes—such as convection—that occur on scales too small to be resolved directly.

Atmospheric GCMs (AGCMs) model the atmosphere and impose sea surface temperatures as boundary conditions. Coupled atmosphere-ocean GCMs (AOGCMs, e.g. HadCM3, EdGCM, GFDL CM2.X, ARPEGE-Climat[43]) combine the two models.

Models can range from relatively simple to quite complex:
  • A simple radiant heat transfer model that treats the earth as a single point and averages outgoing energy
  • this can be expanded vertically (radiative-convective models), or horizontally
  • finally, (coupled) atmosphere–ocean–sea ice global climate models discretise and solve the full equations for mass and energy transfer and radiant exchange.
This is not a full list; for example "box models" can be written to treat flows across and within ocean basins. Furthermore, other types of modelling can be interlinked, such as land use, allowing researchers to predict the interaction between climate and ecosystems.

Other climate models

Earth-system models of intermediate complexity (EMICs)

Depending on the nature of questions asked and the pertinent time scales, there are, on the one extreme, conceptual, more inductive models, and, on the other extreme, general circulation models operating at the highest spatial and temporal resolution currently feasible. Models of intermediate complexity bridge the gap.
One example is the Climber-3 model. Its atmosphere is a 2.5-dimensional statistical-dynamical model with 7.5° × 22.5° resolution and time step of 1/2 a day; the ocean is MOM-3 (Modular Ocean Model) with a 3.75° × 3.75° grid and 24 vertical levels.

Radiative-convective models (RCM)

One-dimensional, radiative-convective models were used to verify basic climate assumptions in the '80s and '90s.[44]

Climate modelers

A climate modeler is a person who designs, develops, implements, tests, maintains or exploits climate models. There are three major types of institutions where a climate modeller may be found:
Posted by at

Navier–Stokes equations

Navier–Stokes equations

From Wikipedia, the free encyclopedia
In physics, the Navier–Stokes equations [navˈjeː stəʊks], named after Claude-Louis Navier and George Gabriel Stokes, describe the motion of fluid substances. These equations arise from applying Newton's second law to fluid motion, together with the assumption that the stress in the fluid is the sum of a diffusing viscous term (proportional to the gradient of velocity) and a pressure term - hence describing viscous flow.

The equations are useful because they describe the physics of many things of scientific and engineering interest. They may be used to model the weather, ocean currents, water flow in a pipe and air flow around a wing. The Navier–Stokes equations in their full and simplified forms help with the design of aircraft and cars, the study of blood flow, the design of power stations, the analysis of pollution, and many other things.
Coupled with Maxwell's equations they can be used to model and study magnetohydrodynamics.

The Navier–Stokes equations are also of great interest in a purely mathematical sense. Somewhat surprisingly, given their wide range of practical uses, it has not yet been proven that in three dimensions solutions always exist (existence), or that if they do exist, then they do not contain any singularity (they are smooth). These are called the Navier–Stokes existence and smoothness problems. The Clay Mathematics Institute has called this one of the seven most important open problems in mathematics and has offered a US$1,000,000 prize for a solution or a counter-example.[1]

Velocity field

The solution of the Navier–Stokes equations is a velocity not position. It is called a velocity field or flow field, which is a description of the velocity of the fluid at a given point in space and time. Once the velocity field is solved for, other quantities of interest, such as flow rate or drag force, may be found. This is different from what one normally sees in classical mechanics, where solutions are typically trajectories of position of a particle or deflection of a continuum. Studying velocity instead of position makes more sense for a fluid; however for visualization purposes one can compute various trajectories.

Properties

Nonlinearity

The Navier–Stokes equations are nonlinear partial differential equations in almost every real situation.[2][3] In some cases, such as one-dimensional flow and Stokes flow (or creeping flow), the equations can be simplified to linear equations. The nonlinearity makes most problems difficult or impossible to solve and is the main contributor to the turbulence that the equations model.

The nonlinearity is due to convective acceleration, which is an acceleration associated with the change in velocity over position. Hence, any convective flow, whether turbulent or not, will involve nonlinearity. An example of convective but laminar (nonturbulent) flow would be the passage of a viscous fluid (for example, oil) through a small converging nozzle. Such flows, whether exactly solvable or not, can often be thoroughly studied and understood.[4]

Turbulence

Turbulence is the time-dependent chaotic behavior seen in many fluid flows. It is generally believed that it is due to the inertia of the fluid as a whole: the culmination of time-dependent and convective acceleration; hence flows where inertial effects are small tend to be laminar (the Reynolds number quantifies how much the flow is affected by inertia). It is believed, though not known with certainty, that the Navier–Stokes equations describe turbulence properly.[5]

The numerical solution of the Navier–Stokes equations for turbulent flow is extremely difficult, and due to the significantly different mixing-length scales that are involved in turbulent flow, the stable solution of this requires such a fine mesh resolution that the computational time becomes significantly infeasible for calculation or direct numerical simulation. Attempts to solve turbulent flow using a laminar solver typically result in a time-unsteady solution, which fails to converge appropriately. To counter this, time-averaged equations such as the Reynolds-averaged Navier–Stokes equations (RANS), supplemented with turbulence models, are used in practical computational fluid dynamics (CFD) applications when modeling turbulent flows. Some models include the Spalart-Allmaras, k-ω (k-omega), k-ε (k-epsilon), and SST models, which add a variety of additional equations to bring closure to the RANS equations. Large eddy simulation (LES) can also be used to solve these equations numerically. This approach is computationally more expensive-in time and in computer memory-than RANS, but produces better results because it explicitly resolves the larger turbulent scales.

Applicability

Together with supplemental equations (for example, conservation of mass) and well formulated boundary conditions, the Navier–Stokes equations seem to model fluid motion accurately; even turbulent flows seem (on average) to agree with real world observations.
The Navier–Stokes equations assume that the fluid being studied is a continuum (it is infinitely divisible and not composed of particles such as atoms or molecules), and is not moving at relativistic velocities. At very small scales or under extreme conditions, real fluids made out of discrete molecules will produce results different from the continuous fluids modeled by the Navier–Stokes equations. Depending on the Knudsen number of the problem, statistical mechanics or possibly even molecular dynamics may be a more appropriate approach.

Another limitation is simply the complicated nature of the equations. Time-tested formulations exist for common fluid families, but the application of the Navier–Stokes equations to less common families tends to result in very complicated formulations and often to open research problems. For this reason, these equations are usually rewritten for Newtonian fluids where the viscosity model is linear; truly general models for the flow of other kinds of fluids (such as blood) do not, as of 2012, exist .[citation needed]

Derivation and description

The derivation of the Navier–Stokes equations begins with an application of Newton's second law: conservation of momentum (often alongside mass and energy conservation) being written for an arbitrary portion of the fluid. In an inertial frame of reference, the general form of the equations of fluid motion is:[6]
Navier–Stokes equations (general) \rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} \right) = -\nabla p + \nabla \cdot\boldsymbol{\mathsf{T}} + \mathbf{f}
where
\mathbf{v} is the flow velocity,
\rho is the fluid density,
p is the pressure,
\boldsymbol{\mathsf{T}} is the (deviatoric) component of the total stress tensor, which has order two,
\mathbf{f} represents body forces (per unit volume) acting on the fluid,
\nabla is the del operator.
This is a statement of the conservation of momentum in a fluid and it is an application of Newton's second law to a continuum; in fact this equation is applicable to any non-relativistic continuum and is known as the Cauchy momentum equation.

This equation is often written using the material derivative Dv/Dt, making it more apparent that this is a statement of Newton's second law:
\rho \frac{D \mathbf{v}}{D t} = -\nabla p + \nabla \cdot\boldsymbol{\mathsf{T}} + \mathbf{f}.
The left side of the equation describes acceleration, and may be composed of time-dependent or convective effects (also the effects of non-inertial coordinates if present). The right side of the equation is in effect a summation of body forces (such as gravity) and divergence of stress (pressure and shear stress).

Convective acceleration


An example of convection. Though the flow may be steady (time-independent), the fluid decelerates as it moves down the diverging duct (assuming incompressible or subsonic compressible flow), hence there is an acceleration happening over position.

A significant feature of the Navier–Stokes equations is the presence of convective acceleration: the effect of time-independent acceleration of a fluid with respect to space. While individual fluid particles indeed experience time dependent-acceleration, the convective acceleration of the flow field is a spatial effect, one example being fluid speeding up in a nozzle.

Regardless of what kind of fluid is being dealt with, convective acceleration is a nonlinear effect. Convective acceleration is present in most flows (exceptions include one-dimensional incompressible flow), but its dynamic effect is disregarded in creeping flow (also called Stokes flow). Convective acceleration is represented by the nonlinear quantity:
\mathbf{v} \cdot \nabla \mathbf{v}
which may be interpreted either as (\mathbf{v}\cdot\nabla)\,\mathbf{v} or as \mathbf{v}\cdot(\nabla\mathbf{v}), with \nabla \mathbf{v} the tensor derivative of the velocity vector \mathbf{v}. Both interpretations give the same result, independent of the coordinate system — provided \nabla is interpreted as the covariant derivative.[7]

Interpretation as (v·∇)v

The convection term is often written as
(\mathbf{v} \cdot \nabla) \mathbf{v}
where the advection operator \mathbf{v} \cdot \nabla is used. Usually this representation is preferred as it is simpler than the one in terms of the tensor derivative \nabla \mathbf{v}.[7]

Interpretation as v·(∇v)

Here \nabla \mathbf{v} is the tensor derivative of the velocity vector, equal in Cartesian coordinates to the component-by-component gradient. Note that the gradient of a vector is being defined as \left[\nabla \mathbf{v}\right]_{mi}=\partial_m v_i, so that \left[\mathbf{v}\cdot\left(\nabla \mathbf{v}\right)\right]_i=v_m \partial_m v_i=\left[(\mathbf{v}\cdot\nabla)\mathbf{v}\right]_i.

In irrotational flow

The convection term may, by a vector calculus identity, be expressed without a tensor derivative:[8][9]
\mathbf{v} \cdot \nabla \mathbf{v} = \nabla \left( \frac{\|\mathbf{v}\|^2}{2} \right) + \left( \nabla \times \mathbf{v} \right) \times \mathbf{v}.
The form has use in irrotational flow, where the curl of the velocity (called vorticity) \omega=\nabla \times \mathbf{v} is equal to zero. Therefore, this reduces to only
\mathbf{v} \cdot \nabla \mathbf{v} = \nabla \left( \frac{\|\mathbf{v}\|^2}{2} \right).

Stresses

The effect of stress in the fluid is represented by the \nabla p and \nabla \cdot\boldsymbol{\mathsf{T}} terms; these are gradients of surface forces, analogous to stresses in a solid. Here \nabla p is called the pressure gradient and arises from the isotropic part of the Cauchy stress tensor, which has order two. This part is given by normal stresses that turn up in almost all situations, dynamic or not. The anisotropic part of the stress tensor gives rise to \nabla \cdot\boldsymbol{\mathsf{T}}, which conventionally describes viscous forces; for incompressible flow, this is only a shear effect. Thus, \boldsymbol{\mathsf{T}} is the deviatoric stress tensor, and the stress tensor is equal to:[10]
\sigma = -p\boldsymbol{\mathsf{I}} + \boldsymbol{\mathsf{T}}
where \boldsymbol{\mathsf{I}} is the 3×3 identity matrix. In the Navier–Stokes equations, only the gradient of pressure matters, not the pressure itself. The effect of the pressure gradient on the flow is to accelerate the fluid in the direction from high pressure to low pressure.

The stress terms p and \boldsymbol{\mathsf{T}} are yet unknown, so the general form of the equations of motion is not usable to solve problems. Besides the equations of motion—Newton's second law—a force model is needed relating the stresses to the fluid motion.[11] For this reason, assumptions based on natural observations are often applied to specify the stresses in terms of the other flow variables, such as velocity and density.
The Navier–Stokes equations result from the following assumptions on the deviatoric stress tensor \boldsymbol{\mathsf{T}}:[12]
  • the deviatoric stress vanishes for a fluid at rest, and, by Galilean invariance, also does not depend directly on the flow velocity, but only on spatial derivatives of the flow velocity
  • the deviatoric stress is expressed as the double-dot product of the tensor gradient \nabla\mathbf{v} of the flow velocity with, a fourth-order viscosity tensor \boldsymbol{\mathsf{A}}, i.e., \boldsymbol{\mathsf{T}} = \boldsymbol{\mathsf{A}} \left( \nabla\mathbf{v} \right)
  • the fluid is assumed to be isotropic, as with gases and simple liquids, and consequently \boldsymbol{\mathsf{A}} is an isotropic tensor; furthermore, since the deviatoric stress tensor is symmetric, it can be expressed in terms of two scalar dynamic viscosities μ and μ": \boldsymbol{\mathsf{T}} = 2 \mu \boldsymbol{\mathsf{E}} + \mu'' \Delta \boldsymbol{\mathsf{I}}, where \boldsymbol{\mathsf{E}}=\tfrac12 \left( \nabla\mathbf{v} \right) + \tfrac12 \left( \nabla\mathbf{v} \right)^\text{T} is the rate-of-strain tensor and \Delta = \nabla\cdot\mathbf{v} is the rate of expansion of the flow
  • the deviatoric stress tensor has zero trace, so for a three-dimensional flow 2μ + 3μ" = 0
As a result, the deviatoric stress tensor has the following form:[12]
\boldsymbol{\mathsf{T}} = 2 \mu \left( \boldsymbol{\mathsf{E}} - \tfrac13 \Delta \boldsymbol{\mathsf{I}} \right)
with the quantity between brackets being the non-isotropic part of the rate-of-strain tensor \boldsymbol{\mathsf{E}}. The dynamic viscosity μ need not be constant – in general, it depends on conditions like temperature and pressure, and, in turbulence modeling, eddy viscosity approximates the average deviatoric stress.
The pressure p is modeled by an equation of state.[13] For the special case of an incompressible flow, the pressure constrains the flow so that the volume of fluid elements is constant: isochoric flow resulting in a solenoidal velocity field with \nabla\cdot\mathbf{v}=0.[14]

Other forces

The vector field \mathbf{f} represents body forces. Typically, these consist of only gravity forces, but may include others, such as electromagnetic forces. In non-inertial coordinate frames, other "forces" associated with rotating coordinates may arise.

Often, these forces may be represented as the gradient of some scalar quantity F\,, with \mathbf{f} = \nabla F , in which case they are called conservative forces. Gravity in the z direction, for example, is the gradient of -\rho g z. Because pressure from such gravitation arises only as a gradient, we may include it in the pressure term as a body force p_m = p - F. The pressure and force terms on the right-hand side of the Navier–Stokes equation become
-\nabla p + \mathbf{f} = -\nabla p + \nabla F = -\nabla \left( p - F \right) = -\nabla p_m.

Other equations

The Navier–Stokes equations are strictly a statement of the conservation of momentum. To fully describe fluid flow, more information is needed, how much depending on the assumptions made. This additional information may include boundary data (no-slip, capillary surface, etc.), conservation of mass, conservation of energy, and/or an equation of state.

Continuity equation

Regardless of the flow assumptions, a statement of the conservation of mass is generally necessary. This is achieved through the mass continuity equation, given in its most general form as:
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0
or, using the substantive derivative:
\frac{D\rho}{Dt} + \rho (\nabla \cdot \mathbf{v}) = 0.

Incompressible flow of Newtonian fluids

The equations are simplified in the case of incompressible flow of a Newtonian fluid. Incompressibility rules out sound propagation or shock waves, so this simplification is not useful if these phenomena are of interest. The incompressible flow assumption typically holds well even with a "compressible" fluid — such as air at room temperature — at low Mach numbers up to about Mach 0.3. With incompressible flow and constant viscosity, the Navier–Stokes equations read[15]
Navier–Stokes equations (Incompressible flow) \rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \mathbf{f}.
in tensor notation:
Navier–Stokes equations (Incompressible flow) \rho\left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j} \right)=-\frac{\partial p}{\partial x_i}+\mu\frac{\partial^2 u_i}{\partial x_j\partial x_j}+f_i
Here f represents "other" body forces (per unit volume), such as gravity or centrifugal force. The shear stress term \nabla \cdot \boldsymbol{\mathsf{T}} becomes \mu \nabla^2 \mathbf{v}, where \nabla^2 is the vector Laplacian.[16]

It's well worth observing the meaning of each term (compare to the Cauchy momentum equation):
\overbrace{\rho \Big( \underbrace{\frac{\partial \mathbf{v}}{\partial t}}_{ \begin{smallmatrix} \text{Unsteady}\\ \text{acceleration} \end{smallmatrix}} + \underbrace{\mathbf{v} \cdot \nabla \mathbf{v}}_{ \begin{smallmatrix} \text{Convective} \\ \text{acceleration} \end{smallmatrix}}\Big)}^{\text{Inertia (per volume)}} = \overbrace{\underbrace{-\nabla p}_{ \begin{smallmatrix} \text{Pressure} \\ \text{gradient} \end{smallmatrix}} + \underbrace{\mu \nabla^2 \mathbf{v}}_{\text{Viscosity}}}^{\text{Divergence of stress}} + \underbrace{\mathbf{f}.}_{ \begin{smallmatrix} \text{Other} \\ \text{body} \\ \text{forces} \end{smallmatrix}}
Only the convective terms are nonlinear for incompressible Newtonian flow. Convective acceleration is caused by a (possibly steady) change in velocity over position, for example, by the speeding up of fluid entering a converging nozzle. Though individual fluid particles are accelerated and thus under unsteady motion, the flow field, a velocity distribution, will not necessarily be time-dependent.

Another important observation is that viscosity is represented by the vector Laplacian of the velocity field, interpreted here as the difference between the velocity at a point and the mean velocity in a small surrounding volume. This implies that – for a Newtonian fluid – viscosity operates as a diffusion of momentum, in much the same way as the diffusion of heat in the heat equation, also expressed with the Laplacian.

If temperature effects are also neglected, the only "other" equation (apart from initial/boundary conditions) needed is the mass continuity equation. Under the assumption of incompressibility, the density of a fluid parcel is constant, and when using the substantive derivative it follows easily that the continuity equation simplifies to:
\nabla \cdot \mathbf{v} = 0.
This is more specifically a statement of the conservation of volume (see divergence and isochoric process).
These equations are commonly used in 3 coordinates systems: Cartesian, cylindrical, and spherical. While the Cartesian equations seem to follow directly from the vector equation above, the vector form of the Navier–Stokes equation involves some tensor calculus which means that writing it in other coordinate systems is not as simple as doing so for scalar equations (such as the heat equation).

Cartesian coordinates

With the velocity vector expanded as \mathbf{v} = (u,v,w), we may write the vector equation explicitly,
\begin{align} \rho \left(\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z}\right) &= -\frac{\partial p}{\partial x} + \mu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right) + \rho g_x \\ \rho \left(\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}+ w \frac{\partial v}{\partial z}\right) &= -\frac{\partial p}{\partial y} + \mu \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial z^2}\right) + \rho g_y \\ \rho \left(\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y}+ w \frac{\partial w}{\partial z}\right) &= -\frac{\partial p}{\partial z} + \mu \left(\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2}\right) + \rho g_z. \end{align}
Note that gravity has been accounted for as a body force, and the values of gx, gy, gz will depend on the orientation of gravity with respect to the chosen set of coordinates.

The continuity equation reads:
{\partial \rho \over \partial t} + {\partial (\rho u ) \over \partial x} + {\partial (\rho v) \over \partial y} + {\partial (\rho w) \over \partial z} = 0.
When the flow is incompressible, \rho does not change for any fluid parcel, and its material derivative vanishes:
{D\rho \over Dt} = 0. The continuity equation is reduced to:
{\partial u \over \partial x} + {\partial v \over \partial y} + {\partial w \over \partial z} = 0.
The velocity components (the dependent variables to be solved for) are typically named u, v, w. This system of four equations comprises the most commonly used and studied form. Though comparatively more compact than other representations, this is still a nonlinear system of partial differential equations for which solutions are difficult to obtain.

Cylindrical coordinates

A change of variables on the Cartesian equations will yield[15] the following momentum equations for r, \phi, and z:
\begin{align} r:\ &\rho \left(\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_r}{\partial \phi} + u_z \frac{\partial u_r}{\partial z} - \frac{u_{\phi}^2}{r}\right) = {}\\ &-\frac{\partial p}{\partial r} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_r}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_r}{\partial \phi^2} + \frac{\partial^2 u_r}{\partial z^2} - \frac{u_r}{r^2} - \frac{2}{r^2}\frac{\partial u_\phi}{\partial \phi} \right] + \rho g_r \\ \phi:\ &\rho \left(\frac{\partial u_{\phi}}{\partial t} + u_r \frac{\partial u_{\phi}}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_{\phi}}{\partial \phi} + u_z \frac{\partial u_{\phi}}{\partial z} + \frac{u_r u_{\phi}}{r}\right) = {}\\ &-\frac{1}{r}\frac{\partial p}{\partial \phi} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_{\phi}}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_{\phi}}{\partial \phi^2} + \frac{\partial^2 u_{\phi}}{\partial z^2} + \frac{2}{r^2}\frac{\partial u_r}{\partial \phi} - \frac{u_{\phi}}{r^2}\right] + \rho g_{\phi} \\ z:\ &\rho \left(\frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_z}{\partial \phi} + u_z \frac{\partial u_z}{\partial z}\right) = {}\\ &-\frac{\partial p}{\partial z} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_z}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_z}{\partial \phi^2} + \frac{\partial^2 u_z}{\partial z^2}\right] + \rho g_z. \end{align}
The gravity components will generally not be constants, however for most applications either the coordinates are chosen so that the gravity components are constant or else it is assumed that gravity is counteracted by a pressure field (for example, flow in horizontal pipe is treated normally without gravity and without a vertical pressure gradient). The continuity equation is:
\frac{\partial\rho}{\partial t} + \frac{1}{r}\frac{\partial}{\partial r}\left(\rho r u_r\right) + \frac{1}{r}\frac{\partial (\rho u_\phi)}{\partial \phi} + \frac{\partial (\rho u_z)}{\partial z} = 0.
This cylindrical representation of the incompressible Navier–Stokes equations is the second most commonly seen (the first being Cartesian above). Cylindrical coordinates are chosen to take advantage of symmetry, so that a velocity component can disappear. A very common case is axisymmetric flow with the assumption of no tangential velocity (u_{\phi} = 0), and the remaining quantities are independent of \phi:
\begin{align} \rho \left(\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + u_z \frac{\partial u_r}{\partial z}\right) &= -\frac{\partial p}{\partial r} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_r}{\partial r}\right) + \frac{\partial^2 u_r}{\partial z^2} - \frac{u_r}{r^2}\right] + \rho g_r \\ \rho \left(\frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + u_z \frac{\partial u_z}{\partial z}\right) &= -\frac{\partial p}{\partial z} + \mu \left[\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_z}{\partial r}\right) + \frac{\partial^2 u_z}{\partial z^2}\right] + \rho g_z \\ \frac{1}{r}\frac{\partial}{\partial r}\left(r u_r\right) + \frac{\partial u_z}{\partial z} &= 0. \end{align}

Spherical coordinates

In spherical coordinates, the r, ϕ, and θ momentum equations are[15] (note the convention used: θ is polar angle, or colatitude,[17] 0 ≤ θ ≤ π):
\begin{align} r:\ &\rho \left(\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{u_{\phi}}{r \sin(\theta)} \frac{\partial u_r}{\partial \phi} + \frac{u_{\theta}}{r} \frac{\partial u_r}{\partial \theta} - \frac{u_{\phi}^2 + u_{\theta}^2}{r}\right) = -\frac{\partial p}{\partial r} + \rho g_r + \\ &\mu \left[\frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u_r}{\partial r}\right) + \frac{1}{r^2 \sin(\theta)^2} \frac{\partial^2 u_r}{\partial \phi^2} + \frac{1}{r^2 \sin(\theta)} \frac{\partial}{\partial \theta}\left(\sin(\theta) \frac{\partial u_r}{\partial \theta}\right) - 2\frac{u_r + \frac{\partial u_{\theta}}{\partial \theta} + u_{\theta} \cot(\theta)}{r^2} - \frac{2}{r^2 \sin(\theta)} \frac{\partial u_{\phi}}{\partial \phi} \right] \\ \phi:\ &\rho \left(\frac{\partial u_{\phi}}{\partial t} + u_r \frac{\partial u_{\phi}}{\partial r} + \frac{u_{\phi}}{r \sin(\theta)} \frac{\partial u_{\phi}}{\partial \phi} + \frac{u_{\theta}}{r} \frac{\partial u_{\phi}}{\partial \theta} + \frac{u_r u_{\phi} + u_{\phi} u_{\theta} \cot(\theta)}{r}\right) = -\frac{1}{r \sin(\theta)} \frac{\partial p}{\partial \phi} + \rho g_{\phi} + \\ &\mu \left[\frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u_{\phi}}{\partial r}\right) + \frac{1}{r^2 \sin(\theta)^2} \frac{\partial^2 u_{\phi}}{\partial \phi^2} + \frac{1}{r^2 \sin(\theta)} \frac{\partial}{\partial \theta}\left(\sin(\theta) \frac{\partial u_{\phi}}{\partial \theta}\right) + \frac{2 \sin(\theta) \frac{\partial u_r}{\partial \phi} + 2 \cos(\theta) \frac{\partial u_{\theta}}{\partial \phi} - u_{\phi}}{r^2 \sin(\theta)^2} \right] \\ \theta:\ &\rho \left(\frac{\partial u_{\theta}}{\partial t} + u_r \frac{\partial u_{\theta}}{\partial r} + \frac{u_{\phi}}{r \sin(\theta)} \frac{\partial u_{\theta}}{\partial \phi} + \frac{u_{\theta}}{r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{u_r u_{\theta} - u_{\phi}^2 \cot(\theta)}{r}\right) = -\frac{1}{r} \frac{\partial p}{\partial \theta} + \rho g_{\theta} + \\ &\mu \left[\frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u_{\theta}}{\partial r}\right) + \frac{1}{r^2 \sin(\theta)^2} \frac{\partial^2 u_{\theta}}{\partial \phi^2} + \frac{1}{r^2 \sin(\theta)} \frac{\partial}{\partial \theta}\left(\sin(\theta) \frac{\partial u_{\theta}}{\partial \theta}\right) + \frac{2}{r^2} \frac{\partial u_r}{\partial \theta} - \frac{u_{\theta} + 2 \cos(\theta) \frac{\partial u_{\phi}}{\partial \phi}}{r^2 \sin(\theta)^2} \right]. \end{align}
Mass continuity will read:
\frac{\partial \rho}{\partial t} + \frac{1}{r^2}\frac{\partial}{\partial r}\left(\rho r^2 u_r\right) + \frac{1}{r \sin(\theta)}\frac{\partial \rho u_\phi}{\partial \phi} + \frac{1}{r \sin(\theta)}\frac{\partial}{\partial \theta}\left(\sin(\theta) \rho u_\theta\right) = 0.
These equations could be (slightly) compacted by, for example, factoring 1/r^2 from the viscous terms. However, doing so would undesirably alter the structure of the Laplacian and other quantities.

Stream function formulation

Taking the curl of the Navier–Stokes equation results in the elimination of pressure. This is especially easy to see if 2D Cartesian flow is assumed (w = 0 and no dependence of anything on z), where the equations reduce to:
\begin{align} \rho \left(\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}\right) &= -\frac{\partial p}{\partial x} + \mu \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + \rho g_x \\ \rho \left(\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}\right) &= -\frac{\partial p}{\partial y} + \mu \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) + \rho g_y. \end{align}
Differentiating the first with respect to y, the second with respect to x and subtracting the resulting equations will eliminate pressure and any conservative force. Defining the stream function \psi through
u = \frac{\partial \psi}{\partial y}; \quad v = -\frac{\partial \psi}{\partial x}
results in mass continuity being unconditionally satisfied (given the stream function is continuous), and then incompressible Newtonian 2D momentum and mass conservation degrade into one equation:
\frac{\partial}{\partial t}\left(\nabla^2 \psi\right) + \frac{\partial \psi}{\partial y} \frac{\partial}{\partial x}\left(\nabla^2 \psi\right) - \frac{\partial \psi}{\partial x} \frac{\partial}{\partial y}\left(\nabla^2 \psi\right) = \nu \nabla^4 \psi
where \nabla^4 is the (2D) biharmonic operator and \nu is the kinematic viscosity, \nu = \frac{\mu}{\rho}. We can also express this compactly using the Jacobian determinant:
\frac{\partial}{\partial t}\left(\nabla^2 \psi\right) + \frac{\partial\left(\psi, \nabla^2\psi \right)}{\partial\left(y,x\right)} = \nu \nabla^4 \psi.
This single equation together with appropriate boundary conditions describes 2D fluid flow, taking only kinematic viscosity as a parameter. Note that the equation for creeping flow results when the left side is assumed zero.

In axisymmetric flow another stream function formulation, called the Stokes stream function, can be used to describe the velocity components of an incompressible flow with one scalar function.

Pressure-free velocity formulation

The incompressible Navier–Stokes equation is a differential algebraic equation, having the inconvenient feature that there is no explicit mechanism for advancing the pressure in time. Consequently, much effort has been expended to eliminate the pressure from all or part of the computational process. The stream function formulation above eliminates the pressure (in 2D) at the expense of introducing higher derivatives and elimination of the velocity, which is the primary variable of interest.

The incompressible Navier–Stokes equation is composite, the sum of two orthogonal equations,
\begin{align} \frac{\partial\mathbf{v}}{\partial t} &= \Pi^S\left(-\mathbf{v}\cdot\nabla\mathbf{v} + \nu\nabla^2\mathbf{v}\right) + \mathbf{f}^S \\ \rho^{-1}\nabla p &= \Pi^I\left(-\mathbf{v}\cdot\nabla\mathbf{v} + \nu\nabla^2\mathbf{v}\right) + \mathbf{f}^I \end{align}
where \Pi^S and \Pi^I are solenoidal and irrotational projection operators satisfying \Pi^S+\Pi^I=1 and \mathbf{f}^S and \mathbf{f}^I are the non-conservative and conservative parts of the body force. This result follows from the Helmholtz Theorem (also known as the fundamental theorem of vector calculus). The first equation is a pressureless governing equation for the velocity, while the second equation for the pressure is a functional of the velocity and is related to the pressure Poisson equation.

The explicit functional form of the projection operator in 3D is found from the Helmholtz Theorem
\Pi^S\,\mathbf{F}(\mathbf{r}) = \frac{1}{4\pi}\nabla\times\int \frac{\nabla^\prime\times\mathbf{F}(\mathbf{r}^\prime)}{|\mathbf{r}-\mathbf{r}^\prime|} d V^\prime, \quad \Pi^I = 1-\Pi^S
with a similar structure in 2D. Thus the governing equation is an integro-differential equation and not convenient for numerical computation.

An equivalent weak or variational form of the equation, proved to produce the same velocity solution as the Navier–Stokes equation,[18] is given by,
\left(\mathbf{w},\frac{\partial\mathbf{v}}{\partial t}\right) = -(\mathbf{w},\mathbf{v}\cdot\nabla\mathbf{v})-\nu(\nabla\mathbf{w}: \nabla\mathbf{v})+(\mathbf{w},\mathbf{f}^S)
for divergence-free test functions \mathbf{w} satisfying appropriate boundary conditions. Here, the projections are accomplished by the orthogonality of the solenoidal and irrotational function spaces. The discrete form of this is imminently suited to finite element computation of divergence-free flow, as we shall see in the next section. There we will be able to address the question, "How does one specify pressure-driven (Poiseuille) problems with a pressureless governing equation?"

The absence of pressure forces from the governing velocity equation demonstrates that the equation is not a dynamic one, but rather a kinematic equation where the divergence-free condition serves the role of a conservation law. This all would seem to refute the frequent statements that the incompressible pressure enforces the divergence-free condition.

Discrete velocity

With partitioning of the problem domain and defining basis functions on the partitioned domain, the discrete form of the governing equation is,
\left(\mathbf{w}_i, \frac{\partial\mathbf{v}_j}{\partial t}\right) = -(\mathbf{w}_i, \mathbf{v}\cdot\nabla\mathbf{v}_j) - \nu(\nabla\mathbf{w}_i: \nabla\mathbf{v}_j) + \left(\mathbf{w}_i, \mathbf{f}^S\right).
It is desirable to choose basis functions which reflect the essential feature of incompressible flow – the elements must be divergence-free. While the velocity is the variable of interest, the existence of the stream function or vector potential is necessary by the Helmholtz Theorem. Further, to determine fluid flow in the absence of a pressure gradient, one can specify the difference of stream function values across a 2D channel, or the line integral of the tangential component of the vector potential around the channel in 3D, the flow being given by Stokes' Theorem. Discussion will be restricted to 2D in the following.

We further restrict discussion to continuous Hermite finite elements which have at least first-derivative degrees-of-freedom. With this, one can draw a large number of candidate triangular and rectangular elements from the plate-bending literature. These elements have derivatives as components of the gradient. In 2D, the gradient and curl of a scalar are clearly orthogonal, given by the expressions,
\nabla\phi = \left[\frac{\partial \phi}{\partial x},\,\frac{\partial \phi}{\partial y}\right]^\mathrm{T}, \quad \nabla\times\phi = \left[\frac{\partial \phi}{\partial y},\,-\frac{\partial \phi}{\partial x}\right]^\mathrm{T}.
Adopting continuous plate-bending elements, interchanging the derivative degrees-of-freedom and changing the sign of the appropriate one gives many families of stream function elements.

Taking the curl of the scalar stream function elements gives divergence-free velocity elements.[19][20] The requirement that the stream function elements be continuous assures that the normal component of the velocity is continuous across element interfaces, all that is necessary for vanishing divergence on these interfaces.

Boundary conditions are simple to apply. The stream function is constant on no-flow surfaces, with no-slip velocity conditions on surfaces. Stream function differences across open channels determine the flow. No boundary conditions are necessary on open boundaries, though consistent values may be used with some problems. These are all Dirichlet conditions.

The algebraic equations to be solved are simple to set up, but of course are non-linear, requiring iteration of the linearized equations.

Similar considerations apply to three-dimensions, but extension from 2D is not immediate because of the vector nature of the potential, and there exists no simple relation between the gradient and the curl as was the case in 2D.

Pressure recovery

Recovering pressure from the velocity field is easy. The discrete weak equation for the pressure gradient is,
(\mathbf{g}_i, \nabla p) = -(\mathbf{g}_i, \mathbf{v}\cdot\nabla\mathbf{v}_j) - \nu(\nabla\mathbf{g}_i: \nabla\mathbf{v}_j) + (\mathbf{g}_i, \mathbf{f}^I)
where the test/weight functions are irrotational. Any conforming scalar finite element may be used. However, the pressure gradient field may also be of interest. In this case one can use scalar Hermite elements for the pressure. For the test/weight functions \mathbf{g}_i one would choose the irrotational vector elements obtained from the gradient of the pressure element.

Compressible flow of Newtonian fluids

There are some phenomena that are closely linked with fluid compressibility. One of the obvious examples is sound. Description of such phenomena requires more general presentation of the Navier–Stokes equation that takes into account fluid compressibility. If viscosity is assumed a constant, one additional term appears, as shown here:[21][22]
\rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \left(\zeta + \frac{\mu}{3}\right) \nabla (\nabla \cdot \mathbf{v}) + \mathbf{f},
where \zeta is the second viscosity.

Application to specific problems

The Navier–Stokes equations, even when written explicitly for specific fluids, are rather generic in nature and their proper application to specific problems can be very diverse. This is partly because there is an enormous variety of problems that may be modeled, ranging from as simple as the distribution of static pressure to as complicated as multiphase flow driven by surface tension.

Generally, application to specific problems begins with some flow assumptions and initial/boundary condition formulation, this may be followed by scale analysis to further simplify the problem.

Visualization of a) parallel flow and b) radial flow.

a) Assume steady, parallel, one dimensional, non-convective pressure-driven flow between parallel plates, the resulting scaled (dimensionless) boundary value problem is:
\frac{d^2 u}{d y^2} = -1; \quad u(0) = u(1) = 0.
The boundary condition is the no slip condition. This problem is easily solved for the flow field:
u(y) = \frac{y - y^2}{2}.
From this point onward more quantities of interest can be easily obtained, such as viscous drag force or net flow rate.

b) Difficulties may arise when the problem becomes slightly more complicated. A seemingly modest twist on the parallel flow above would be the radial flow between parallel plates; this involves convection and thus non-linearity. The velocity field may be represented by a function f(z) that must satisfy:
\frac{d^2 f}{d z^2} + R f^2 = -1; \quad f(-1) = f(1) = 0.
This ordinary differential equation is what is obtained when the Navier–Stokes equations are written and the flow assumptions applied (additionally, the pressure gradient is solved for). The nonlinear term makes this a very difficult problem to solve analytically (a lengthy implicit solution may be found which involves elliptic integrals and roots of cubic polynomials). Issues with the actual existence of solutions arise for R > 1.41 (approximately; this is not the square root of 2), the parameter R being the Reynolds number with appropriately chosen scales.[23] This is an example of flow assumptions losing their applicability, and an example of the difficulty in "high" Reynolds number flows.[23]

Exact solutions of the Navier–Stokes equations

Some exact solutions to the Navier–Stokes equations exist. Examples of degenerate cases — with the non-linear terms in the Navier–Stokes equations equal to zero — are Poiseuille flow, Couette flow and the oscillatory Stokes boundary layer. But also more interesting examples, solutions to the full non-linear equations, exist; for example the Taylor–Green vortex.[24][25][26] Note that the existence of these exact solutions does not imply they are stable: turbulence may develop at higher Reynolds numbers.

Some of the flow lines along a Hopf fibration.

A nice steady-state example with no singularities comes from considering the flow along the lines of a Hopf fibration. Let r be a constant radius to the inner coil. One set of solutions is given by:[28]
\begin{align} \rho(x, y, z) &= \frac{3B}{r^2 + x^2 + y^2 + z^2} \\ p(x, y, z) &= \frac{-A^2B}{(r^2 + x^2 + y^2 + z^2)^3} \\ \mathbf{v}(x, y, z) &= \frac{A}{(r^2 + x^2 + y^2 + z^2)^2}\begin{pmatrix} 2(-ry + xz) \\ 2(rx + yz) \\ r^2 - x^2 - y^2 + z^2 \end{pmatrix} \\ g &= 0 \\ \mu &= 0 \end{align}
for arbitrary constants A and B. This is a solution in a non-viscous gas (compressible fluid) whose density, velocities and pressure goes to zero far from the origin. (Note this is not a solution to the Clay Millennium problem because that refers to incompressible fluids where \rho is a constant, neither does it deal with the uniqueness of the Navier–Stokes equations with respect to any turbulence properties.) It is also worth pointing out that the components of the velocity vector are exactly those from the Pythagorean quadruple parametrization. Other choices of density and pressure are possible with the same velocity field:
Wyld diagrams

Wyld diagrams are bookkeeping graphs that correspond to the Navier–Stokes equations via a perturbation expansion of the fundamental continuum mechanics. Similar to the Feynman diagrams in quantum field theory, these diagrams are an extension of Keldysh's technique for nonequilibrium processes in fluid dynamics. In other words, these diagrams assign graphs to the (often) turbulent phenomena in turbulent fluids by allowing correlated and interacting fluid particles to obey stochastic processes associated to pseudo-random functions in probability distributions.[29]

Navier–Stokes equations use in games

The Navier–Stokes equations are used extensively in video games in order to model a wide variety of natural phenomena. Simulations of small-scale gaseous fluids, such as fire and smoke are often based on the seminal paper "Real-Time Fluid Dynamics for Games"[30] by Jos Stam, which elaborates one of the methods proposed in Stam's earlier, more famous paper "Stable Fluids"[31] from 1999. Stam proposes stable fluid simulation using a Navier–Stokes solution method from 1968, coupled with an unconditionally stable semi-Lagrangian advection scheme, as first proposed in 1992.

More recent implementations based upon this work run on the GPU as opposed to the CPU and achieve a much higher degree of performance.[32][33] Many improvements have been proposed to Stam's original work, which suffers inherently from high numerical dissipation in both velocity and mass.

An introduction to interactive fluid simulation can be found in the 2007 ACM SIGGRAPH course, Fluid Simulation for Computer Animation.[34]
Posted by at
Subscribe to: Posts (Atom)