text:
Solar System, technical/Venus

It is known as a common fact, Venus is the second planet in terms of proximity towards the Sun. Venus is also known as the Morning or the Evening star due to the fact that it is the brightest right before sunrise or a little after sunset. Venus is also known as the Earth's 'sister planet' because of its diameter, size, and bulk composition as it is nearly same as the earth's. Venus shines very well in night sky as its atmosphere is very thick and it reflects sunlight very brightly. Venus is also the hottest planet in our solar system. This is also because of its atmosphere. Its atmosphere mainly consists of sulphuric acid and carbon dioxide (CO2) that traps heat of the Sun and do not let it go back. Venus has no moons. It has a surface temperature of about 450℃~470℃. Venus is the only planet named after a female figure. Statistics and other important informations about Venus are given below.

Some information about the planet

 * Day (time taken to spin round its orbit): 243 Earth days
 * Year (time taken to cover its orbit): 224.7 Earth days
 * Diameter (at equator) : 12,102 km
 * Temperature: 450℃~470℃.
 * Distance from Sun: 107 million km~109 million km
 * Spacecraft from earth: Magellan orbiter (which reached Venus in 1980s)

Orbital characteristics

 * Aphelion: 108,942,109 km (0.72823128 AU)
 * Perihelion: 107,476,259 km (0.71843270 AU)
 * Semi-major axis: 108,208,930 km (0.723332 AU)
 * Eccentricity: 	0.0068
 * Orbital period: 224.70069 day (0.6151970 yr)
 * Synodic period: 583.92 days
 * Avg. orbital speed: 35.02 km/s
 * Inclination: 3.39471° (3.86° to Sun's equator)
 * Longitude of ascending node: 76.67069°
 * Argument of perihelion: 54.85229°
 * Satellites: None

Physical characteristics
Kelvin: 735K Celsius: 461.85℃
 * Mean radius: 6051.8 ± 1.0 km (0.9499 Earths)
 * Flattening: < 0.0002
 * Surface area: 	4.60×108 km² (0.902 Earths)
 * Volume: 9.38×1011 km³ (0.857 Earths)
 * Mass: 4.8685×1024 kg (0.815 Earths)
 * Mean density: 5.204 g/cm³
 * Equatorial surface gravity: 8.87 m/s2 (0.904 g)
 * Escape velocity: 10.46 km/s
 * Sidereal rotation period: −243.0185 day
 * Rotation velocity at equator: 6.52 km/h
 * Axial tilt: 177.36°
 * Right ascension of North pole: 18 hours, 11 minutes, 2 seconds (272.76°)
 * Declination of North pole: 67.16°
 * Albedo: 0.65
 * Surface temp.:
 * Apparent magnitude: up to -4.6
 * Angular diameter: 9.7" — 66.0"
 * Adjectives: Venusian or (rarely) Cytherean

Atmosphere

 * Surface pressure: 9.3 MPa

 * Composition:
 * 1) ≈ 96.5% Carbon dioxide
 * 2) ≈ 3.5% Nitrogen
 * 3) 0.015% Sulphur dioxide
 * 4) 0.007% Argon
 * 5) 0.002% Water vapor
 * 6) 0.0017% Carbon monoxide
 * 7) 0.0012% Helium
 * 8) 0.0007% Neon
 * 9) Trace Carbonyl sulfide
 * 10) Trace Hydrogen chloride
 * 11) Trace Hydrogen fluoride

Meteor astronomy
In visual astronomy "[a]lmost no variation or detail can be seen in the clouds. ... The surface is obscured by a thick blanket of clouds. ... Venus is shrouded by an opaque layer of highly reflective clouds of sulfuric acid, preventing its surface from being seen from space in visible light. ... [It has] thick clouds of sulfur dioxide ... [There are] lower and middle cloud layers ... [The] thick clouds consisting mainly of sulfur dioxide and sulfuric acid droplets. These clouds reflect and scatter about 90% of the sunlight that falls on them back into space, and prevent visual observation of the Venusian surface. The permanent cloud cover means that although Venus is closer than Earth to the Sun, the Venusian surface is not as well lit."

"Strong 300 km/h winds at the cloud tops circle the planet about every four to five earth days. Venusian winds move at up to 60 times the speed of the planet's rotation, while Earth's fastest winds are only 10% to 20% rotation speed. "

"Average cloud-top wind speeds on Venus rose 33 percent between 2006 and 2012, jumping from 186 mph (300 km/h) to 249 mph (400 km/h), observations by Europe's Venus Express orbiter show."

"This is an enormous increase in the already high wind speeds known in the atmosphere, ... Such a large variation has never before been observed on Venus, and we do not yet understand why this occurred."

The wind speeds are determined "by studying images captured by Venus Express between [50° N and S latitude, and tracking] the movements of tens of thousands of feature in the cloud tops some [70 km] above the planet's surface."

"Our analysis of cloud motions at low latitudes in the southern hemisphere showed that over the six years of study the velocity of the winds changed by up 70 km/h over a time scale of 255 Earth days — slightly longer than a year on Venus".

"Sometimes clouds took 3.9 days to zip all the way around Venus, for example, while on other occasions the journey required 5.3 days."

"Although there is clear evidence that the average global wind speeds have increased, further investigations are needed in order to explain what drives the atmospheric circulation patterns that are responsible, and to explain the changes seen in localized areas on shorter timescales ... The atmospheric super-rotation of Venus is one of the great unexplained mysteries of the solar system ... These results add more mystery to it, as Venus Express continues to surprise us with its ongoing observations of this dynamic, changing planet."

"Venus has a super-rotating atmosphere that whips around the planet once every four Earth days; Venus itself takes 243 Earth days to complete one rotation."

X-ray astronomy
The first ever X-ray image of Venus is shown at right. The "half crescent is due to the relative orientation of the Sun, Earth and Venus. The X-rays from Venus are produced by fluorescent radiation from oxygen and other atoms in the atmosphere between 120 and 140 kilometers above the surface of the planet. In contrast, the optical light from Venus is caused by the reflection from clouds 50 to 70 kilometers above the surface. Solar X-rays bombard the atmosphere of Venus, knock electrons out of the inner parts of atoms, and excite the atoms to a higher energy level. The atoms almost immediately return to their lower energy state with the emission of a fluorescent X-ray. A similar process involving ultraviolet light produces the visible light from fluorescent lamps."

Ultraviolet astronomy
When imaged in the ultraviolet (article top right), Venus appears like a gas dwarf object rather than a rocky object.

The image on the right has been re-processed through the clear, blue, and UV filters of Mariner 10 from the image taken of Venus by Mariner 10 on May 5, 1974, to show greater detail.

Visual astronomy

When imaged in visible light (left) Venus appears like a gas dwarf rather than a rocky body. The same image result occurs when it is viewed in the ultraviolet (right at page top).

Violets
Violet photographs of the planet Venus taken in 1927 “recorded two nebulous bright streaks, or bands, running ... approximately at right angles to the terminator” that may be from the upper atmosphere.

Green astronomy
"Venus at times shows [the oxygen] green line emission with an intensity equal to terrestrial values [Slanger et al., 2001]. Furthermore, the intensity is quite variable, as is true for the much stronger O2( a-X) 1.27 μ emission."

"In 1999, observations of the Venus nightglow with the Keck I telescope showed that the 5577 Å oxygen green line was a significant feature, comparable in intensity to the terrestrial green line. Subsequent measurements have been carried out at the Apache Point Observatory (APO) and again at Keck I, confirming the presence of the line with substantially varying intensity."

"Ground-based studies suggest that the [557.7 nm oxygen green line] emission is correlated with the solar cycle."

Infrareds
"The Herzberg II system of O2 has been a known feature of Venus' nightglow since the Venera 9 and 10 orbiters detected its c(0)-X(v″) progression more than 3 decades ago."

"Spectroscopic observations of the differential Doppler shift in a CO2 absorption line on Venus show that the upper atmospheric wind near the equator appears to have both a retrograde motion of about -85 ± 10 m s-1 ... and ... a periodically varying component, with an amplitude of about 40 ± 14 m s-1 and a period of 4.3 ± 0.2 days."

At right is a false-color near-infrared image of the lower-level clouds on the night side of Venus, obtained by the Near Infrared Mapping Spectrometer aboard the Galileo spacecraft as it approached the planet's night side on February 10, 1990.

"Bright slivers of sunlit high clouds are visible above and below the dark, glowing hemisphere. The spacecraft is about 100,000 kilometers (60,000 miles) above the planet. An infrared wavelength of 2.3 microns (about three times the longest wavelength visible to the human eye) was used. The map shows the turbulent, cloudy middle atmosphere some 50-55 kilometers (30- 33 miles) above the surface, 10-16 kilometers or 6-10 miles below the visible cloudtops. The red color represents the radiant heat from the lower atmosphere (about 400 degrees Fahrenheit) shining through the sulfuric acid clouds, which appear as much as 10 times darker than the bright gaps between clouds. This cloud layer is at about -30 degrees Fahrenheit, at a pressure about 1/2 Earth's surface atmospheric pressure. Near the equator, the clouds appear fluffy and blocky; farther north, they are stretched out into East-West filaments by winds estimated at more than 150 mph, while the poles are capped by thick clouds at this altitude."

Radio astronomy
"The first un-ambiguous detection of Venus was made by [the] Jet Propulsion Laboratory (JPL) on 10 March 1961. A correct measurement of the AU soon followed."

"The advantages of radar in planetary astronomy result from (1) the observer's control of all the attributes of the coherent signal used to illuminate the target, especially the wave form's time/frequency modulation and polarization; (2) the ability of radar to resolve objects spatially via measurements of the distribution of echo power in time delay and Doppler frequency; (3) the pronounced degree to which delay-Doppler measurements constrain orbits and spin vectors; and (4) centimeter-to-meter wavelengths, which easily penetrate optically opaque planetary clouds and cometary comae, permit investigation of near-surface macrostructure and bulk density, and are sensitive to high concentrations of metal or, in certain situations, ice."

When viewed using radio astronomy, the resulting radar image, at left, shows that just beneath the cloud layers is a rocky object.

Gaseous objects
Venus has been detected as a gaseous object using X-ray through red astronomy.

Atmospheric astronomy
"During a rare period of very low density solar outflow, the ionosphere of Venus was observed to become elongated downstream, rather like a long-tailed comet. ... Under normal conditions, the solar wind has a density of 5 - 10 particles per cubic cm at Earth's orbit, but occasionally the solar wind almost disappears, as happened in May 1999. ... A rare opportunity to examine what happens when a tenuous solar wind arrives at Venus came 3 - 4 August 2010, following a series of large coronal mass ejections on the Sun. NASA's STEREO-B spacecraft, orbiting downstream from Venus, observed that the solar wind density at Earth's orbit dropped to the remarkably low figure of 0.1 particles per cubic cm and persisted at this value for an entire day."

"The observations show that the night side ionosphere moved outward to at least 15 000 km from Venus' centre over a period of only a few hours," said Markus Fraenz, also from the Max Planck Institute for Solar System Research, who was the team leader and a co-author of the paper. "It may possibly have extended for millions of kilometres, like an enormous tail."

"Although we cannot determine the full length of the night-side ionosphere, since the orbit of Venus Express provides limited coverage, our results suggest that Venus' ionosphere resembled the teardrop-shaped ionosphere found around comets, rather than the more symmetrical, spherical shape which usually exists."

"The side of Venus' ionosphere that faces away from the sun can billow outward like the tail of a comet, while the side facing the star remains tightly compacted, researchers said. ... "As this significantly reduced solar wind hit Venus, Venus Express saw the planet’s ionosphere balloon outwards on the planet’s ‘downwind’ nightside, much like the shape of the ion tail seen streaming from a comet under similar conditions," ESA officials said in a statement today (Jan. 29). It only takes 30 to 60 minutes for the planet's comet-like tail to form after the solar wind dies down. Researchers observed the ionosphere stretch to at least 7,521 miles (12,104 kilometers) from the planet, said Yong Wei, a scientist at the Max Planck Institute in Katlenburg, Germany who worked on this research."

"[B]ecause of the lack of a planetary magnetic field, the free hydrogen has been swept into interplanetary space by the solar wind. "

"The clouds of Venus are capable of producing lightning much like the clouds on Earth. The existence of lightning had been controversial since the first suspected bursts were detected by the Soviet Venera probes. In 2006–07 Venus Express clearly detected whistler mode waves, the signatures of lightning. Their intermittent appearance indicates a pattern associated with weather activity. The lightning rate is at least half of that on Earth. In 2007 the Venus Express probe discovered that a huge double atmospheric vortex exists at the south pole of the planet. "

"Another discovery made by the Venus Express probe in 2011 is that an ozone layer exists high in the atmosphere of Venus. "

"Venus has an extremely dense atmosphere, which consists mainly of carbon dioxide and a small amount of nitrogen. The atmospheric mass is 93 times that of Earth's atmosphere, while the pressure at the planet's surface is about 92 times that at Earth's surface—a pressure equivalent to that at a depth of nearly 1 kilometer under Earth's oceans. The density at the surface is 65 kg/m³ (6.5% that of water)."

Magnetohydrodynamics
"In 1967, Venera-4 found the Venusian magnetic field is much weaker than that of Earth. This magnetic field is induced by an interaction between the ionosphere and the solar wind, ... Venus's small induced magnetosphere provides negligible protection to the atmosphere against cosmic radiation. This radiation may result in cloud-to-cloud lightning discharges. "

"The weak magnetosphere around Venus means the solar wind is interacting directly with the outer atmosphere of the planet. Here, ions of hydrogen and oxygen are being created by the dissociation of neutral molecules from ultraviolet radiation. The solar wind then supplies energy that gives some of these ions sufficient velocity to escape the planet's gravity field. This erosion process results in a steady loss of low-mass hydrogen, helium, and oxygen ions, while higher-mass molecules, such as carbon dioxide, are more likely to be retained."

Meteorites
"While there is little or no water on Venus, there is a phenomenon which is quite similar to snow. The Magellan probe imaged a highly reflective substance at the tops of Venus's highest mountain peaks which bore a strong resemblance to terrestrial snow. This substance arguably formed from a similar process to snow, albeit at a far higher temperature. Too volatile to condense on the surface, it rose in gas form to cooler higher elevations, where it then fell as precipitation. The identity of this substance is not known with certainty, but speculation has ranged from elemental tellurium to lead sulfide (galena). "

Crater astronomy
"The absence of evidence of lava flow accompanying any of the visible caldera remains an enigma. The planet has few impact craters".

"After the Venera missions were completed, the prime meridian was redefined to pass through the central peak in the crater Ariadne. "

"Almost a thousand impact craters on Venus are evenly distributed across its surface. ... On Venus, about 85% of the craters are in pristine condition. ... Venusian craters range from 3 km to 280 km in diameter. No craters are smaller than 3 km, because of the effects of the dense atmosphere on incoming objects. Objects with less than a certain kinetic energy are slowed down so much by the atmosphere, they do not create an impact crater. Incoming projectiles less than 50 meters in diameter will fragment and burn up in the atmosphere before reaching the ground. "

Classical planets
The day of the week for Venus is Friday and its color is white.

5102 b2k, -3102 or 3102 BC, is the historical year assigned to a Hindu table of planets that does not include the classical planet Venus. "Vénus seule ne s'y trouvait pas." "Venus alone is not found there." "Babylonian astronomy, too, had a four-planet system. In ancient prayers the planets Saturn, Jupiter, Mars, and Mercury are invoked; the planet Venus is missing; and one speaks of "the four-planet system of the ancient astronomers of Babylonia." "

“That the planet Venus is missing will not startle anybody who knows the eminent importance of the four-planet system in the Babylonian astronomy” “Weidner supposes that Venus is missing in the list of planets because “she belongs to a triad with the moon and the sun.””

3581 b2k: "The Venus tablet of Ammisaduqa, dated 1581 BC, shows that the Babylonians understood that the two were a single object, referred to in the tablet as the "bright queen of the sky," and could support this view with detailed observations. "

"The Greeks thought of the two as separate stars, Phosphorus and Hesperus, until the time of Pythagoras in the sixth century BC. "

"The Romans designated the morning aspect of Venus as Lucifer, literally "Light-Bringer", and the evening aspect as Vesper".

title: Solar System, technical/Venus


text: File:FEcoll.jpg

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

title: File:FEcoll.jpg


text: File:CollSubExact.jpg

title: File:CollSubExact.jpg


text: File:AllWResExact.jpg

title: File:AllWResExact.jpg


text: File:PlotWRSols.png

title: File:PlotWRSols.png


text:
Nonlinear finite elements/Linear heat equation

Finite element methods are used to solve boundary value problems (BVP) or initial boundary value problems (IBVP) in engineering. BVPs are mathematical models of real-life situations. Such situations can be physical, biological, economic, and so on.

We will explore the problem of heat conduction and see how we build a finite element model and solve this problem. The first step will be to build a model. The model is the IBVP in the form of a partial differential equation or a variational problem. We will discuss whether the problem is well-posed. Then we will try to construct an approximate solution to the problem. Finally, we will discuss how good the solution is.

Construction of a Model
Figure 1 shows a region $$\Omega \in \mathbb{R}^{3}$$ through which heat is flowing. Points in the medium are represented by $$\mathbf{x}$$ with components $$(x_1, x_2, x_3)$$ with respect to the basis $$(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$$ and the origin. The heat flux into (or out of) the medium is $$\mathbf{q}(\mathbf{x},t)$$ where $$t$$ is the time. The heat flux takes place across part of the boundary of the medium ($$\Gamma_q$$). The temperature on the remainder of the boundary ($$\Gamma_T$$) has a known value $$\overline{T}$$. Heat sources inside the medium (for example, chemical reactions and plastic deformation) are given as a function $$s(\mathbf{x},t)$$.

The goal is to find the temperature distribution in the medium ($$T(\mathbf{x},t)$$) as time evolves.

A realistic model for this problems needs two components:

 * 1) A balance (or conservation) equation for energy.
 * 2) A constitutive equation for the medium.

Balance of energy
Let us assume that there are no sources of energy other than thermal energy. Then, the balance of energy states that:

 * Rate of change of thermal energy = heat generated by internal sources + heat flow into the body.

We have to convert this statement into mathematical form. To do this, consider an arbitrary part of the body ($$\Omega_1$$) with boundary $$\Gamma_1$$. Let $$\mathbf{n}$$ be the outward unit normal to the boundary.

The total thermal energy in $$\Omega_1$$ at a particular time is given by the heat capacity of the body. The heat capacity $$C_v$$ is the amount of energy needed to raise the temperature (of a unit mass of the body) by one unit.

Let the mass density be $$\rho(\mathbf{x})$$. Then, the total thermal energy in $$\Omega_1$$ at time $$t$$ is
 * $$\text{(2)} \qquad

\int_{\Omega_1} C_v(\mathbf{x})~\rho(\mathbf{x})~T(\mathbf{x},t)~dV~. $$

The total heat generated by internal sources is
 * $$\text{(3)} \qquad

\int_{\Omega_1} s(\mathbf{x},t)~dV~. $$

And the total heat flux into $$\Omega_1$$ is
 * $$\text{(4)} \qquad

- \int_{\Gamma_1} \mathbf{q}(\mathbf{x},t)\bullet\mathbf{n}~dA~. $$

Apply the divergence theorem to (4) to get
 * $$\text{(5)} \qquad

- \int_{\Gamma_1} \mathbf{q}(\mathbf{x},t)\bullet\mathbf{n}~dA = - \int_{\Omega_1} \boldsymbol{\nabla} \bullet \mathbf{q}(\mathbf{x},t)~dV~. $$

Put (2), (3), and (5) together, and apply the balance of energy to get
 * $$\text{(6)} \qquad

\cfrac{d}{dt}\int_{\Omega_1} C_v(\mathbf{x})~\rho(\mathbf{x})~T(\mathbf{x},t)~dV = \int_{\Omega_1} s(\mathbf{x},t)~dV - \int_{\Omega_1} \boldsymbol{\nabla} \bullet \mathbf{q}(\mathbf{x},t)~dV~. $$

The limits of integration are fixed. So we can write (6) as
 * $$\text{(7)} \qquad

\int_{\Omega_1} C_v(\mathbf{x})~\rho(\mathbf{x})~\frac{\partial T(\mathbf{x},t)}{\partial t}~dV = \int_{\Omega_1} s(\mathbf{x},t)~dV - \int_{\Omega_1} \boldsymbol{\nabla} \bullet \mathbf{q}(\mathbf{x},t)~dV~. $$

Since $$\Omega_1$$ is arbitrary, if the functions that appear in (7) are smooth enough, the equation is equivalent to
 * $$\text{(8)} \qquad

{ C_v(\mathbf{x})~\rho(\mathbf{x})~\frac{\partial T(\mathbf{x},t)}{\partial t} + \boldsymbol{\nabla} \bullet \mathbf{q}(\mathbf{x},t) = s(\mathbf{x},t)~. } $$ We can write equation (8) in index notation as

{ C_v~\rho~\dot{T} + q_{i,i} = s } $$ where we have assumed that the components are with respect to the Cartesian basis ($$\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3$$). In the following, whenever we use index notation, the components are assumed to be with respect to that Cartesian basis.

Constitutive equation
There are two unknowns in equation (8). These are $$T(\mathbf{x},t)$$ and $$\mathbf{q}(\mathbf{x},t)$$ and only one equation. Therefore, we need another equation that characterizes the material.

One possibility is the Fourier law which states that:

 *  the heat flux is linearly related to the temperature gradient.

In mathematical form,
 * $$\text{(9)} \qquad

{ \mathbf{q}(\mathbf{x},t) = -\boldsymbol{\kappa}(\mathbf{x})\bullet\boldsymbol{\nabla} T(\mathbf{x},t)~; \boldsymbol{\kappa} = \boldsymbol{\kappa}^T~. } $$ The quantity $$\boldsymbol{\kappa}$$ is a second-order tensor called the  thermal conductivity tensor. The minus sign shows that heat flows from hot to cold. Recall that a second-order tensor takes a vector to another vector (in this case a temperature gradient to a heat flux).

In index notation, we can write equation (9) as

{ q_i = -\kappa_{ij}~T_{,j} ~;\kappa_{ij} = \kappa_{ji} ~. } $$ If the region $$\Omega$$ is  homogeneous then $$\boldsymbol{\kappa}$$ is constant. If the region $$\Omega$$ is  isotropic, the thermal conductivity tensor takes the form

\kappa_{ij} (\mathbf{x}) = \kappa(\mathbf{x}) \delta_{ij} $$ where $$\kappa$$ is a scalar.

Heat equation
To get the heat equation, combine (9) with (8) to get
 * $$\text{(10)} \qquad

{ C_v(\mathbf{x})~\rho(\mathbf{x})~\frac{\partial T(\mathbf{x},t)}{\partial t} - \boldsymbol{\nabla} \bullet [\boldsymbol{\kappa}(\mathbf{x})\bullet\boldsymbol{\nabla} T(\mathbf{x},t)] = s(\mathbf{x},t)~. } $$ Rearrange to get
 * $$\text{(11)} \qquad

{ \frac{\partial T(\mathbf{x},t)}{\partial t} - \frac{1}{C_v(\mathbf{x})~\rho(\mathbf{x})} \boldsymbol{\nabla} \bullet [\boldsymbol{\kappa}(\mathbf{x})\bullet\boldsymbol{\nabla} T(\mathbf{x},t)] = Q(\mathbf{x},t) } $$ where $$Q(\mathbf{x},t) := s(\mathbf{x},t)/[C_v(\mathbf{x})~\rho(\mathbf{x})]$$.

In index notation, equation (11) is

{ \dot{T} - \cfrac{1}{C_v~\rho}~(\kappa_{ij}~T_{,j})_{,i} = Q ~. } $$ In expanded form, equation (11) reads
 * $$\text{(12)} \qquad

{ \begin{align} \frac{\partial T}{\partial t} - \frac{1}{C_v~\rho}& \left[\frac{\partial }{\partial x_1}\left( \kappa_{11} \frac{\partial T}{\partial x_1} + \kappa_{12} \frac{\partial T}{\partial x_2} +  \kappa_{13} \frac{\partial T}{\partial x_3} \right) + \right. \\ & ~\frac{\partial }{\partial x_2}\left( \kappa_{12} \frac{\partial T}{\partial x_1} + \kappa_{22} \frac{\partial T}{\partial x_2} +  \kappa_{23} \frac{\partial T}{\partial x_3} \right) + \\ & \left.~\frac{\partial }{\partial x_3}\left( \kappa_{13} \frac{\partial T}{\partial x_1} + \kappa_{23} \frac{\partial T}{\partial x_2} +  \kappa_{33} \frac{\partial T}{\partial x_3} \right)\right] = Q ~. \end{align} } $$ Equation (12) is the transient, inhomogeneous, heat equation.

Boundary conditions
Boundary conditions (BCs) are needed to make sure that we get a unique solution to equation (12).

The temperature is prescribed on $$\Gamma_T$$. Prescribed boundary conditions are also called  Dirichlet BCs or  essential BCs. In this case
 * $$\text{(13)} \qquad

{ T = \overline{T}(\mathbf{x}, t) ~\text{on}\Gamma_T~. } $$ The heat flux is given on the remainder of the boundary ($$\Gamma_q$$). Such flux boundary conditions are also known as  Neumann BCs or  natural BCs. The flux boundary condition is
 * $$\text{(14)} \qquad

{ \mathbf{q}(\mathbf{x},t)\bullet\mathbf{n}(\mathbf{x}) = \overline{q}(\mathbf{x},t) ~\text{on}\Gamma_q~. } $$ In index notation, the essential boundary condition is

{ q_i~n_i = \overline{q} \qquad\text{on}~\Gamma_q ~. } $$ Plug equation (9) into (14). We get
 * $$\text{(15)} \qquad

-(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\mathbf{n} = \overline{q}~. $$

If the region $$\Omega$$ is isotropic with thermal conductivity $$\kappa$$, we can define $$g := -\overline{q}/\kappa$$ and $$\partial T/\partial n := \boldsymbol{\nabla} T\bullet\mathbf{n}$$ (also called the normal derivative). Then the flux BC simplifies to

{ \frac{\partial T}{\partial n} = g ~\text{on} \Gamma_q~. } $$

If the boundary $$\Gamma_q$$ is insulated, then $$\overline{q}$$ = 0 = $$g$$.

Initial conditions
If the problem depends on time, we also need an initial condition for the temperature in the body,

{ T(\mathbf{x}, 0) = T_0(\mathbf{x})~. } $$

The complete model
The model initial boundary value problem (IBVP) for heat conduction is

$$\text{(16)} \qquad { \begin{align} & &\mathsf{The~initial~boundary~value~problem~for~heat~conduction}\\ & & \\ & \text{PDE:} & \frac{\partial T}{\partial t} - \frac{1}{C_v~\rho} \boldsymbol{\nabla} \bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T) = Q \text{in}\Omega, t > 0\quad\\ & \text{BCs:} & T = \overline{T}(\mathbf{x},t)\text{on}\Gamma_T \text{and} \mathbf{q}\bullet\mathbf{n} = \overline{q} \text{on}\Gamma_q\quad\\ & \text{IC:}~ & T(\mathbf{x},0) = T_0(\mathbf{x})\text{in}\Omega\quad \end{align} } $$


text:
Category:Heat transfer

This category contains all learning materials related to heat transfer, i.e., conduction, convection, radiation.

title: Category:Heat transfer


text: File:HeatCondProb.png

title: File:HeatCondProb.png


text:
Nonlinear finite elements/Steady state heat conduction

Steady state heat conduction
If the problem does not depend on time and the material is isotropic, we get the boundary value problem for steady state heat conduction. $$ { \begin{align} & &\mathsf{ The~ boundary~ value~ problem ~for~ steady~ heat ~conduction}\\ & & \\ & \text{PDE:} &- \frac{1}{C_v~\rho} \boldsymbol{\nabla} \bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla T)} = Q \text{in}\Omega\quad\\ & \text{BCs:} & T = \overline{T}(\mathbf{x})\text{on}\Gamma_T \text{and} \frac{\partial T}{\partial n} = g(\mathbf{x})\text{on}\Gamma_q\quad\\ \end{align} } $$

Poisson's equation
If the material is homogeneous the density, heat capacity, and the thermal conductivity are constant. Define the  thermal diffusivity as

k := \frac{\kappa}{C_v~\rho} $$ Then, the boundary value problem becomes $$ { \begin{align} & &\mathsf{ Poisson's~ equation}\\ & & \\ & \text{PDE:} &- k \nabla^2 T = Q \text{in}\Omega\quad\\ & \text{BCs:} & T = \overline{T}(\mathbf{x})\text{on}\Gamma_T \text{and} \frac{\partial T}{\partial n} = g(\mathbf{x})\text{on}\Gamma_q\quad\\ \end{align} } $$ where $$\nabla^2 T$$ is the  Laplacian

\nabla^2 T := \boldsymbol{\nabla} \bullet \boldsymbol{\nabla T} $$

Laplace's equation
Finally, if there is no internal source of heat, the value of $$Q$$ is zero, and we get Laplace's equation. $$ { \begin{align} & &\mathsf{ Laplace's~ equation}\\ & & \\ & \text{PDE:} &\nabla^2 T = 0 \text{in}\Omega\quad\\ & \text{BCs:} & T = \overline{T}(\mathbf{x})\text{on}\Gamma_T \text{and} \frac{\partial T}{\partial n} = g(\mathbf{x})\text{on}\Gamma_q\quad\\ \end{align} } $$

The Analogous Membrane Problem
The thin elastic membrane problem is another similar problem. See Figure 1 for the geometry of the membrane.

The membrane is thin and elastic. It is initially planar and occupies the 2D domain $$\Omega$$. It is fixed along part of its boundary $$\Gamma_u$$. A transverse force $$\mathbf{f}$$ per unit area is applied. The final shape at equilibrium is nonplanar. The final displacement of a point $$\mathbf{x}$$ on the membrane is $$\mathbf{u}(\mathbf{x})$$. There is no dependence on time.

The goal is to find the displacement $$\mathbf{u}(\mathbf{x})$$ at equilibrium.

It turns out that the equations for this problem are the same as those for the heat conduction problem - with the following changes:
 * The time derivatives vanish.
 * The balance of energy is replaced by the balance of forces.
 * The constitutive equation is replaced by a relation that states that the vertical force depends on the displacement gradient ($$\boldsymbol{\nabla} \mathbf{u}$$).

If the membrane if inhomogeneous, the boundary value problem is: $$ { \begin{align} & &\mathsf{ The~boundary~ value~ problem~ for~ membrane~ deformation}\\ & & \\ & \text{PDE:} &- \boldsymbol{\nabla} \bullet (E~\boldsymbol{\nabla \mathbf{u})} = Q \text{in}\Omega\quad\\ & \text{BCs:} & \mathbf{u} = \bar{\mathbf{u}}(\mathbf{x})\text{on}\Gamma_u \text{and} \frac{\partial u}{\partial n} = g(\mathbf{x})\text{on}\Gamma_t\quad\\ \end{align} } $$

For a homogeneous membrane, we get $$ { \begin{align} & &\mathsf{ Poisson's~ Equation}\\ & & \\ & \text{PDE:} &- E \nabla^2 \mathbf{u} = Q \text{in}\Omega\quad\\ & \text{BCs:} & \mathbf{u} = \bar{\mathbf{u}}(\mathbf{x})\text{on}\Gamma_u \text{and} \frac{\partial u}{\partial n} = g(\mathbf{x})\text{on}\Gamma_t\quad\\ \end{align} } $$

Note that the membrane problem can be formulated in terms of a problem of minimization of potential energy.


text: File:MembraneProb.png

title: File:MembraneProb.png


text:
Nonlinear finite elements/Weak form of Poisson equation

The Equivalent Variational Problem
So far, we have seen mathematical models in PDE form. That is not the only way of developing a model even though it is commonly used.

Many finite element models are based on an alternative form called the  variational problem. The problem is formulated as one where the goal is to minimize a functional.

The variational form and the partial differential equation form are equivalent for the same problem. Each can be derived from the other.

Let us look at a simplified problem - the Poisson equation. Simplify it further by assuming that the temperature is prescribed to be zero on the entire boundary.

Functional
The first stage is to find a functional that applies to the problem.This is the most difficult part of the process. One way of finding the right functional is to try out a number of different forms and see if they correspond to known differential equations. This process is relatively easy for linear problems but not for nonlinear problems.

Recall that a functional converts a function ($$v(\mathbf{x})$$) to a real number.

For the Poisson equation, we can find a functional by thinking in terms of the membrane problem. In that problem, the total potential energy is stored energy of the membrane minus the work done by the force. The functional can be written in the form
 * $$\text{(18)} \qquad

I(v) = \frac{1}{2}\int_{\Omega} |\boldsymbol{\nabla} v|^2~dV - \int_{\Omega} f~v~dV~. $$ where $$ v \equiv v(\mathbf{x})$$. The functional in equation (18) represents the total potential energy of the membrane. The first term of the functional represents the stored energy and the second term represents the work done by the applied forces. For the membrane the solution is the one that minimizes potential energy.

Since the same functional also works for the heat conduction problem (same PDE), we can minimize the functional to get at the temperature field that we seek.

Minimization problem
The Poisson equation for the heat conduction problem can then be written as

\text{Find a function}~T~\text{such that}~I(T) \le I(v)~\text{for all admissible}~v. $$

What are  admissible functions? In brief, these are functions that satisfy the following requirements (in this particular case)
 * They are continuously differentiable on $$\bar{\Omega}$$ where $$\bar{\Omega} := \Omega \cup \Gamma$$.
 * They satisfy the boundary condition $$v = 0$$ on $$\Gamma_T$$.

To summarize, the minimization problem is: $$ { \begin{align} &\mathsf{ Minimization~ Problem~ for~ Poisson's~ Equation}\\ & \\ &\text{Find a function}~T~\text{in}~\mathcal{X}~\text{that satisfies}\\ & \\ & I(T) \le I(v) ~\text{for all functions}~v~\text{in}~\mathcal{X} \\ & \\ &\text{where}~ I(v) = \frac{1}{2}\int_{\Omega} |\boldsymbol{\nabla} v|^2~dV - \int_{\Omega} f~v~dV\\ & \\ &\mathcal{X} = \{\text{Continuously differentiable functions on}~\bar{\Omega}~ \text{that vanish on}~\Gamma_T\} \end{align} } $$

Finding a minimum
Let us take a simple function of one variable, say $$h(x)$$. See Figure 1(a). If a minimum exists, we can find it by checking where the slope $$h^{'}(x) = 0$$.

A similar idea is used when we try to minimize a functional. See Figure 1(b).

For our problem, we have to find a minimizing function $$T(\mathbf{x})$$ from the space of all admissible functions $$v(\mathbf{x})$$.

The process of finding this minimizing function has the following steps:

 * Assume that a minimizing function exists.
 * Let this minimizing function be $$T(\mathbf{x})$$.
 * Choose an arbitrary function $$v(\mathbf{x}) \in \mathcal{X}$$.
 * Form a new function $$T(\mathbf{x}) + \epsilon v(\mathbf{x})$$.
 * If we consider only the dependence of $$\epsilon$$, we can treat the functional $$I(T + \epsilon v)$$ as a function only of $$\epsilon$$, that is, $$I(T + \epsilon v) \equiv F(\epsilon)$$.
 * Now, the minimum is achieved when $$\epsilon = 0$$. Hence,
 * $$\text{(19)} \qquad

\left. \cfrac{d}{d\epsilon} F(\epsilon) \right|_{\epsilon = 0} = 0 ~\Leftrightarrow~ \left. \cfrac{d}{d\epsilon} I(T + \epsilon v) \right|_{\epsilon = 0} = 0 $$
 * Recall that the functional to be minimized is
 * $$\text{(20)} \qquad

I(v) = \frac{1}{2}\int_{\Omega} |\boldsymbol{\nabla} v|^2~dV - \int_{\Omega} f~v~dV $$
 * Plug (20) into (19) and you will get
 * $$\text{(21)} \qquad

{ \int_{\Omega} \boldsymbol{\nabla} T\bullet\boldsymbol{\nabla} v~dV = \int_{\Omega} f~v~dV } $$
 * (Show that this is true.)

 * Equation (21) can be solved for $$T(\mathbf{x})$$ and is called the  variational boundary value problem.

The variational boundary value problem for the Poisson's equation form of the heat equation is $$ { \begin{align} &\mathsf{ Variational ~BVP ~for ~Poisson's ~Equation}\\ & \\ &\text{Find a function}~T~\text{in}~\mathcal{X}~\text{that satisfies}\\ & \\ & \int_{\Omega} \boldsymbol{\nabla} T\bullet\boldsymbol{\nabla} v~dV = \int_{\Omega} f~v~dV \quad \text{for all}~v \in \mathcal{X} \end{align} } $$

Do I NEED to know the functional?
If the PDE form of the BVP is known, then we can derive the variational form directly from the PDE!

This process has the following steps:

 * Start with the PDE BVP. Let's look at the Poisson equation for heat conduction as an example. The PDE is
 * $$\text{(23)} \qquad

- k \nabla^2 T = Q\text{or,} -\nabla^2 T = f $$
 * where $$f := Q/k$$.

 * Choose an arbitrary function $$v(\mathbf{x}) \in \mathcal{X}$$.
 * Multiply (23) by $$v$$ and integrate over $$\Omega$$. We get
 * $$\text{(24)} \qquad

- \int_{\Omega} (\nabla^2 T)~v~ dV = \int_{\Omega} f~v~ dV \text{or,} - \int_{\Omega} [\boldsymbol{\nabla} \bullet (\boldsymbol{\nabla T})]~v~ dV = \int_{\Omega} f~v~ dV~. $$
 * The left hand side of the equation has second derivatives. We want simplify things so that only first derivatives are involved. To do this, we use the identity

\boldsymbol{\nabla} \bullet (a\mathbf{b}) = (\boldsymbol{\nabla} a)\bullet\mathbf{b} + a (\boldsymbol{\nabla} \bullet b) $$
 * where $$a$$ is a scalar valued function ($$v$$ in our case), and $$\mathbf{b}$$ is a vector valued function (like $$\boldsymbol{\nabla T}$$). Apply this identity to the first term in equation (24). We get

- \int_{\Omega} [\boldsymbol{\nabla} \bullet (v\boldsymbol{\nabla T)} - (\boldsymbol{\nabla} v)\bullet(\boldsymbol{\nabla} T)] ~ dV = \int_{\Omega} f~v~ dV~. $$
 * Rearrange to get
 * $$\text{(25)} \qquad

- \int_{\Omega} \boldsymbol{\nabla} \bullet (v\boldsymbol{\nabla T)}~dV + \int_{\Omega} (\boldsymbol{\nabla} T)\bullet(\boldsymbol{\nabla} v)~dV = \int_{\Omega} f~v~ dV~. $$
 * The first term is (25) still has second derivatives. Apply the Gauss divergence theorem

\int_{\Omega} \boldsymbol{\nabla} \bullet \mathbf{b}~dV = \int_{\Gamma} \mathbf{n}\bullet\mathbf{b}~dA $$
 * where $$\mathbf{b}$$ is a vector valued function and $$\mathbf{n}$$ is the normal to the boundary $$\Gamma$$. The first term of (25) becomes

\int_{\Omega} \boldsymbol{\nabla} \bullet (v\boldsymbol{\nabla T)}~dV= \int_{\Gamma} \mathbf{n}\bullet(v\boldsymbol{\nabla} T)~dA = \int_{\Gamma} v(\boldsymbol{\nabla} T)\bullet\mathbf{n}~dA = \int_{\Gamma} v\frac{\partial T}{\partial n}~dA~. $$
 * The final form of (25) is called the  Weak Form of the PDE and is written as
 * $$\text{(26)} \qquad

{ - \int_{\Gamma} v\frac{\partial T}{\partial n}~dA + \int_{\Omega} (\boldsymbol{\nabla} T)\bullet(\boldsymbol{\nabla} v)~dV = \int_{\Omega} f~v~ dV~. } $$

At this point, we can break the integral over the boundary into two parts - an integral over $$\Gamma_T$$ and another over $$\Gamma_q$$, that is

- \int_{\Gamma} v\frac{\partial T}{\partial n}~dA = - \int_{\Gamma_T} v\frac{\partial T}{\partial n}~dA - \int_{\Gamma_q} v\frac{\partial T}{\partial n}~dA~. $$ Recall that the function $$v \in \mathcal{X}$$ has to meet the requirement that $$v = 0$$ on $$\Gamma_T$$. Therefore,

\int_{\Gamma_T} v\frac{\partial T}{\partial n}~dA= 0~. $$ We also have the boundary condition that the boundary is insulated ($$\partial T/\partial n = 0$$) on $$\Gamma_q$$. This implies that

\int_{\Gamma_q} v\frac{\partial T}{\partial n}~dA= 0~. $$ Then, the weak form becomes
 * $$\text{(27)} \qquad

{ \int_{\Omega} (\boldsymbol{\nabla} T)\bullet(\boldsymbol{\nabla} v)~dV = \int_{\Omega} f~v~ dV~. } $$ which is the same as the variational boundary value problem (21). We have arrived at the variational BVP without knowing the functional to be minimized!

Remarks

 * Note that the variational BVP implicitly contains the flux boundary condition and we do not have to apply it separately when solving the BVP. This type of BC is called a  natural boundary condition.
 * The prescribed boundary condition is not satisfied by the variational BVP and has to be applied explicitly when solving the BVP. This type of BC is called an  essential boundary condition.

You will often see the variational BVP (27) written in the form

{ a(T,v) = \langle f,v \rangle } $$ where
 * $$\begin{align}

a(T,v) & = \int_{\Omega} (\boldsymbol{\nabla} T)\bullet(\boldsymbol{\nabla} v)~dV \\ \langle f,v\rangle & = \int_{\Omega} f~v~ dV~. \end{align}$$ Here $$\langle f,v \rangle$$ is a linear functional $$f$$ acting on $$v$$ and $$a(.,.)$$ is an example of a  symmetric bilinear form. A bilinear form is an operator that maps a pair of elements to the real numbers, and which is linear in each of its slots. Let $$c_1$$ and $$c_2$$ be constants, and let $$u$$, $$v$$, and $$w$$ be functions. Then a bilinear form has the properties that
 * $$\begin{align}

a(c_1~u + c_2~v, w) & = c_1~a(u, w) + c_2~a(v,w) \\ (c_1~u + c_2~v, w) & = c_1~(u, w) + c_2~(v,w) \end{align}$$ A  symmetric bilinear form has the additional property that

a(u,v) = a(v,u) ~; ~ (u,v) = (v,u) $$

(Show that equation (27) is actually a symmetric bilinear form.)


text: File:MinFunction.png

title: File:MinFunction.png


text: File:MinFunctional.png

title: File:MinFunctional.png


text:
Nonlinear finite elements/Weak form of heat equation

Weak form of the heat equation
Let us now derive the weak form of the model of heat conduction in equations (16). It is more convenient to use the form of the governing equation given in equation (10). The equation is

\rho~C_v~\frac{\partial T}{\partial t} - \boldsymbol{\nabla} \bullet [\boldsymbol{\kappa}\boldsymbol{\nabla T]} = s~. $$

Let $$\mathcal{V}$$ be the space of weighting functions (or test functions). Then any $$w(\mathbf{x}) \in \mathcal{V}$$ ($$\mathbf{x} \in \bar{\Omega}$$) is continuously differentiable. The weighting functions also satisfy $$w(\mathbf{x}) = 0$$ on $$\Gamma_T$$. Recall that we called a similar set of functions $$\mathcal{X}$$ in our discussion of the Poisson problem. We can write

\mathcal{V} = \{ w(\mathbf{x}) ~|~ \mathbf{x} \in \bar{\Omega},~ w(\mathbf{x}) \in \mathcal{H}^1(\Omega), ~\text{and}~ w(\mathbf{x}) = 0~\text{on}~ \Gamma_T\} ~. $$

Let $$\mathcal{S}_t$$ be the set of trial solutions. Then any trial function $$T(\mathbf{x},t) \in \mathcal{S}_t$$ has to satisfy the essential boundary conditions on $$\Gamma_T$$. This is written as

\mathcal{S}_t = \{ T(\mathbf{x},t) ~|~ \mathbf{x} \in \bar{\Omega},~ t > 0, ~ T(\mathbf{x},t) \in \mathcal{H}^1(\Omega), ~\text{and}~ T(\mathbf{x},t) = \overline{T}(\mathbf{x},t) ~\text{on}~ \Gamma_T\} ~. $$

To get the weak form, we multiply the governing equation by the weighting function and integrate over the volume to get

 * $$\text{(28)} \qquad

\int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV - \int_{\Omega} \left[\boldsymbol{\nabla} \bullet (\boldsymbol{\kappa}\boldsymbol{\nabla T)}\right]~w~dV = \int_{\Omega} s~w~dV~. $$

The second term in the equation has second-order derivatives. We will convert these into first order derivatives using the divergence theorem and the identity

\boldsymbol{\nabla} \bullet (\varphi \mathbf{u}) = \varphi(\boldsymbol{\nabla} \bullet \mathbf{u}) + (\boldsymbol{\nabla} \varphi)\bullet\mathbf{u} ~. $$

From the identity we get

 * $$\text{(29)} \qquad

w\left[\boldsymbol{\nabla} \bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right] = \boldsymbol{\nabla} \bullet \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right] - (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T) ~. $$

Substitute (29) into the second term in (28) to get

 * $$\text{(30)} \qquad

\int_{\Omega} \left[\boldsymbol{\nabla} \bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]~w~dV = \int_{\Omega} \boldsymbol{\nabla} \bullet \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]~dV - \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV ~. $$

Apply the divergence theorem to the first term on the right hand side of (30). You will get

 * $$\text{(31)} \qquad

\int_{\Omega} \boldsymbol{\nabla} \bullet \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]~dV= \int_{\Gamma} \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]\bullet\mathbf{n}~dA $$

where $$\mathbf{n}$$ is the unit outward normal to the boundary $$\Gamma$$.

Since $$w(\mathbf{x}) = 0$$ on $$\Gamma_T$$, equation (31) becomes

 * $$\text{(32)} \qquad

\int_{\Omega} \boldsymbol{\nabla} \bullet \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]~dV= \int_{\Gamma_q} \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]\bullet\mathbf{n}~dA~. $$

Substitute (31) into (30) and (30) back into (28). You will get

 * $$\text{(33)} \qquad

\int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV - \int_{\Gamma_q} \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]\bullet\mathbf{n}~dA + \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV= \int_{\Omega} s~w~dV~. $$

After rearrangement, we get the  exact weak form of the heat equation
 * $$\text{(34)} \qquad

{ \int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV + \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV= \int_{\Omega} s~w~dV + \int_{\Gamma_q} \left[w~(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\right]\bullet\mathbf{n}~dA~. } $$

Recall that equation (15) gives us

-(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)\bullet\mathbf{n} = \overline{q} \text{on} \Gamma_q ~. $$

Therefore equation (34) can be written as

 * $$\text{(35)} \qquad

{ \int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV + \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV= \int_{\Omega} s~w~dV - \int_{\Gamma_q} w~\overline{q}~dA~. } $$

In more compact notation

 * $$\text{(36)} \qquad

{ \langle w,~\rho~C_v~\dot{T\rangle} + a(w,~T) = \langle w,~s\rangle - \langle w,~\overline{q}\rangle_{\Gamma} ~. } $$

Following the same process for the  initial condition, we get

 * $$\text{(37)} \qquad

{ \int_{\Omega} w~\left[\rho~C_v~T(0)\right]~dV = \int_{\Omega} w~\left[\rho~C_v~T_0\right]~dV ~. } $$ In compact notation,

{ \langle w,~\rho~C_v~T(0)\rangle = \langle w,~\rho~C_v~T_0\rangle ~. } $$

The variation initial boundary value problem for heat conduction can then be stated as follows.

$$ { \begin{align} & \mathsf{ Variational~ BVP ~for~ the ~Heat~ Equation}\\ & \\ \text{Find a function} & ~T(t) \in \mathcal{S}_t, t \in [0,\tau] ~\text{such that for all}~ w \in \mathcal{V}\\ & \\ & \int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV + \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV= \int_{\Omega} s~w~dV - \int_{\Gamma_q} w~\overline{q}~dA \\ & \int_{\Omega} w~\rho~C_v~T(0)~dV = \int_{\Omega} w~\rho~C_v~T_0~dV ~. \end{align} } $$

Well-posedness of the boundary value problem
Unless $$\Omega$$ has a simple geometry - square, spherical or cylindrical, it is very difficult to solve the BVP in closed form (using separation of variables, for instance). For some BVPs, it may not be possible  at all to get a closed form solution. In fact, it may not even be obvious that a solution exists or is unique or that the solution depends continuously on the data.

A  well-posed problem is one that satisfied the three conditions :
 * a solution exists.
 * the solution is unique.
 * the solution depends continuously on the data (that is, smallchanges in the data do not cause wild fluctuations in the solution).

BVPs can be used to get reliable results only when they are well-posed.

Let us look at an example. Recall the BVP for the Poisson equation.

$$ { \begin{align} & &\mathsf{ Poisson's~ equation}\\ & & \\ & \text{PDE:} &- \nabla^2 T = f \text{in}\Omega\quad\\ & \text{BCs:} & T = \overline{T}(\mathbf{x})\text{on}\Gamma_T \text{and} \frac{\partial T}{\partial n} = g(\mathbf{x})\text{on}\Gamma_q\quad\\ \end{align} } $$

Assume that $$\Gamma = \Gamma_q$$. That means that the heat flux is given on the entire boundary.

Suppose $$T$$ is a solution. Then $$T+c$$ is also a solution, where $$c$$ is a constant. To see why, compute the Laplacian of $$T+c$$.

\nabla^2 (T+c) = \boldsymbol{\nabla} \bullet [\boldsymbol{\nabla (T+c)]} = \boldsymbol{\nabla} \bullet (\boldsymbol{\nabla T)} = \nabla^2 T~. $$

So the solution $$T+c$$ satisfies the PDE.

How about the boundary conditions?Plug in $$T+c$$ into the boundary condition to get

\frac{\partial (T+c)}{\partial n} = [\boldsymbol{\nabla} (T+c)]\bullet\mathbf{n} = (\boldsymbol{\nabla} T)\bullet\mathbf{n} = \frac{\partial T}{\partial n} $$

So the boundary conditions are also satisfied. That means that if $$\Gamma_T = \varnothing$$, then the solution is  not unique. We can add any constant temperature to the body and the solution will be consistent with the governing equations.

We can also check the conditions under which a solution will  exist for this problem. Integrate the Poisson equation over $$\Omega$$ to get

\int_{\Omega} \nabla^2 T~dV = - \int_{\Omega} f~dV. $$

Apply the divergence theorem to the Laplacian. We get

\int_{\Omega} \nabla^2 T~dV = \int_{\Omega} \boldsymbol{\nabla} \bullet (\boldsymbol{\nabla T)}~dV = \int_{\Gamma} \mathbf{n}\bullet\boldsymbol{\nabla} T~dA = \int_{\Gamma} \frac{\partial T}{\partial n}~dA. $$

From the boundary condition,

\int_{\Gamma} \frac{\partial T}{\partial n}~dA = \int_{\Gamma} g~dA~. $$

Therefore, we have

\int_{\Gamma} g~dA = - \int_{\Omega} f~dV~\text{or,}~ { \int_{\Gamma} g~dA + \int_{\Omega} f~dV = 0~.} $$

The boxed equation is called a  compatibility condition and a solution does not exist unless this condition is satisfied.

In the context of the steady heat conduction problem, the compatibility condition says that the heat generated in the body must equal the heat flux. A similar (but more complicated) exercise can be used to show the existence and uniqueness of solutions for the full heat equation.


text:
Nonlinear finite elements/Solution of Poisson equation

Construction of Approximate Solutions
If we know that the problem is well-posed but does not have a closed form solution, we can go ahead and try to get an approximate solution. The finite element method is one way of getting at approximate solutions (among many other numerical methods).

The finite element method starts off with the variational form (or the weak form) of the BVP. The method is a special case of a class of methods called  Galerkin methods.

Finite element solution for the Poisson equation
Recall the variational boundary value problem for the Poisson equation: $$ { \begin{align} &\mathsf{ Variational~ BVP ~for ~Poisson's~ Equation}\\ & \\ &\text{Find a function}~T~\text{in}~\mathcal{X}~\text{that satisfies}\\ & \\ & \int_{\Omega} \boldsymbol{\nabla} T\bullet\boldsymbol{\nabla} v~dV = \int_{\Omega} f~v~dV \quad \text{for all}~v \in \mathcal{X} & \\ &\mathcal{X} = \{\text{Continuously differentiable functions on}~\bar{\Omega}~ \text{that vanish on}~\Gamma_T\} \end{align} } $$

The space $$\mathcal{X}$$ is continuous and an infinite number of functions could be chosen from this space of functions. In the finite element method, we choose a trial function from the space of approximate solutions $$\mathcal{X}_h$$ where $$\mathcal{X}_h \subset \mathcal{X}$$. A defining feature of these approximate trial solutions is that they are associated with a  mesh or  discretization of the domain $$\Omega$$. These functions also have the feature that they are finite dimensional with each dimension being associated with a  node on the mesh.

Assume that we are given $$\mathcal{X}_h$$. Let us choose a weighting function $$v_h \in \mathcal{X}_h$$ that satisfies $$v_h = 0$$ on $$\Gamma_T$$. We can choose another function $$T_h \in \mathcal{X}_h$$ as our trial solution. Since the boundary condition on $$\Gamma_T$$ is $$T = 0$$, both $$v_h$$ and $$T_h$$ can have the same form. In the next section, we will look at the general form of the heat equation where $$T \ne 0$$ on the boundary.

In finite element methods we choose trial solutions $$T_h$$ of the form

T_h(\mathbf{x}) = T_1 N_1(\mathbf{x}) + T_2 N_2(\mathbf{x}) + \dots + T_n N_n(\mathbf{x}) = \sum_{i=1}^n T_i N_i(\mathbf{x}) ~. $$ where $$T_1$$, $$T_2$$, $$\dots$$, $$T_n$$ are nodal temperatures which are constant on $$\bar{\Omega}$$. The functions $$N_1, N_2, \dots, N_n$$ form a basis that spans the subspace $$\mathcal{X}_h$$ and are known as  basis functions or  shape functions. Note that $$n$$ is the total number of nodes minus the number of nodes on $$\Gamma_T$$ where $$T$$ is specified.

Since the functions $$v_h$$ come from the same space of functions, we can represent them as

v_h(\mathbf{x}) = b_1 N_1(\mathbf{x}) + b_2 N_2(\mathbf{x}) + \dots + b_n N_n(\mathbf{x}) = \sum_{i=1}^n b_i N_i(\mathbf{x}) ~. $$ where $$b_1$$, $$b_2$$, $$\dots$$, $$b_n$$ are arbitrary constant on $$\bar{\Omega}$$ with the restriction that $$v_h = 0$$ on $$\Gamma_T$$.

If we plug in these finite dimensional forms of $$v$$ and $$T$$ into the variational BVP, we get an approximate form of the variational BVP which can be stated as: $$ { \begin{align} &\mathsf{Finite~ Element ~Variational~ BVP~ for~ Poisson's~ Equation}\\ & \\ \text{Find a function}& ~T_h~\text{in}~\mathcal{X}_h~\text{that satisfies}\\ & \\ & \int_{\Omega} \boldsymbol{\nabla} T_h\bullet\boldsymbol{\nabla} v_h~dV = \int_{\Omega} f~v_h~dV \quad \text{for all}~v_h \in \mathcal{X}_h ~. \end{align} } $$

After substituting the expressions for $$v_h$$ and $$T_h$$ in the variational BVP we get
 * $$\begin{align}

0 & = \int_{\Omega} \boldsymbol{\nabla} T_h\bullet\boldsymbol{\nabla} v_h~dV - \int_{\Omega} f~v_h~dV \\ & = \int_{\Omega} \boldsymbol{\nabla} (T_1 N_1+\dots+T_n N_n)\bullet \boldsymbol{\nabla} (b_1 N_1+\dots+b_n N_n)~dV - \int_{\Omega} f~(b_1 N_1+\dots+b_n N_n)~dV \\ & = \int_{\Omega} (T_1\boldsymbol{\nabla} N_1+\dots+T_n\boldsymbol{\nabla} N_n)\bullet (b_1\boldsymbol{\nabla} N_1+\dots+b_n\boldsymbol{\nabla} N_n)~dV - \int_{\Omega} f~(b_1 N_1+\dots+b_n N_n)~dV \\ & = \sum_{i,j=1}^n K_{ij} T_i b_j - \sum_{j=1}^n f_j b_j \\ & = \sum_{j=1}^n b_j \left[\sum_{i=1}^n K_{ij} T_i - f_j\right] \\ \end{align}$$ where,

{ K_{ij} = \int_{\Omega} \boldsymbol{\nabla} N_i\bullet\boldsymbol{\nabla} N_j~dV ~\text{and} f_{j} = \int_{\Omega} f N_j~dV~. } $$ In matrix form, we have
 * $$\text{(38)} \qquad

\mathbf{b}^{T} \left[\mathbf{K} \mathbf{T} - \mathbf{f}\right] = \mathbf{0} $$ where $$\mathbf{b}^T = [b_1, b_2, \dots, b_n]$$, $$\mathbf{K}$$ is a $$n \times n$$ symmetric matrix, $$\mathbf{T} = [T_1, T_2, \dots, T_n]$$ is a $$n \times 1$$ vector, and $$\mathbf{f}$$ is a $$n \times 1$$ vector.

Since $$\mathbf{b}$$ can be arbitrary, equation (38) can be further simplified to the form
 * $$\text{(39)} \qquad

{ \mathbf{K} \mathbf{T} = \mathbf{f} } $$ This system of equations has a solution since $$\mathbf{K}$$ is positive-definite and therefore has an inverse. Once the $$T_i$$s are known, the approximate solution can be found using

T_h(\mathbf{x}) = T_1 N_1(\mathbf{x}) + T_1 N_2(\mathbf{x}) + \dots + T_n N_n(\mathbf{x}) ~. $$ The functions $$N_1, \dots, N_n$$ have special forms in the finite element method that have the property that the quality of the approximation improves with an increase in the dimension $$n$$ of the basis.


text:
Nonlinear finite elements/Solution of heat equation

Finite element solution for the Heat equation
Let us now try to create a finite element approximation for the variational initial boundary value problem for the heat equation. This problem can be stated as $$ { \begin{align} & \mathsf{ Variational~ IBVP ~for~ the~ Heat~ Equation}\\ & \\ \text{Find a function} & ~T(t) \in \mathcal{S}_t, t \in [0,\tau] ~\text{such that for all}~ w \in \mathcal{V}\\ & \\ & \int_{\Omega} \left(\rho~C_v~\frac{\partial T}{\partial t}\right)~w~dV + \int_{\Omega} (\boldsymbol{\nabla} w)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T)~dV= \int_{\Omega} s~w~dV - \int_{\Gamma_q} w~\overline{q}~dA \\ & \int_{\Omega} w~\rho~C_v~T(0)~dV = \int_{\Omega} w~\rho~C_v~T_0~dV ~. \end{align} } $$

Let $$\mathcal{V}_h$$ be a finite dimensional approximation to $$\mathcal{V}$$ and let $$\mathcal{S}_h$$ be a finite dimensional approximation to $$\mathcal{S}$$. Let the weighting functions $$w_h \in \mathcal{V}_h$$ be such that $$w_h = 0$$ on $$\Gamma_T$$. Similarly, any other function $$v_h \in \mathcal{V}_h$$ also goes to zero on the boundary $$\Gamma_T$$.

Assume that the trial solutions $$T_h \in \mathcal{S}_h$$ can be represented as

T_h = v_h + \overline{T}_h \text{where}~ v_h \in \mathcal{V}_h ~. $$ The additional quantity $$\overline{T}_h$$ results in the satisfaction of the boundary condition $$T = \overline{T}$$ on $$\Gamma_T$$.

If we substitute these quantities into the variational initial boundary value problem, we get the Galerkin formulation. The steps in this process are shown below.

Approximate IBVP
The approximate form of the variational IBVP is

\int_{\Omega} \left(\rho~C_v~\frac{\partial T_h}{\partial t}\right)~w_h~dV + \int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} T_h)~dV= \int_{\Omega} s~w_h~dV - \int_{\Gamma_q} w_h~\overline{q}~dA ~. $$ Replacing $$T_h$$ with $$(v_h + \overline{T}_h)$$ we get

\int_{\Omega} \left(\rho~C_v~\frac{\partial (v_h + \overline{T}_h)}{\partial t}\right)~w_h~dV + \int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa} \bullet\boldsymbol{\nabla} (v_h + \overline{T}_h))~dV= \int_{\Omega} s~w_h~dV - \int_{\Gamma_q} w_h~\overline{q}~dA ~. $$ After expanding and rearranging terms, we get
 * $$\begin{align}

\int_{\Omega} w_h~\rho~C_v~\frac{\partial v_h}{\partial t}~dV + \int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} v_h)~dV = & \int_{\Omega} w_h~s~dV - \int_{\Gamma_q} w_h~\overline{q}~dA \\ & - \int_{\Omega} w_h~\rho~C_v~\frac{\partial \overline{T}_h}{\partial t}~dV - \int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} \overline{T}_h)~dV ~. \end{align}$$ In compact notation

\langle w_h,~\rho~C_v~\dot{v}\rangle_h + a(w_h,~v_h) = \langle w_h,~s\rangle - \langle w_h,~\overline{q}\rangle_{\Gamma} - \langle w_h,~\rho~C_v~\dot{\overline{T}}\rangle_h - a(w_h,~\overline{T}_h) ~. $$ Similarly, the initial condition can be written as

\int_{\Omega} w_h~\rho~C_v~T_h(0)~dV = \int_{\Omega} w_h~\rho~C_v~T_0~dV ~. $$ Substituting for $$T_h$$ and expanding out terms, we have

\int_{\Omega} w_h~\rho~C_v~v_h(0)~dV = \int_{\Omega} w_h~\rho~C_v~T_0~dV - \int_{\Omega} w_h~\rho~C_v~\overline{T}_h(0)~dV ~. $$ In compact form,

\langle w_h,~\rho~C_v~v_h(0)\rangle = \langle w_h,~\rho~C_v~T_0\rangle - \langle w_h,~\rho~C_v~\overline{T}_h(0)\rangle~. $$ Therefore the approximate form of the variational BVP can be written as $$\text{(40)} \qquad { \begin{align} & \mathsf{ Approximate~ Variational~ IBVP~ for~ the~ Heat~ Equation}\\ & \\ \text{Find a function} & ~T_h(t) = v_h + \overline{T}_h \in \mathcal{S}_h, t \in [0,\tau] ~\text{such that for all}~ w_h \in \mathcal{V}_h\\ & \\ \int_{\Omega} w_h~\rho~C_v~\frac{\partial v_h}{\partial t}~dV & + \int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} v_h)~dV = \int_{\Omega} w_h~s~dV - \int_{\Gamma_q} w_h~\overline{q}~dA \\ & - \int_{\Omega} w_h~\rho~C_v~\frac{\partial \overline{T}_h}{\partial t}~dV -\int_{\Omega} (\boldsymbol{\nabla} w_h)\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} \overline{T}_h)~dV ~. \\ \int_{\Omega} w_h~\rho~C_v~v_h(0)~dV & = \int_{\Omega} w_h~\rho~C_v~T_0~dV- \int_{\Omega} w_h~\rho~C_v~\overline{T}_h(0)~dV ~. \end{align} } $$ Note that the derivatives with respect to time are still continuous in this approximate variational BVP.

Finite element approximation
We now discretize our domain into elements $$\Omega_e$$ where $$1 < e < E$$ and $$E$$ is the total number of elements. Nodes may exist anywhere within the element domains. But the most common locations are at the element vertices and along the interelement boundaries.

Let the set of global node numbers be

\eta = \{1, 2, \dots, N\} ~. $$

Let the set nodes at which a temperature (essential BC) is prescribed be

\eta^T \subset \eta ~. $$

Then the set of nodes at which a temperature is to be determined is the complement

\eta - \eta^T ~. $$ Let $$n$$ be the number of nodes in $$\eta - \eta^T$$. Then the number of equations to be solved is also $$n$$. See Figure 1 to get an idea about what these sets of nodes refer to.

As we have seen before, a typical weighting function $$w_h \in \mathcal{V}_h$$ is assumed to have the form
 * $$\text{(42)} \qquad

{ w_h(\mathbf{x}) = \sum_{i=1}^n b_i N_i(\mathbf{x}) \equiv \sum_{i\in\eta-\eta^T} b_i N_i(\mathbf{x})~. } $$ Note that $$w_h = 0$$ only if $$b_i = 0$$ for every $$i \in \eta - \eta^T$$. On the other hand $$w_h = 0$$ for every node in $$\eta^T$$.

Similarly, a typical trial solution $$v_h \in \mathcal{V}_h$$ can be written as
 * $$\text{(43)} \qquad

{ v_h(\mathbf{x},t) = \sum_{i\in\eta-\eta^T} T_i(t) N_i(\mathbf{x}) } $$ where $$T_i(t)$$ is the unknown temperature at node $$i$$ at time $$t$$. We also often represent the approximate form of the essential BCs in a similar way, i.e.,
 * $$\text{(44)} \qquad

{ \overline{T}_h(\mathbf{x},t) = \sum_{i\in \eta^T} \overline{T}_i(t) N_i(\mathbf{x}), \qquad \overline{T}_i(t) = \overline{T}(\mathbf{x}_i,t) ~. } $$

Substitute equations (42), (43), and (44) into the first of equations (40). You will get
 * $$\begin{align}

\int_{\Omega} \left(\sum_j b_j N_j\right)~ & \rho~C_v~\left(\sum_i \frac{\partial T_i}{\partial t} N_i\right)~dV + \int_{\Omega} \left(\sum_j b_j \boldsymbol{\nabla} N_j\right)\bullet \left[\sum_i T_i \left(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_i\right)\right]~dV =\\ & \int_{\Omega} \left(\sum_j b_j N_j\right)~s~dV - \int_{\Gamma_q} \left(\sum_j b_j N_j\right)~\overline{q}~dA\\ & - \int_{\Omega} \left(\sum_j b_j N_j\right)~ \rho~C_v~\left(\sum_k \frac{\partial \overline{T}_k}{\partial t} N_k\right)~dV \\ & - \int_{\Omega} \left(\sum_j b_j \boldsymbol{\nabla} N_j\right)\bullet \left[\sum_k \overline{T}_k \left(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_k\right)\right]~dV , \qquad i,j \in \eta - \eta^T, \quad k \in \eta^T~. \end{align}$$ Assuming that $$\Omega$$ does not change with time, we can take the constants and time-dependent parameters outside the integrals to get
 * $$\begin{align}

\sum_j b_j & \left[ \sum_i \frac{\partial T_i}{\partial t}\left(\int_{\Omega}\rho~C_v~N_i~N_j~dV\right) + \sum_i T_i\left(\int_{\Omega}\boldsymbol{\nabla} N_j \bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_i)~dV \right)\right] =\\ & \sum_j b_j \left[ \int_{\Omega} N_j~s~dV - \int_{\Gamma_q} N_j~\overline{q}~dA - \sum_k \frac{\partial \overline{T}_k}{\partial t} \left(\int_{\Omega} \rho~C_v~N_j~N_k~dV\right)\right. \\ & - \left.\sum_k \overline{T}_k \left(\int_{\Omega} \boldsymbol{\nabla} N_j\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_k)~dV \right)\right] , \qquad i,j \in \eta - \eta^T, \quad k \in \eta^T~. \end{align}$$ Using the usual argument about the arbitrary nature of $$b_j$$, we get

{ \begin{align} \sum_i & \frac{\partial T_i}{\partial t}\left(\int_{\Omega}\rho~C_v~N_i~N_j~dV\right) + \sum_i T_i\left(\int_{\Omega}\boldsymbol{\nabla} N_j\bullet (\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_i)~dV \right) = \\ & \int_{\Omega} N_j~s~dV - \int_{\Gamma_q} N_j~\overline{q}~dA - \sum_k \frac{\partial \overline{T}_k}{\partial t} \left(\int_{\Omega} \rho~C_v~N_j~N_k~dV\right) \\ & - \sum_k \overline{T}_k \left(\int_{\Omega} \boldsymbol{\nabla} N_j\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_k)~dV \right) , \qquad i,j \in \eta - \eta^T, \quad k \in \eta^T ~. \end{align} } $$ In compact notation,

{ \sum_i \langle N_j,~\rho C_v N_i\rangle~\dot{T_i} + \sum_i a(N_j, N_i)~T_i = \langle N_j,~s\rangle - \langle N_j,~\overline{q}\rangle_{\Gamma} - \sum_k \langle N_j,~\rho C_v N_k\rangle~\dot{\overline{T}}_k - \sum_k a(N_j, N_k)~\overline{T}_k } $$ where $$i, j \in \eta - \eta^T$$, and $$k \in \eta^T$$.

Let us rename parts of the above equation as follows:
 * $$\begin{align}

M_{ij} & := \langle N_i,~\rho C_v N_j\rangle \\ K_{ij} & := a(N_i, N_j) \\ f_{j} & := \langle N_j,~s\rangle - \langle N_j,~\overline{q}\rangle_{\Gamma} - \sum_k \langle N_j,~\rho C_v N_k\rangle~\dot{\overline{T}}_k - \sum_k a(N_j, N_k)~\overline{T}_k ~. \end{align}$$ Then we get

{ \sum_i M_{ji}~\dot{T_i} + \sum_i K_{ji}~T_i = f_j ~. } $$ In matrix form,
 * $$\text{(45)} \qquad

{ \mathbf{M}~\dot{\mathbf{T}} + \mathbf{K}~\mathbf{T} = \mathbf{f} ~. } $$ As before, we can break up the global matrices into sums over elements. Thus,
 * $$\begin{align}

\mathbf{M} & = \sum_{e=1}^E \mathbf{M}^e \\ \mathbf{K} & = \sum_{e=1}^E \mathbf{K}^e \\ \mathbf{f} & = \sum_{e=1}^E \mathbf{f}^e \end{align}$$ where
 * $$\begin{align}

M^e_{ij} & = \int_{\Omega^e} \rho~C_v~N_i~N_j~dV~,\\ K^e_{ij} & = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_j)~dV \equiv \int_{\Omega^e} \mathbf{B}^T_i~\mathbf{D}~\mathbf{B}_j~dV ~, \\ f^e_j & = \int_{\Omega^e} N_j~s~dV - \int_{\Gamma^e_q} N_j~\overline{q}~dA - \sum_{k=1}^{n^e} \left[ M^e_{jk} \dot{\overline{T}}_k - K^e_{jk} \overline{T}_k \right] \end{align}$$ and $$n^e$$ is the number of nodes in an element.

We can do the same for the initial conditions. However, a more common approach is to specify the values of $$\mathbf{T}_0$$ directly at the nodes instead of going through an assembly process.

Computing M, K, f
The finite element system of equations at time $$t$$:
 * $$\begin{align}

M^e_{ij} & = \int_{\Omega^e} \rho~C_v~N_i~N_j~dV~,\\ K^e_{ij} & = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet({\boldsymbol{\kappa}\bullet{\boldsymbol{\nabla} N_j})}~dV \equiv \int_{\Omega^e} \mathbf{B}^T_i~\mathbf{D}~\mathbf{B}_j~dV ~, \\ f^e_j & = \int_{\Omega^e} N_j~s~dV - \int_{\Gamma^e_q} N_j~\overline{q}~dA - \sum_{k=1}^{n^e} \left[ M^e_{jk} \dot{\overline{T}}_k - K^e_{jk} \overline{T}_k \right] \end{align}$$

Coordinate Transformation

 * $$\begin{align}

x & = x_1^e~N_1^e(\xi,\eta) + x_2^e~N_2^e(\xi,\eta) + x_3^e~N_3^e(\xi,\eta) + x_4^e~N_4^e(\xi,\eta) \\ & \\ y & = y_1^e~N_1^e(\xi,\eta) + y_2^e~N_2^e(\xi,\eta) + y_3^e~N_3^e(\xi,\eta) + y_4^e~N_4^e(\xi,\eta) \end{align}$$

Integrating Stiffness Matrix
Stiffness matrix terms:

K^e_{ij} = \int_{\Omega^e}\boldsymbol{\nabla} N_i\bullet(\boldsymbol{\kappa}\bullet\boldsymbol{\nabla} N_j)~dV $$

Index notation:

K^e_{ij} = \int_{\Omega^e} N_{i,m}~\kappa_{mn}~N_{j,n}~dV $$

Expanded out (in 2D) (assume $$\kappa_{12} = 0$$):

K^e_{ij} = \int_{\Omega^e} \left( \frac{\partial N_i}{\partial x}\kappa_{11}\frac{\partial N_j}{\partial y} + \frac{\partial N_i}{\partial x}\kappa_{22}\frac{\partial N_j}{\partial y}\right)~dx~dy $$

Transformation
Chain Rule:
 * $$\begin{align}

\frac{\partial N_i}{\partial \xi} & = \frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \xi} + \frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial N_i}{\partial \eta} & = \frac{\partial N_i}{\partial x}\frac{\partial x}{\partial \eta} + \frac{\partial N_i}{\partial y}\frac{\partial y}{\partial \eta} \\ \end{align}$$

Matrix notation:

\begin{bmatrix} \frac{\partial N_i}{\partial \xi} \\ \\ \frac{\partial N_i}{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & & \frac{\partial y}{\partial \xi} \\ \\ \frac{\partial x}{\partial \eta} & & \frac{\partial y}{\partial \eta} \end{bmatrix} \begin{bmatrix} \frac{\partial N_i}{\partial x}\\ \\ \frac{\partial N_i}{\partial y} \end{bmatrix} = \mathbf{J}^e \begin{bmatrix} \frac{\partial N_i}{\partial x}\\ \\ \frac{\partial N_i}{\partial y} \end{bmatrix} $$

Inverse Transformation

\begin{bmatrix} \frac{\partial N_i}{\partial x} \\ \\ \frac{\partial N_i}{\partial y} \end{bmatrix} = (\mathbf{J}^e)^{-1} \begin{bmatrix} \frac{\partial N_i}{\partial \xi}\\ \\ \frac{\partial N_i}{\partial \eta} \end{bmatrix} ~; \mathbf{J}^{*} = (\mathbf{J}^e)^{-1}~\text{exists if}~\det(\mathbf{J}^e) \ne 0 ~. $$

with
 * $$\begin{align}

\frac{\partial x}{\partial \xi} & = x_1^e~\frac{\partial N_1^e}{\partial \xi} + x_2^e~\frac{\partial N_2^e}{\partial \xi} + x_3^e~\frac{\partial N_3^e}{\partial \xi} + x_4^e~\frac{\partial N_4^e}{\partial \xi} \\ \frac{\partial y}{\partial \xi} & = y_1^e~\frac{\partial N_1^e}{\partial \xi} + y_2^e~\frac{\partial N_2^e}{\partial \xi} + y_3^e~\frac{\partial N_3^e}{\partial \xi} + y_4^e~\frac{\partial N_4^e}{\partial \xi} \\ \frac{\partial x}{\partial \eta} & = x_1^e~\frac{\partial N_1^e}{\partial \eta} + x_2^e~\frac{\partial N_2^e}{\partial \eta} + x_3^e~\frac{\partial N_3^e}{\partial \eta} + x_4^e~\frac{\partial N_4^e}{\partial \eta} \\ \frac{\partial y}{\partial \eta} & = y_1^e~\frac{\partial N_1^e}{\partial \eta} + y_2^e~\frac{\partial N_2^e}{\partial \eta} + y_3^e~\frac{\partial N_3^e}{\partial \eta} + y_4^e~\frac{\partial N_4^e}{\partial \eta} \end{align}$$

Returning to the integration of the stiffness matrix terms:

K^e_{ij} = \int_{\Omega^e} \left( \frac{\partial N_i}{\partial x}\kappa_{11}\frac{\partial N_j}{\partial y} + \frac{\partial N_i}{\partial x}\kappa_{22}\frac{\partial N_j}{\partial y}\right)~dx~dy $$

In parent element coordinates:
 * $$\begin{align}

K^e_{ij} & = \int_{-1}^1 \left( \kappa_{11} \left( J^{*}_{11}\frac{\partial N_i}{\partial \xi} + J^{*}_{12}\frac{\partial N_i}{\partial \eta}\right) + \left( J^{*}_{11}\frac{\partial N_j}{\partial \xi} + J^{*}_{12}\frac{\partial N_j}{\partial \eta}\right) + \right. \\ & \left. \kappa_{22} \left( J^{*}_{21}\frac{\partial N_i}{\partial \xi} + J^{*}_{22}\frac{\partial N_i}{\partial \eta}\right) + \left( J^{*}_{21}\frac{\partial N_j}{\partial \xi} + J^{*}_{22}\frac{\partial N_j}{\partial \eta}\right) \right) ~d\xi~d\eta \end{align}$$

Gaussian Quadrature

\int_{-1}^1 F_{ij}(\xi,\eta)~d\xi~d\eta = \sum_{I=1}^M\sum_{J=1}^NF_{IJ}(\xi_I,\eta_J)~W_I~W_J $$

Gaussian Quadrature Example
Find the constants Co, C1, and X so that the quadrature formula

\int_{0}^1 f (x)\,dx = Co f(0)+ C1 f(x_1). $$

This has  the  highest  possible  degree  of  precision.

Solution.

Since there  are  three  unknowns,  Co,  C1  and  X1,  we  will  expect  the  formula  to  be  exact  for

f (x) = 1,  x,  and \,     \ x^2

$$ Thus

f (x)= 1,\ \int_{0}^1 f (x)\,dx= 1 = C_0 + C_1\  \qquad                  \ Equation 1,   \qquad

f (x)= x,\ \int_{0}^1 f (x)\,dx= \frac{1}{2} = C_1x_1\  \qquad             \ Equation 2. \qquad

f (x)= x^2,\ \int_{0}^1 f (x)\,dx= \frac{1}{3} = C_1x_1^2. $$

Equation 2 and 3 will yield.

\frac{c_1x_1}{c_1x_1^2} = x_1 = \frac{2}{3}. \qquad  c_1=\frac{3}{4}. \qquad c_0= \frac{1}{4}. $$ Hence

\int_{0}^1 f (x)\,dx = \frac{1}{4}f(0)+ \frac{3}{4}f(\frac{2}{3}) $$ Now,

f(x)=x^3.\qquad \int_{0}^1 x^3\,dx = \frac{1}{4} $$ And

\frac{1}{4}(0)^3 + \frac{3}{4}(\frac{2}{3})^3=\frac{2}{9}. $$ Thus the degree of the precision is 2

Robin boundary conditions
A similar approach can be used when Robin boundary conditions are needed to be applied. Let us assume that the Robin boundary conditions take the form

\alpha~T + \beta~(\boldsymbol{\nabla} T \cdot \mathbf{n}) = h \qquad \qquad \text{on}~\Gamma ~. $$

Recall from the previous section that the boundary term in the weak form of the heat equation is

\int_{\Gamma} w~\overline{q}~dA $$ where

\overline{q} = (\boldsymbol{\kappa}\cdot\boldsymbol{\nabla}T)\cdot\mathbf{n} ~. $$. Let us assume, for simplicity, that the material is isotropic. In that case $$\boldsymbol{\kappa}$$ becomes a scalar quantity $$\kappa\,$$ and we can write

\overline{q} = \kappa~(\boldsymbol{\nabla}T\cdot\mathbf{n}) ~. $$

Now, we can write the Robin boundary condition as

\boldsymbol{\nabla} T \cdot \mathbf{n} = \cfrac{h}{\beta} - \cfrac{\alpha}{\beta}~T \equiv \hat{h} - \hat{\alpha}~T \qquad \qquad \text{on}~\Gamma ~. $$ Using this expression we can get of the flux term in $$\overline{q}$$ to get

\overline{q} = \kappa~(\hat{h} - \hat{\alpha}~T) ~. $$

Then the boundary term takes the form

\int_{\Gamma} w~\overline{q}~dA = \int_{\Gamma} w~\kappa~(\hat{h} - \hat{\alpha}~T)~dA ~. $$ If we ignore the contributions to the trial function from the Dirichlet boundary conditions we can write

w_h(\mathbf{x}) = \sum_{j=1}^n b_j N_j(\mathbf{x}) ~; T_h = v_h(\mathbf{x},t) = \sum_{i=1}^n T_i(t) N_i(\mathbf{x}) $$ Plugging these expressions into the boundary term leads to

\begin{align} \sum_j b_j \int_{\Gamma} N_j~\overline{q}~dA & = \sum_j b_j \int_{\Gamma} N_j~\kappa~(\hat{h} - \hat{\alpha}~\sum_i T_i N_i)~dA \\ & = \sum_j b_j \left(\int_{\Gamma} \kappa~\hat{h}~N_j~dA -                           \sum_i T_i~\int_{\Gamma} \kappa~\hat{\alpha}~N_i~N_j~dA \right) \end{align} $$ Now recall that spatially discretized equation for transient heat conduction can be expressed as (in simplified form; after getting rid of contributions due to Dirichlet boundary conditions and with a scalar $$\kappa\,$$)
 * $$\begin{align}

\sum_j b_j & \left[ \sum_i \frac{\partial T_i}{\partial t}\left(\int_{\Omega}\rho~C_v~N_i~N_j~dV\right) + \sum_i \kappa~T_i\left(\int_{\Omega}\boldsymbol{\nabla} N_j \bullet \boldsymbol{\nabla} N_i~dV \right)\right] =\\ & \sum_j b_j \left[ \int_{\Omega} N_j~s~dV - \int_{\Gamma} N_j~\overline{q}~dA\right] \end{align}$$ Substituting the last term on the right hand side with the new expression for the boundary term, we have
 * $$\begin{align}

\sum_j b_j & \left[ \sum_i \frac{\partial T_i}{\partial t}\left(\int_{\Omega}\rho~C_v~N_i~N_j~dV\right) + \sum_i \kappa~T_i\left(\int_{\Omega}\boldsymbol{\nabla} N_j \bullet \boldsymbol{\nabla} N_i~dV \right)\right] =\\ & \sum_j b_j \left[ \int_{\Omega} N_j~s~dV - \kappa~\int_{\Gamma} \hat{h}~N_j~dA + \sum_i \kappa~T_i~\int_{\Gamma} \hat{\alpha}~N_i~N_j~dA \right] \end{align}$$ Invoking the arbitrariness of $$b_j$$ and rearranging, we get
 * $$\begin{align}

\sum_i \frac{\partial T_i}{\partial t}\left(\int_{\Omega}\rho~C_v~N_i~N_j~dV\right) & + \sum_i \kappa~T_i\left(\int_{\Omega}\boldsymbol{\nabla} N_j \bullet \boldsymbol{\nabla} N_i~dV - \int_{\Gamma}\hat{\alpha}~N_i~N_j~dA \right) = \\ & \int_{\Omega} s~N_j~dV - \int_{\Gamma} \kappa~\hat{h}~N_j~dA ~. \end{align}$$ Define

\begin{align} M_{ji} & := \int_{\Omega}\rho~C_v~N_j~N_i~dV \\ K_{ji} & := \int_{\Omega}\boldsymbol{\nabla} N_j \bullet \boldsymbol{\nabla} N_i~dV - \int_{\Gamma}\hat{\alpha}~N_j~N_i~dA \\ f_j & := \int_{\Omega} s~N_j~dV - \int_{\Gamma} \kappa~\hat{h}~N_j~dA ~. \end{align} $$ Then the finite element system of equations is

\sum_i M_{ji}~\dot{T}_i + \sum_i K_{ji}~T_i = f_j ~. $$ or, in matrix form

\mathbf{M}~\dot{\mathbf{T}} + \mathbf{K}~\mathbf{T} = \mathbf{f} ~. $$ Note that the $$\mathbf{K}$$ matrix and the $$\mathbf{f}$$ vector are different from the case where there are no Robin boundary conditions. The boundary terms in the discretized system of equations can be computed in the usual manner and will add terms to the diagonal of the $$\mathbf{K}$$ matrix.


text: File:HeatCondProbFEM.png

title: File:HeatCondProbFEM.png


text:
Nonlinear finite elements/Heat equation time integration

Solving the Matrix Problem
The matrix equations for the Poisson problem do not involve time and can be solved directly using either direct or iterative methods for solving systems of equations. You have done that in your introductory course on finite elements.

For the time-dependent heat equation, a few extra steps are needed. This is because the equations we have developed so far still have continuous time derivatives which need to be approximated.

Recall equations (45)
 * $$\text{(45)} \qquad

\mathbf{M}~\dot{\mathbf{T}} + \mathbf{K}~\mathbf{T} = \mathbf{f} ~. $$

These equations are a coupled system of first-order ordinary differential equations rather than a system of algebraic equations.

One way of solving the system of differential equations (45) is to use the  generalized midpoint rule.

Generalized midpoint rule
Consider the following initial value problem

\begin{align} \dot{T}(t) & = g(T(t)) \qquad \text{for all} \qquad t \in [0, \tau] \\ T(0) & = T_n ~, \end{align} $$ where $$g : \mathbb{R}^{} \rightarrow \mathbb{R}^{}$$ is a smooth function.

The generalized midpoint rule can be used to solve this differential equation in an approximate manner.

Let us discretize the domain $$[0, \tau]$$ into increments of $$\Delta t$$. Let $$T_{n+1} := T(t_{n+1})$$ be the approximation to the exact value of $$T(t)$$ at time $$t_{n+1} = t_n + \Delta t$$.

The integration rule for the generalized midpoint rule is

{ \begin{align} T_{n+1}& = T_n + \Delta t~g(T_{n+\theta}) \\ g(T_{n+\theta}) & := (1- \theta)~g(T_n) + \theta~g(T_{n+1}) \qquad \theta \in [0,1] ~. \end{align} } $$

When we choose $$\theta = 0$$, we get

{ T_{n+1} = T_n + \Delta t~g(T_n);\qquad \theta = 0 \implies \text{Forward Euler - fully explicit.} } $$

When we choose $$\theta = 1$$, we get

{ T_{n+1} = T_n + \Delta t~g(T_{n+1});\qquad \theta = 1 \implies \text{Backward Euler - fully implicit.} } $$

When we choose $$\theta = 1/2$$, we get

{ T_{n+1} = T_n + \Delta t~\left(\cfrac{g(T_n) + g(T_{n+1})}{2}\right); \qquad \theta = \frac{1}{2} \implies \text{Midpoint rule, Crank-Nicolson.} } $$

Generalized midpoint rule for Heat equation
If we apply the generalized midpoint rule to the system of differential equations

\mathbf{M}~\dot{\mathbf{T}} + \mathbf{K}~\mathbf{T} = \mathbf{f} $$ we get
 * $$\text{(46)} \qquad

{ \begin{align} \mathbf{M}~\dot{\mathbf{T}}_{n+1} + \mathbf{K}~\mathbf{T}_{n+1} &= \mathbf{f}_{n+1} \\ \mathbf{T}_{n+1} & = \mathbf{T}_n + \Delta t ~\dot{\mathbf{T}}_{n + \theta} \\ \dot{\mathbf{T}}_{n+\theta} & = (1 - \theta)~\dot{\mathbf{T}}_n + \theta~\dot{\mathbf{T}}_{n+1} \end{align} } $$ The computational problem is to find $$\mathbf{T}_{n+1}$$ and $$\dot{\mathbf{T}}_{n+1}$$ given $$\mathbf{T}_n$$ and $$\dot{\mathbf{T}}_n$$.

We start at $$t = 0$$ at which time the initial condition $$\mathbf{T}_0$$ is given. In that case, the value of $$\dot{\mathbf{T}}_0$$ can be calculated using

\dot{\mathbf{T}}_0 = \mathbf{M}^{-1}(\mathbf{f}_0 - \mathbf{K}~\mathbf{T}_0) ~. $$ Implementation of the algorithm for subsequent times may take various forms.

Let us combine the second and third equations in (46) to get
 * $$\text{(47)} \qquad

\mathbf{T}_{n+1} = \mathbf{T}_n + (1 - \theta)~\Delta t~\dot{\mathbf{T}}_n + \theta~\Delta t~\dot{\mathbf{T}}_{n+1} ~. $$ If we collect the known quantities at time $$t_n$$, we can write equation (47) as
 * $$\text{(48)} \qquad

{ \mathbf{T}_{n+1} = \tilde{\mathbf{T}}_{n+1} + \theta~\Delta t~\dot{\mathbf{T}}_{n+1} } $$ where

\tilde{\mathbf{T}}_{n+1} := \mathbf{T}_n + (1 - \theta)~\Delta t~\dot{\mathbf{T}}_n ~. $$ The quantity $$\tilde{\mathbf{T}}_{n+1}$$ is often called the  predictor value of $$\mathbf{T}_{n+1}$$.

At this stage we can proceed in two ways (or more).

The v-form of the Generalized Midpoint Rule.
Substitute equation (48) into the first equation in (46) to get

\mathbf{M}~\dot{\mathbf{T}}_{n+1} + \mathbf{K}~\left(\tilde{\mathbf{T}}_{n+1} + \theta~\Delta t~\dot{\mathbf{T}}_{n+1}\right) = \mathbf{f}_{n+1} ~. $$ Collect terms containing $$\dot{\mathbf{T}}_{n+1}$$ and rearrange to get
 * $$\text{(49)} \qquad

{ \left(\mathbf{M} + \theta~\Delta t~\mathbf{K}\right)~\dot{\mathbf{T}}_{n+1} = \mathbf{f}_{n+1} - \mathbf{K}~\tilde{\mathbf{T}}_{n+1} ~. } $$ We can solve equation (49) for $$\dot{\mathbf{T}}_{n+1}$$. Substitute this solution into equation (48) to get $$\mathbf{T}_{n+1}$$.

This approach is called the $$\mathbf{v}$$-form because the "velocity" or rate of change of the unknown quantity is calculated first followed by the actual quantity.

The d-form of the Generalized Midpoint Rule.
Substitute equation (48) into the first equation in (46) to get

\mathbf{M}~\left(\cfrac{\mathbf{T}_{n+1} - \tilde{\mathbf{T}}_{n+1}}{\theta~\Delta t}\right) + \mathbf{K}~\mathbf{T}_{n+1} = \mathbf{f}_{n+1} ~. $$ Collect terms containing $$\mathbf{T}_{n+1}$$ and rearrange to get
 * $$\text{(50)} \qquad

{ \left(\mathbf{M} + \theta~\Delta t~\mathbf{K}\right)\mathbf{T}_{n+1} = \theta~\Delta t~\mathbf{f}_{n+1} + \mathbf{M}~\tilde{\mathbf{T}}_{n+1}~. } $$ The system of equations (50) can be solved for $$\mathbf{T}_{n+1}$$. Once we know $$\mathbf{T}_{n+1}$$, we can calculate $$\dot{\mathbf{T}}_{n+1}$$ using equation (48).

This approach is called the $$\mathbf{d}$$-form because the "displacement" or the actual unknown quantity is calculated first followed by its rate.

Explicit foward Euler method
FE System of ODEs:

\mathbf{M}~\dot{\mathbf{T}} + \mathbf{K}~\mathbf{T} = \mathbf{f} $$

We know:

\Delta t~; \qquad \mathbf{T}_n~; \qquad \dot{\mathbf{T}}_n ~; \qquad \mathbf{f}_{n+1} $$

Apply generalized midpoint rule ($$\theta = 0$$):

\begin{align} \mathbf{M}~\dot{\mathbf{T}}_{n+1} + \mathbf{K}~\mathbf{T}_{n+1} &= \mathbf{f}_{n+1} \\ \mathbf{T}_{n+1} & = \mathbf{T}_n + \Delta t ~\dot{\mathbf{T}}_n \\ \end{align} $$ Combine:

\mathbf{M}~\dot{\mathbf{T}}_{n+1} + \mathbf{K}~\left(\mathbf{T}_n + \Delta t~\dot{\mathbf{T}}_n\right) = \mathbf{f}_{n+1} $$ Solve for $$\dot{\mathbf{T}}_{n+1}$$.

Many other techniques are available for solving systems of time-dependent ODEs. We will look into some of them later.

Quality of Approximate Solutions
The usual engineer's approach is to stop after a solution has been obtained and assume that this solution is adequate. However, it is quite important to have some information about the quality of the approximation. Unless such information is available, the finite element solution is essentially useless because it could have little resemblance to the actual solution.

 Verification is the process of determining if the numerical approximation is an accurate representation of the mathematical model. The first step in the process is to obtain a qualitative estimate of the error in the approximation. Functional analysis plays a vital role in determining these estimates of error.

The next step in the verification process is to obtain some information about whether an approximate solution  converges to the exact solution as the mesh is refined. We can also determine what the  rate of convergence or  order of accuracy for a particular approach is. We will not get into the details of error estimation in this course except for a few specific cases.

The final step in the verification process involves comparisons of numerical results with known exact solutions and experimental results of well-characterized benchmark problems.

We also need to  validate our models. Validation is the process of determining the degree to which our mathematical model represents physical reality (as far as the intended use of the model is concerned).

Later, we will discuss some aspects of verification and validation in the context of multi-physics problems.


text:
Talk:C++/Templates 1

Vitor 22:17, 14 August 2007 (UTC): You can help improving the lesson by saying what you had difficulty to understand or what you want to see in it.
 * Is it possible to come up with an example with a linked list, where the items can be different objects? --HappyCamper 03:21, 15 August 2007 (UTC)

Vitor 08:52, 15 August 2007 (UTC): Yes, thanks for asking. I will do that as soon as I have time. Thanks, HappyCamper. You mean item of different types or a list like std::list?
 * List like std::list. Maybe we can have something about the standard template library too? --HappyCamper 16:18, 15 August 2007 (UTC)
 * Vitor 16:38, 15 August 2007 (UTC): Will do. There's a topic about the STL in C++ page. Edit: Done.

"you should only define default arguments for templates with more than one type"
I disagree. This is valid:

template class Useless { public: T1 value; }; Useless<> value; /* Valid! */

--Ropez 12:58, 19 December 2007 (UTC)

Indeed, my bad. Fixed already, thanks. Vitor 12:16, 15 March 2008 (UTC)

Not well formed
Your sample program containing a List implementation is not well formed. Any name starting with _ followed by a capitol letter is reserved for implementations, as in the compiler and its library. Though you can get away with this in practice much of the time, it shouldn't be allowed in instructional material.

Hello, unknown guy. I suppose you missed the part where I said it was an immitation of the standard libraries. Vitor 12:15, 15 March 2008 (UTC)

title: Talk:C++/Templates 1


text:
Filmmaking Basics/Matte Painting/SpaceStation Intro

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto"
 * style="background-color: #ffd; border: 1px solid #777777; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px; height: 60px;" colspan="2" |

Easy for kids

 * Tux Paint
 * Tux Paint is a free painting program for kids which has special advantages for doing matte paintings.

 * Tux Paint has starter pages. These are pages that can be used by many students in a classroom situation where kids share the same computer.

 * Tux Paint has protected areas. To create protected areas, the starter page has two parts -- the protected area mask and the background picture. Kids will find it very easy to create a matte painting.

 * Mistakes are impossible
 * For this assignment, the completed part of the movie set is protected so kids will not paint over the good parts of the shot. Therefore, it is extremely easy for kids to paint a matte painting since they can never make mistakes.   Kids can paint all they want and they will not accidently cover up the protected area.

 * }

Overview

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto"
 * style="width: 100%; background-color: #eef; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" |

Last shot of "Out For A Walk"

 * Matte painting to end this scene
 * Here is the final shot of the scene called "Out For A Walk" which you saw in the lesson on film dailies.

 * The movie set is not complete. It is supposed to be a huge space station in deep space.  But really, it is just a movie set in Burbank, California.  You need to create a matte painting.

 * The Camera Movement
 * In this shot, the camera moves. Because this scene takes place in outer space, this cannot be a real camera.  Therefore we make our matte painting very large and create the effect of the camera moving.   Using a compositing program such as Adobe AfterEffects, you can use a matte painting to create the effect of zooming out into space.

 * So for the final shot, the scene starts with a shot of just the window. Then the camera moves back to show the window on the huge space ship.  Wow!

 * Therefore, the matte painting must be large -- as large as possible. The size of the matte painting determines how much of a zoom effect we get.

→ → →   This is the zoom out shot.

 * }

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto"
 * style="width: 50%; background-color: #ffffe0; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" rowspan="2" |

Quick Start

 * For this lesson, you need to
 * Down load the free (and fun) Tux Paint from http://www.tuxpaint.org/download/

 * Download the correct pair of pictures (starter) and put them into the STARTERS folder in Tux Paint.

 * Create matte painting for both scene.

 * style="width: 50%; background-color: #f0f0ff; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" |

Space Station Images
The images you use depends on the size of your computer screen. The larger, the better!

Tux Paint has a setup program which allows you to specify the window size of Tux Paint. You need to make it as large as will fit on your computer screen.
 * style="width: 50%; background-color: cornsilk; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px; height: 200px;" |
 * style="width: 50%; background-color: cornsilk; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px; height: 200px;" |

Storyboard movie set
On the next page, you can select the correct set of pictures for your computer. Click on the images to see the original. Download this pictures from Wikiversity and place the the correct folder in Tux Paint.
 * }

Basic information

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto"

 * style="width: 70%; background-color: #ffffe0; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" |

Only half a movie set

 * Half real ; half matte painting
 * A matte painting is a painting which completes a movie set. Only half of the movie set is built.  This is the half of the movie set where the actors are.  The other half (where the actors never go) is simply painted.  That is the matte painting.

 * Using a matte painting is not the same as a blue screen or green screen where all of the movie set is replaced by a painting, photograph, or 3D model. Instead, a matte painting extends the physical movie set beyond the area where the actors are.  The painting completes the movie set.

 * The second cheapest optical illusion in filmmaking
 * This means that a matte painting is an optical illusion.  You build half the movie set and then you paint the rest of the movie set.  Therefore, this illusion is very inexpensive yet the actors are able to interact with a physical movie set.

 * Example
 * Look to the right at the half-finished movie set. Part of the movie set is gorgeous.  The rest is just scrap lumber and empty space.  Your job is to replace the ugly part (the unfinished part) with a painting.  Your painting should be a space station in outer space.

The physical movie set
 * style="width: 30%; background-color: white; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px; height: 200px;" |

A sample frame
Here is a sample frame from a the video clip which is 320 by 240.

You can see the film clip on the disk from the Star Movie Shop called "24 Unedited Scenes". This scene is called "Out for a Walk".

 * }

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto" Actors in deep space
 * style="width: 100%; background-color: BlanchedAlmond; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" |

The matte painting
Here is the final matte painting. This is the movie set combined with the matte painting.
 * The complete movie set

All around the actual movie set, you see the imaginary world of the matte painting. You see still see the original movie set but now it is inside a painted movie set which, in this case, is an entire exterior of a space station.

 * }

{| cellpadding="20" cellspacing="5" style="width: 100%; background-color: inherit; margin-left: auto; margin-right: auto"
 * style="width: 100%; background-color: LemonChiffon; border: 1px solid #777777; vertical-align: top; -moz-border-radius-topleft: 8px; -moz-border-radius-bottomleft: 8px; -moz-border-radius-topright: 8px; -moz-border-radius-bottomright: 8px;" rowspan="2" |

The next page

 * For your next page, you can get a quick start using Tux Paint (a fun kids program) to create a matte painting.

 * Starter page sets pre-built for Tux Paint [[Image:Crystal Clear action forward.png|24px]]

 * }


text:
User:Ricardozuco

Ricardo Zucolotto is a 31 years-old Bachelors Degree in Information Technology. Currently living in Sydney/Australia, this professional holds approximately 10 years of experience on IT field, most of them focused in Network and Systems troubleshooting, implementations and design evolving Linux, Cisco and Microsoft infrastructure.

title: User:Ricardozuco


text:
Nonlinear finite elements/Nonlinearities in solid mechanics

Nonlinearity is natural in physical problems. In fact the linear assumptions we make are only valid in special circumstances and usually involve some measure of "smallness", for example, small strains, small displacements, small rotations, small changes in temperature, and so on.

We use linear approximations not because they are more correct but because
 * Linear solutions are easier to compute.
 * The computational cost is smaller.
 * Solutions can be superposed on each other.

However, linear analysis is not adequate and nonlinear analysis is necessary when
 * Designing high performance components.
 * Establishing the causes of failure.
 * Simulating true material behavior.
 * Trying to gain a better understanding of physical phenomena.

The choice of analysis depends on the problem at hand. The usual questions an engineer asks when deciding on the type of analysis to perform are A cost-benefit analysis is usually necessary before embarking on the nonlinear analysis of a problem.
 * What is the goal of the analysis?
 * What are the acceptable amounts of error?

Types of Nonlinearity
Nonlinearities are classified into three major types:
 * 1) Geometric nonlinearities.
 * 2) Material nonlinearities.
 * 3) Contact nonlinearities.

Geometric nonlinearities
Geometric nonlinearities involve nonlinearities in kinematic quantities such as the strain-displacement relations in solids. Such nonlinearities can occur due to large displacements, large strains, large rotations, and so on. Contact can also be classified as a geometric nonlinearity because the area of contact is a function of the deformation (some authors puts contact in another class called nonlinear boundary conditions).

An example of a problem that involves a geometric nonlinearity is the torsional spring problem (see Figure 1).

We want to find the relationship between the rotation $$\theta$$ of the bar and the applied force $$F$$.

A balance of moments on the bar gives us

M = F~l~\cos\theta $$ where $$M$$ is the moment at the pinned end of the bar.

We can relate the moment $$M$$ to the rotation of the spring by the constitutive equation for a torsional spring

M = k_T~\theta $$ where $$k_T$$ is the torsional spring constant.

Therefore, we can write

{     F = \cfrac{k_T~\theta}{l~\cos\theta} \qquad (\text{nonlinear in}~\theta)~. } $$ If the angular displacement $$\theta$$ is small (i.e., $$\theta \rightarrow 0$$), we have $$\cos\theta \rightarrow 1$$. Therefore, the force-rotation relationship becomes

{     F = \cfrac{k_T~\theta}{l} \qquad (\text{linear in}~\theta)~. } $$

Figure 2 shows how the linear approximation compares with the nonlinear solution. The nonlinear and linear solutions are quite different beyond a rotation of around 30 degrees.

Material nonlinearities
Material nonlinearities occur when the stress-strain or force-displacement law is not linear, or when material properties change with the applied loads.

For the torsional spring, we can introduce material nonlinearities if we assume that the constitutive equation for the spring is

M = (k_0 + k_1~\theta)~\theta ~. $$ In that case, the force-rotation relation becomes

{     F = \cfrac{(k_0 + k_1~\theta)~\theta}{l~\cos\theta} \qquad (\text{both material and geometric nonlinearity})~. } $$ If we assume small rotations, we get

{     F = \cfrac{(k_0 + k_1~\theta)~\theta}{l} \qquad (\text{only material nonlinearity})~. } $$ Figure 3 shows how the inclusion of material nonlinearities affects the solution.


text: File:TorsionSpring.png

title: File:TorsionSpring.png


text: File:TorsionSpringSol.png

title: File:TorsionSpringSol.png


text: File:TorsionSpringSolMGNL.png

title: File:TorsionSpringSolMGNL.png


text:
Nonlinear finite elements/Nonlinear axially loaded bar

Non-linear axial bar
Consider an axial bar of length $$L$$, the displacement of the bar is denoted by $$u$$. We assume the material of the bar to follow a non-linear constitutive rule. Figure 1. shows the axial bar with distributed body force $$q$$. If we denote the axial force in the bar as $$\sigma$$, then $$\sigma$$ depends on the deformation of the bar, this fact is consciously written as $$\sigma(u)$$. The equilibrium equation describing the axial bar is given as,

\frac{d \sigma}{dx}+q=0 $$ One should note that above equation is valid irrespective of the material behaviour of the axial bar. Taking into account that a rigid body motion do not produce axial force in the bar, on can conclude, the $$\sigma$$ can depend on the gradient of the displacement but not the displacement as such.

\sigma=f\left(\frac{du}{dx}\right) $$ In the above equation $$f$$ is the constitutive function which relates the gradient of the displacement with axial force, for a general three dimensional continua this relation will be replace by a non-linear tensorial relation. For the sake of concreteness, let us assume the $$\sigma$$ depends quadratically on the displacement gradient.

\sigma=AE\left(1+\frac{du}{dx}\right)\frac{du}{dx} $$ $$AE$$ is a constant, generally inferred as the linear axial stiffness. Using the above mentioned constitutive relation the equilibrium equation can now be written in terms of the displacements as,

\frac{d}{dx}\left(AE\left(1+\frac{du}{dx}\right)\frac{du}{dx}\right)+q=0 $$ The equilibrium equations must be supplemented with additional boundary conditions for the problem to be complete. The above equation admits two kinds of boundary conditions,

1. Dirichlet boundary $$u(x)=g$$, $$g$$ is a prescribed function defined only on the boundary $$x\in\{0,L\}$$

2. Newman boundary $$\sigma|_x=AE\left(1+\frac{du}{dx}\right)\frac{du}{dx}|_x=n$$, $$n$$ describes the traction condition of the bar at the boundary.

Although the above discussed model for a materially non-linear axial bar is really simple, it contains most of the essential features of a small deformation materially non-linear solid continua.


text: File:SigEpsNLBar.png

title: File:SigEpsNLBar.png