{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example of a model with a compressible material\n",
"\n",
"An instantaneous 2D extension model. A compressible material is subject to lateral extension boundary conditions along the vertical walls.\n",
"\n",
"The compressible stokes flow equations are based on the compressible elasticity formulation: see Hughes, sec4.3, The Finite Element Method, 1987\n",
"\n",
"-----\n",
"The momentum equtaion\n",
"\n",
"$\n",
" \\sigma_{i,j} = -\\mathbf{p}\\delta_{i,j} + 2 \\eta \\dot\\epsilon_{i,j} = f_{i}\n",
"$\n",
"\n",
"The continuity equation\n",
"\n",
"$\n",
" \\mathbf{v}_{i,i} + \\frac{\\mathbf{p}}{\\lambda} = 0\n",
"$ \n",
"where:\n",
" * $ \\dot \\epsilon_{i,j} = \\frac{1}{2}\\left[ \\mathbf{v}_{i,j} + \\mathbf{v}_{j,i} \\right ]$\n",
" * $\\mathbf{v}$ is the velocity field\n",
" * $\\mathbf{p}$ is the pressure like variable\n",
" * $\\eta$ is the isotropic shear viscosity\n",
" * $\\lambda$ is the bulk viscosity\n",
" * $f$ is the body force"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import underworld as uw\n",
"import numpy as np\n",
"import math\n",
"import os\n",
"from underworld import function as fn\n",
"import glucifer"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"minX = 0.0\n",
"maxX = 2.0\n",
"maxY = 1.0\n",
"resX = 64\n",
"resY = 32\n",
"elementType=\"Q1/dQ0\"\n",
"\n",
"mesh = uw.mesh.FeMesh_Cartesian( elementType = (elementType), \n",
" elementRes = ( resX, resY), \n",
" minCoord = ( minX, 0.), \n",
" maxCoord = ( maxX, maxY),\n",
" periodic = [False, False] ) \n",
"\n",
"vField = uw.mesh.MeshVariable( mesh=mesh, nodeDofCount=mesh.dim )\n",
"pField = uw.mesh.MeshVariable( mesh=mesh.subMesh, nodeDofCount=1 )\n",
"\n",
"vField.data[:] = [0.,0.]\n",
"pField.data[:] = 0."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# model parameters\n",
"viscosityFn = 1.0 # isoviscous\n",
"vel_extend = 0.5 # simple extension velocity\n",
"oneonlambdaFn = 1.0e3 # 1/(bulk viscosity)\n",
"buoyancyFn = ( 0.0, 0.0 ) # the body force"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# first define strain rate tensor\n",
"strainRateFn = fn.tensor.symmetric( vField.fn_gradient )\n",
"strainRate_2ndInvariantFn = fn.tensor.second_invariant(strainRateFn)\n",
"velmag = fn.math.sqrt( fn.math.dot(vField, vField) )"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"leftWall = mesh.specialSets[\"MinI_VertexSet\"]\n",
"rightWall = mesh.specialSets[\"MaxI_VertexSet\"]\n",
"\n",
"bottomWall = mesh.specialSets[\"MinJ_VertexSet\"]\n",
"topWall = mesh.specialSets[\"MaxJ_VertexSet\"]\n",
"\n",
"iWalls = leftWall + rightWall\n",
"jWalls = bottomWall + topWall"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def buildVelocityField( mesh, velField, extV ):\n",
" '''\n",
" Build the extension velocity field. Only extending the incomp beam.\n",
" '''\n",
" import math\n",
" # set the all nodes on the vertical wall to extend extension\n",
" for index in leftWall: # velocity to the left\n",
" ycoord = mesh.data[index][1]\n",
" velField.data[index] = [-1.0*extV, 0.]\n",
"# velField.data[index] = [-1.0*extV*math.sin(ycoord*2*math.pi), 0.]\n",
" for index in rightWall: # velocity to the right\n",
" ycoord = mesh.data[index][1]\n",
" velField.data[index] = [extV, 0.]\n",
"\n",
"buildVelocityField( mesh, vField, vel_extend)\n",
"bcs_1 = uw.conditions.DirichletCondition( variable = vField, \n",
" indexSetsPerDof = ( iWalls, jWalls) ) "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# fn_bodyforce is 0. So ONLY dynamic pressure is produced in this model.\n",
"# For the incompressible material d. pressure is the mesh variable\n",
"# For the compressibly material d. pressure is -lambda*div(vField)\n",
"stokes = uw.systems.Stokes( velocityField = vField, \n",
" pressureField = pField, \n",
" conditions = bcs_1,\n",
" fn_viscosity = viscosityFn,\n",
" fn_bodyforce = buoyancyFn,\n",
" fn_one_on_lambda = oneonlambdaFn )\n",
"\n",
"solver = uw.systems.Solver( stokes )"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# can't use the AugmentedLagrangian with the penaly method yet\n",
"# solver.set_penalty(1.0e6) \n",
"solver.solve( nonLinearIterate=False, nonLinearTolerance=1e-2 )\n",
"# solver.print_stats()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# analytics\n",
"gradV = vField.fn_gradient\n",
"divV = gradV[0] + gradV[3] # du_dx + dv_dy\n",
"shouldBeZero = fn.math.abs(divV + pField*oneonlambdaFn)\n",
"\n",
"errorInt = uw.utils.Integral( shouldBeZero, mesh)\n",
"volInt = uw.utils.Integral(1.0, mesh)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.52124312125e-15\n"
]
}
],
"source": [
"# error across domain\n",
"tol = 1e-6\n",
"error = errorInt.evaluate()[0]/volInt.evaluate()[0]\n",
"if error > tol:\n",
" raise RuntimeError(\"Error: The continuity equation isn't solving within a volume averaged\" +\n",
" \" tolerance of {} - it's value is {}\".format(tol, error))\n",
"print(error)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(1.1477448914757307e-14, 1.6428302536547724e-16)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"stokes.eqResiduals"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"figVel = glucifer.Figure()\n",
"figVel.append( glucifer.objects.VectorArrows(mesh, vField, scaling=.25, arrowHead=0.2) )\n",
"figVel.show()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# # The pressure using the compute field\n",
"# figP = glucifer.Figure( **figVel )\n",
"# figP.append( glucifer.objects.Surface(mesh, pField*oneonlambdaFn, onMesh=True) )\n",
"# figP.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}