Skip to content

Instantly share code, notes, and snippets.

@computational-sediment-hyd

computational-sediment-hyd/1D Unsteady Flow using jupyter.ipynb

Last active October 6, 2017 16:15
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save computational-sediment-hyd/2ff266ef913bcbf3a142a39991e9df81 to your computer and use it in GitHub Desktop.
Save computational-sediment-hyd/2ff266ef913bcbf3a142a39991e9df81 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
1D Unsteady Flow using jupyter.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib.animation as animation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Continuous Equation\n",
"\n",
"$$\n",
" \\begin{align}\n",
" \\frac{\\partial A}{\\partial t}+\\frac{\\partial Q}{\\partial x} = 0\n",
" \\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def conEq(A, Q, dt,dx, hmin):\n",
" nmax=len(A)\n",
" An = np.array([0.0 for i in range(nmax)])\n",
"\n",
"# for i in range(1, nmax-1):\n",
" for i in range(nmax):\n",
" \n",
" if i == nmax-1:\n",
" Qfp = 0.0\n",
" else:\n",
" Qc, Qip = Q[i], Q[i+1]\n",
" Qfp = Qc if Qc > 0.0 and Qip > 0.0 else \\\n",
" (Qip if Qc < 0.0 and Qip < 0.0 else 0.5*Qc+0.5*Qip )\n",
"\n",
" if i == 0:\n",
" Qfm = 0.0\n",
" else:\n",
" Qc, Qim = Q[i], Q[i-1]\n",
" Qfm = Qim if Qc > 0.0 and Qim > 0.0 else \\\n",
" (Qc if Qc < 0.0 and Qim < 0.0 else 0.5*Qc+0.5*Qim )\n",
"\n",
" An[i] = A[i] - dt * (Qfp - Qfm) / dx\n",
" if An[i] < hmin : An[i] = hmin\n",
"\n",
" # An[0] = A[0] - dt * (Q[1] - 0.0) / dx\n",
" # An[nmax-1] = A[nmax-1] - dt * (0.0 - Q[nmax-2]) / dx\n",
" \n",
" return An"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Momentum Equation\n",
"\n",
"$$\n",
" \\begin{align}\n",
" \\frac{\\partial Q}{\\partial t}+\\frac{\\partial }{\\partial x}\\left(\\frac{Q^2}{A}\\right)+gA\\frac{\\partial H}{\\partial x}+gAi_e = 0\n",
" \\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def momentEq(A, Q, An, zb, dt,dx, hmin):\n",
" gravity = 9.8\n",
" cManing = 0.03\n",
" nmax=len(Q)\n",
"\n",
" Qn = np.array([0.0 for i in range(nmax)])\n",
"\n",
" for i in range(1, nmax-1):\n",
" if An[i] <= hmin*2 :\n",
" Qn[i] = 0.0\n",
" else:\n",
" Qc, Qip, Qim ,Qihp ,Qihm = Q[i], Q[i+1], Q[i-1], 0.5*(Q[i]+Q[i+1]), 0.5*(Q[i]+Q[i-1])\n",
" Ac, Aip, Aim ,Aihp ,Aihm = A[i], A[i+1], A[i-1], 0.5*(A[i]+A[i+1]), 0.5*(A[i]+A[i-1])\n",
" Hnc, Hnip, Hnim ,Hnihp ,Hnihm = An[i]+zb[i], An[i+1]+zb[i+1], An[i-1]+zb[i-1], \\\n",
" 0.5*(An[i]+zb[i]+An[i+1]+zb[i+1]), 0.5*(An[i]+zb[i]+An[i-1]+zb[i-1]) \n",
"\n",
" Q2byAfp = Qc**2/Ac if Qc > 0.0 and Qip > 0.0 else \\\n",
" (Qip**2/Aip if Qc < 0.0 and Qip < 0.0 else \\\n",
" Qihp**2/Aihp)\n",
"\n",
" Q2byAfm = Qim**2/Aim if Qc > 0.0 and Qim > 0.0 else \\\n",
" (Qc**2/Ac if Qc < 0.0 and Qim < 0.0 else \\\n",
" Qihm**2/Aihm) \n",
"\n",
" if Qc > 0.0 and Qip > 0.0 and Qim > 0.0: \n",
" Cr1, Cr2 = 0.5*(abs(Qc/Ac)+abs(Qip/Aip))*dt/dx, 0.5*(abs(Qc/Ac)+abs(Qim/Aim))*dt/dx\n",
" dHdx1,dHdx2 = (Hnip-Hnc)/dx, (Hnc-Hnim)/dx \n",
" elif Qc < 0.0 and Qip < 0.0 and Qim < 0.0:\n",
" Cr1, Cr2 = 0.5*(abs(Qc/Ac)+abs(Qim/Aim))*dt/dx, 0.5*(abs(Qc/Ac)+abs(Qip/Aip))*dt/dx\n",
" dHdx1,dHdx2 = (Hnc-Hnim)/dx, (Hnip-Hnc)/dx \n",
" else:\n",
" Cr1 = Cr2 = 0.5*(abs(0.5*(Qc/Ac+Qip/Aip))+abs(0.5*(Qc/Ac+Qim/Aim)))*dt/dx\n",
" dHdx1 = dHdx2 = (Hnihp-Hnihm)/dx\n",
"\n",
"# 1:down,2:up\n",
"# if Qc > 0.0 and Qip > 0.0 : \n",
"# Cr1 = 0.5*(abs(Qc/Ac)+abs(Qip/Aip))*dt/dx\n",
"# dHdx1 = (Hnip-Hnc)/dx \n",
"# elif Qc < 0.0 and Qip < 0.0 : \n",
"# Cr1 = 0.5*(abs(Qc/Ac)+abs(Qim/Aim))*dt/dx\n",
"# dHdx1 = (Hnc-Hnim)/dx\n",
"# else :\n",
"# Cr1 = 0.5*(abs(0.5*(Qc/Ac+Qip/Aip))+abs(0.5*(Qc/Ac+Qim/Aim)))*dt/dx\n",
"# dHdx1 = (Hnihp-Hnihm)/dx\n",
" \n",
"# if Qc > 0.0 and Qim > 0.0 : \n",
"# Cr2 = 0.5*(abs(Qc/Ac)+abs(Qim/Aim))*dt/dx\n",
"# dHdx2 = (Hnc-Hnim)/dx \n",
"# elif Qc < 0.0 and Qim < 0.0 : \n",
"# Cr2 = 0.5*(abs(Qc/Ac)+abs(Qip/Aip))*dt/dx\n",
"# dHdx2 = (Hnip-Hnc)/dx\n",
"# else :\n",
"# Cr2 = 0.5*(abs(0.5*(Qc/Ac+Qip/Aip))+abs(0.5*(Qc/Ac+Qim/Aim)))*dt/dx\n",
"# dHdx2 = (Hnihp-Hnihm)/dx\n",
" \n",
" w1, w2 = 1-Cr1, Cr2\n",
" Qn[i] = Q[i] - dt * (Q2byAfp - Q2byAfm) / dx \\\n",
" - dt * gravity * An[i] * (w1 * dHdx1 + w2 * dHdx2) \\\n",
" - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
"\n",
" i = 0\n",
" Qc, Qip ,Qihp = Q[i], Q[i+1], 0.5*(Q[i]+Q[i+1])\n",
" Ac, Aip ,Aihp = A[i], A[i+1], 0.5*(A[i]+A[i+1])\n",
" Hnc, Hnip ,Hnihp = An[i]+zb[i], An[i+1]+zb[i+1], 0.5*(An[i]+zb[i]+An[i+1]+zb[i+1]) \n",
"\n",
" Q2byAfp = Qc**2/Ac if Qc > 0.0 and Qip > 0.0 else \\\n",
" (Qip**2/Aip if Qc < 0.0 and Qip < 0.0 else \\\n",
" Qihp**2/Aihp)\n",
"\n",
"# Qn[i] = Q[i] - dt * (2*Q[i+1]*Q[i+1]/A[i+1]) / dx \\\n",
" Qn[i] = Q[i] - dt * (2*Q2byAfp) / dx \\\n",
" - dt * gravity * An[i] * ((Hnip - Hnc)/dx) \\\n",
" - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
" \n",
" i = nmax-1\n",
" Qc, Qim, Qihm = Q[i], Q[i-1], 0.5*(Q[i]+Q[i-1])\n",
" Ac, Aim, Aihm = A[i], A[i-1], 0.5*(A[i]+A[i-1])\n",
" Hnc, Hnim, Hnihm = An[i]+zb[i], An[i-1]+zb[i-1], 0.5*(An[i]+zb[i]+An[i-1]+zb[i-1]) \n",
"\n",
" Q2byAfm = Qim**2/Aim if Qc > 0.0 and Qim > 0.0 else \\\n",
" (Qc**2/Ac if Qc < 0.0 and Qim < 0.0 else \\\n",
" Qihm**2/Aihm) \n",
" \n",
" Qn[i] = Q[i] - dt * (-2*Q2byAfm) / dx \\\n",
" - dt * gravity * An[i] * ((Hnc - Hnim)/dx) \\\n",
" - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
" \n",
"# Qn[0] = 0.0 #-Qn[1]\n",
"# Qn[nmax-1] =0.0#-Qn[nmax-2]\n",
" \n",
" return Qn"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"def simulation(A, Q, zb, dt,dx, hmin):\n",
" An = conEq(A, Q, dt,dx, hmin)\n",
" Qn = momentEq(A, Q, An, zb, dt,dx, hmin)\n",
" return An, Qn"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Main routine"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"#%%timeit -r1\n",
"%matplotlib nbagg\n",
"#%matplotlib inline\n",
"\n",
"hmin = 10.0**(-4)\n",
"nmax = 100\n",
"ng = 50\n",
"hu = 1.0\n",
"hd = 0.5 #hmin\n",
"dx = 0.1\n",
"dt = 0.001\n",
"dtout= 0.1\n",
"tmax = 20.0 #5\n",
"\n",
"Q = np.array([0.0 for i in range(nmax)])\n",
"zb = np.array([0.0 for i in range(nmax)])\n",
"Lx = np.array([i*dx for i in range(nmax)])\n",
"\n",
"# Initial Condition\n",
"A = np.array([hd for i in range(nmax)])\n",
"# A[:ng] = hu\n",
"A[50] =0.1\n",
"#A[49] = 0.3\n",
"#A[51] = 0.3\n",
"\n",
"#A[30:40] = hu\n",
"#A[60:70] = hu\n",
"\n",
"#gragh \n",
"fig, ax = plt.figure(figsize=(4, 3), dpi=80), plt.subplot()\n",
"\n",
"def timeGenerater():\n",
" it, itout = 0, 1\n",
" isTiter, isOut = True, False\n",
" yield it*dt, isTiter, isOut\n",
"\n",
" while True:\n",
" it += 1\n",
" if it*dt >= itout*dtout:\n",
" itout += 1\n",
" isOut = True\n",
" else:\n",
" isOut = False\n",
"\n",
" if it*dt >= tmax: isTiter = False\n",
"\n",
" yield it*dt, isTiter, isOut\n",
"\n",
"tgen = timeGenerater() \n",
"t, isTiter, isOut = next(tgen)\n",
"\n",
"\n",
"\n",
"ims=[]\n",
"while isTiter:\n",
" A, Q = simulation(A, Q, zb, dt,dx, hmin)\n",
" t, isTiter, isOut = next(tgen)\n",
" \n",
" if isOut:\n",
"# plt.title('t = ' + str(t))\n",
"# title = ax.set_title('t = ' + str(t))\n",
" fr = ax.plot(Lx, A, \"r\")\n",
" ims.append(fr)\n",
"\n",
" \n",
"ani = animation.ArtistAnimation(fig, ims)\n",
"print(Q)\n",
"\n",
"\n",
"#plt.show()\n",
"#ani.save('out.gif', writer='imagemagick', fps=5) #出力かなり時間がかかる"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment