1 Atmospheric flow
Consider the steady-state form of the linearized equations in [1] restricted to 2D (that is, no v velocity component)
∂u ∂π
U ∂x = − ∂x
∂w ∂π
U ∂x = − ∂z + b
U ∂b + N 2w = g q
∂x ∂u ∂w + ∂x ∂z
cpT0 = 0
Here we have assumed a nonrotating, inviscid, Boussinesq airflow system and have linearised around a base state of an adiabatic hydrostatic atmosphere with uniform horizontal velocity U . The forcing term here is the diabatic heating rate q = q(x, z), which has units of power per unit mass. However, unlike [1], we have not made the hydrostatic assumption and so in the vertical momentum equation we have now included the inertial term U∂w/∂x.
If we use the continuity equation to replace ∂u in the first equation by − ∂w , differentiate the first equation
with respect to z and the second equation with respect to x and subtract, we notice that the common ∂2π/∂x∂z
This then becomes a single equation to determine the vertical velocity perturbation w. This could be done in FEniCS. To complete the problem we need some boundary conditions on w. Perhaps the simplest boundary condition would be to suppose w = 0 on all of the boundary of the computational region shown in the figure below. On the lower boundary, this seems an obvious boundary condition. If we think of the forcing term q as being localised to a small region within the computational domain, then it seems plausible that the resulting w will also be localised. So if the region where q is non zero is small compared to L and is a long way from the right and left hand boundaries , then setting w = 0 on these boundaries seems reasonable. The boundary condition on the upper boundary seems a little less clear. Setting w = 0 again seems to suggest that the effect of the forcing will only be slight at the upper boundary if the region where q is nonzero is small compared to H and a long way from the upper boundary.
Having found w = w(x, z), we can then use the continuity equation (the fourth equation) to find u, and the energy equation (third equation) to find b. These equations only involve the first partial derivative of u and b with respect to x, so only a single integration with respect to x is required, starting from the upwind boundary. Again, if we assume that the disturbance is only local, we can use u = 0 and b = 0 as the starting values for each of these integrations at the upwind boundary. (This is saying that the forcing will have little effect upwind.)
Tasks
Select the geometric parameters, L and W . I would suggest something like L = 200km and W = For this project, since the region is just a rectangle, you will not need to use gmsh to generate the geometry and the mesh, you can instead simply use the in-built FEniCS function RectangleMesh(....). Use a rather fine mesh, say 694 x 64
Select some values for the base state parameters, the wind speed U and N . Now
where g/cp is the adiabatic temperature gradient (lapse rate). So setting N = 0 corresponds to an adiabatic base state. Maybe this would be a good case to consider.
Select a thermal forcing disturbance q. One possibility could be to choose a gaussian disturbance centred at x = x0, z = z0,q(x, z) =Q
so as to ensure that q(x, z) dxdz = Q. (The value of A can be easily found using FEniCS.) Here the integrals are taken over the entire region. Note that q and Q are expressed per unit mass, and that for the purpose of calculating A we suppose that the density is constant over the region.
You might want to consider a number of different cases of thermal forcing, corresponding to different
locations x0, y0 and different values of σ. (Choose 4 or 5 cases in total).
x0=100, z0=2 , z0=5, sigma=5, sigma=10
For each case, use FEniCS to solve the equations for w, u and b, and view the results using
Because we are not that confident about the boundary condition at the top surface z = H, it would be useful to pick one of the cases above and, say, double H to see what effect it may
By looking at the results of your case studies, draw some conclusions about the effect of the location of the and size of the thermal forcing on the flow pattern. Also comment on the effect of the height H
1.4 Calculation of u, b and π using FEniCS
Once you have calculated the vertical velocity w using FEniCS, you next need to calculate the horizontal velocity u and the buoyancy b. As noted above these are simple one-dimensional x-integrations starting from the upwind boundary where we have supposed that u = 0 and b = 0. If we want to do this in FEniCS we need to do it a little indirectly. We need to first express it in a weak, or integral formulation that FEniCS can handle.
From the continuity equation (the fourth equation at the top), we have
∂u ∂w
∂x = − ∂z .
Multiply this by an arbitrary test function φ that is zero on the upwind (left vertical) boundary and integrate over the entire computational region R:
Now, unlike we usually do, we now do not use integration by parts to transform any of these integrals, we just leave them as they are. Our weak formulation is simply that we want to find u such that the u = 0 on the upwind boundary and the above holds for all test functions φ which are zero on the upwind boundary.
We can do a similar thing to find the buoyancy perturbation b from w. In this case the energy equation (the third equation at the top) becomes
U ∂b = g q − N 2w,
∂x cpT0
which after multiplying by a test function φ and integrating gives
Again, we leave this unchanged and seek a function b such that the b = 0 on the upwind boundary and the above holds for all test functions φ which are also zero on the upwind boundary.
Finally, we might want to calculate the pressure perturbation π; to do this, we use the x-momentum equation (the first equation above) to conclude that π = −Uu + f
where f (z) depends on z alone. Since our equations only determine π up to an additive constant, we may suppose that f (0) = 0. Thus, for z = 0, we have π(x, 0) = Uu(x, 0). Now from the z-momentum equation (the second equation) we have
∂π ∂w
∂z = b − U ∂x .
So, multiplying by a test function φ which is zero on the lower boundary z = 0 and integrating over the region
R gives
Our variational formulation then becomes to find π such that π = Uu on the lower boundary z = 0 and the above equation holds for all test function φ which are zero on the lower boundary z = 0.
We would expect that at the far upwind boundary, the effect of the disturbance is negligible, so π 0 there. From (1), it would then follow that f (z) = 0, since u = 0 is our upwind bc. Thus, π Uu should be a good approximation throughout the region R.
References
Han J-Y, Baik J-J, 2008Theoretical Studies of Convectively Forced Mesoscale Flows in Three Part 1: Uniform Basic-State Flow, J. Atmos. Sci., 66, 947–965
DescriptionIn this final assignment, the students will demonstrate their ability to apply two ma
Path finding involves finding a path from A to B. Typically we want the path to have certain properties,such as being the shortest or to avoid going t
Develop a program to emulate a purchase transaction at a retail store. Thisprogram will have two classes, a LineItem class and a Transaction class. Th
1 Project 1 Introduction - the SeaPort Project series For this set of projects for the course, we wish to simulate some of the aspects of a number of
1 Project 2 Introduction - the SeaPort Project series For this set of projects for the course, we wish to simulate some of the aspects of a number of