Engineers and scientists are typically interested in how certain quantities vary within space. A scalar field is a mapping of values to positions within some coordinate system. In physical problems it is common to map values to three dimensions or fewer. For example, scalar fields may be generated from temperature measurements within a chemical reactor, pressure data over the Earth's surface, or x-ray intensity from telescope observations. Any mapping of a numerical quantity to a spatial coordinate is a scalar field.
Materials engineers interested in phase transformations typically deal with scalar fields corresponding to concentration as a function of position. This could be indicated using the symbols $c(x,y)$ where the $c$ corresponds to the field of interest and the $(x,y)$ indicating the dimensions of and position within the field.
There are three ways that functions can be combined: addition, multiplication and compostion (one function placed as an argument into another function). The product and chain rule are the result of computing differentials for products and compositions. The sum rule is stated as:
$$ {\frac {d(af+bg)}{dx}}=a{\frac {df}{dx}}+b{\frac {dg}{dx}} $$The product rule is:
$$ {\dfrac{d}{dx}}(u\cdot v)={\dfrac {du}{dx}}\cdot v+u\cdot {\dfrac {dv}{dx}} $$Assuming that $z(y)$ and $y(x)$ then the chain rule is stated as:
$$ \frac{dz}{dx}={\frac{dz}{dy}}\cdot {\frac{dy}{dx}} $$Reviewing these rules in the context of geometry will help build a better intuition for why these are the results.
Note that this function takes an $(x,y)$ position as input and returns a single scalar value.
In [ ]:
%matplotlib notebook
import sympy as sp
from sympy.plotting import plot3d
import numpy as np
import matplotlib.pyplot as plt
x, y = sp.symbols('x y')
We use plot3d() to visualize this field so that we use the height of the surface as a proxy for the magnitude of the scalar field.
In [ ]:
# Plot our scalar function over a specified range.
plot3d(sp.sin(2*x)*sp.sin(y), (x, -3, 3), (y, -2, 2));
An alternative to the three dimensional plot is to project the scalar values onto a two dimensional surface. Two common options for this are contour plots and heat maps. Rather than using a third dimension to represent the scalar values a contour plot traces single-valued lines through the domain whereas a heat-map uses colors to represent the value of the scalar field.
As previously stated, iso-lines are plotted within the domain and the color, annotations and position of the lines quantify the scalar field. Topographic maps use this representation to help the reader understand the locations and incline of mountains and valleys. matplotlib has a .contour() method that will generate contour plots. We demonstrate this next.
In [ ]:
delta = 0.025
xnp = np.arange(-3.0, 3.0, delta)
ynp = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(xnp, ynp)
Z = np.sin(2*X)*np.sin(Y)
In [ ]:
contours = 10
plt.figure()
CS = plt.contour(X, Y, Z, contours)
plt.clabel(CS, inline=1, fontsize=10)
plt.title(r'A Simple Contour Plot of Your Scalar Field')
plt.show()
Color choice requires consideration of the medium in which you distribute your data as well as the capabilities of the reader. An essay has been written on the topic and provides some guidance for preparation of figures. matplotlib has a module called cm that provides some standard color maps. User specified maps are possible.
In [ ]:
# Your code here.
Scalar fields have associated vector fields. One such vector field is known as the gradient and is indicated with the symbol $\nabla$ and is called nabla. The gradient is a type of directional derivative and is a vector quantity. In a physical problem of three dimensions the gradient can be thought of as pointing in the direction of the maximum spatial rate of change. Each basis vector magnitude is multiplied by the partial derivative of the field with respect to the basis vector's coordinate.
You can visualize the gradient operator as a vector with components:
$$\overrightarrow{\nabla} = \frac{\partial}{\partial x} \hat{i} + \frac{\partial}{\partial y} \hat{j} + \frac{\partial}{\partial z} \hat{k} $$When applied to a scalar field the result is a vector field - the gradient. Geometrically, the gradient is a vector field that "points uphill". The following illustrations are meant to convey some of that geometric intuition.
One way to visualize a vector field in a physical problem is to use small oriented arrow to indicate direction and length to indicate magnitude. These are known as quiver plots and they can be used in conjunction with heat maps to help visualize vector fields. In the example below the scalar field is sampled with two different point densities. A higher point density (100 points) for the filled contour plot and a lower point density (30 points) for the quivers.
In [ ]:
x0, x1 = (-3,3)
y0, y1 = (-2,2)
# Read the docstrings to understand why the numbers are given
# as complex quantities. use: ?np.mgrid
Y, X = np.mgrid[y0:y1:100j, x0:x1:100j]
Y1, X1 = np.mgrid[y0:y1:25j, x0:x1:25j]
Z = np.sin(2*X)*np.sin(Y)
# u and v here are the results of applying the gradient operation
# to our scalar field. Probably wise to check this in a seperate
# code block.
u = (2*np.sin(Y1)*np.cos(2*X1))
v = (np.sin(2*X1)*np.cos(Y1))
In [ ]:
fig, ax = plt.subplots()
filled_contour = plt.contourf(X,Y,Z,5)
fig.colorbar(filled_contour, ax=ax)
plt.quiver(X1,Y1,u,v, color='white')
plt.show()
If we interpret the gradient field as a measurement of the flow of some quantity through space, then the divergence of that vector field is a measurement of the accumulation of that quantity. To compute the divergence we compute the dot product of nabla and the vector field. The resulting quantity is a scalar quantity:
$$ \overrightarrow{\nabla} \cdot \overrightarrow{F} = \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z} $$Using the diffusion of mass as an example:
$$\overrightarrow{J} = -M \overrightarrow{\nabla} \cdot {\mu} $$and then computing the accumulation based on the vector field:
$$\frac{\partial X(x,t)}{\partial t} = - \overrightarrow{\nabla} \cdot \overrightarrow{J} $$where the minus sign indicates that accumulation occurs antiparallel to the gradient. You may have seen
In this section we keep track of the vector components by hand using the mathematical definitions. Writing this out we have a scalar field:
$$ F(x,y) $$the flux written as minus the gradient of the field (with $M=1$):
$$ \overrightarrow{J} = -\nabla F(x,y) = -\left(\frac{\partial F}{\partial x} \hat{i} + \frac{\partial F}{\partial y}\hat{j} \right) $$and the accumulation as minus the divergence of the flux:
$$ A = -\nabla \cdot \overrightarrow{J} = \left(\frac{\partial^2 F}{\partial x^2} + \frac{\partial^2 F}{\partial y^2}\right) $$
In [ ]:
x, y = sp.symbols('x y')
# Note NEW field function:
concentrationFunction = sp.sin(sp.pi*x)*sp.cos(sp.pi*y)
fluxX = -sp.diff(concentrationFunction,x)
fluxY = -sp.diff(concentrationFunction,y)
accumulationFunction = sp.diff(concentrationFunction,x,2) + sp.diff(concentrationFunction,y,2)
# We use lambdify to permit the functions to take arguments and vectorize the computations.
myConcentrationFunction = sp.lambdify((x,y), concentrationFunction, 'numpy')
myFluxX = sp.lambdify((x,y), fluxX, 'numpy')
myFluxY = sp.lambdify((x,y), fluxY, 'numpy')
myAccumulationFunction = sp.lambdify((x,y), accumulationFunction, 'numpy')
In [ ]:
concentrationFunction
In [ ]:
fluxX
In [ ]:
fluxY
We use our quiver plotting capability to plot:
zuv
In [ ]:
import numpy as np
x0, x1 = (-1,1)
y0, y1 = (-1,1)
plotResolution = 200
Y, X = np.mgrid[y0:y1:200j, x0:x1:200j]
# Quivers are on a seperate grid since they clutter things up.
Y1, X1 = np.mgrid[y0:y1:20j, x0:x1:20j]
Z = myConcentrationFunction(X,Y)
u = myFluxX(X1,Y1)
v = myFluxY(X1,Y1)
In [ ]:
fig, ax = plt.subplots()
from matplotlib import cm
plt.contourf(X,Y,Z,20, cmap=cm.coolwarm)
plt.colorbar()
plt.quiver(X1,Y1,u,v, color='white')
plt.show()
In [ ]:
accumulationFunction
In this example we use the contour plot to show locations of high (red) and low (blue) concentrations. From this concentration field we compute the flux and show that as a quiver plot overlaid on the contour plot. Examination of this figure confirms our physical understanding that mass flows from high to low concentration areas when the chemical potential is proportional to the concentration of species.
The actual RATE of accumulation at this particular INSTANT in time is given by the accumulationFunction. This needs to be recomputed for every increment of time to be practical. I'll show you how to do that in the upcoming lectures.
The Gibbs free energy is an auxiliary function re-cast in variables that are more appropriate for experimental observation.
$$ G = G(T,P,n_i,...) $$The flux is the divergence of the chemical potential. The chemical potential is the derivative of the Gibbs free energy (per mole, hence the switch to mole fraction) with respect to composition. Review the chain rule before you start this problem. The ideal solution model in terms of mole fraction for a two component system is:
$$ G(X_B) = (1-X_B)G_A + X_B G_B + RT(X_B \ln X_B + (1-X_B) \ln (1-X_B)) $$The ideal solution model requires that the activity of a species be proportional to the mole fraction of that species in solution and that the heat of mixing be zero.
Compute the divergence of the gradient of the derivative (w.r.t. concentration) of this function. Note that the concentration is a function of spatial coordinates. In this problem assume that the mole fraction $X$ is only a function of $x$, i.e. $X(x)$. $G_A$ and $G_B$ are constants that depend on the melting temperature and heat of fusion, $R$ is the gas constant and $T$ is the absolute temperature.
There are multiple ways to interact with Python and get at the gradient of a function. In the first instance we can use the coordinate system capabilities of sympy so that we can access the built in method .gradient(). We start by defining a coordinate system and then calling gradient on our scalar function. Scalars and vectors are objects of the coordinate system. See this page for more information on the vector module.
For the purposes of numeric computing, NumPy has functions and operators that work as you would expect for vector operations. What is presented below is for convenience, only.
In [ ]:
import sympy as sp
import sympy.vector as spv
C = spv.CoordSys3D('C')
spv.gradient(sp.sin(2*C.x)*sp.sin(C.y))
# The gradient function should return something that looks like u+v
# where u is a vector in the x direction and v is a vector in the
# y direction.
In [ ]:
# Define your example scalar field (a concentration like
# variable C(x,y).
exampleField = sp.sin(sp.pi*C.x)*sp.cos(sp.pi*C.y)
exampleField
We then use the built in sympy function gradient to compute the gradient:
In [ ]:
# Compute the gradient.
gradientOfField = -spv.gradient(exampleField)
gradientOfField
We then compute the divergence of the gradient. Note the absence of C.i and C.j in the answer indicating that these are not components of a vector. (Compare this to the last slide.)
In [ ]:
# Compute the divergence.
accumulation = -spv.divergence(gradientOfField)
accumulation