Skip to content

Commit

Permalink
add notebook for flat free surface
Browse files Browse the repository at this point in the history
  • Loading branch information
@tduretz
tduretz committed Aug 31, 2023
1 parent 618d4cc commit dd35af1
Show file tree
Hide file tree
Showing 2 changed files with 168 additions and 13 deletions.
34 changes: 21 additions & 13 deletions examples/Wave2D_Makie_Day6.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ function MainSource()
# Time domain
c_eff = sqrt((K₀*(1+Fb_b)+4/3*G₀)/ρ₀)
Δt = min(1e10, 0.3*Δ.x/c_eff, 0.3*Δ.y/c_eff ) # Courant criteria from wavespeed
Nt = 1000
Nt = 400
Nout = 100
t = -t₀

Expand Down Expand Up @@ -166,18 +166,26 @@ function MainSource()
@.. ε̇.j.yz = 1//2*(L.j.zy)

# Stress update
@.. τ.i.xx = θs(devi...)*(ε̇.i.xx) + χs(devi...)*τ.i.xx
@.. τ.j.xx = θs(devj...)*(ε̇.j.xx) + χs(devj...)*τ.j.xx
@.. τ.i.yy = θs(devi...)*(ε̇.i.yy) + χs(devi...)*τ.i.yy
@.. τ.j.yy = θs(devj...)*(ε̇.j.yy) + χs(devj...)*τ.j.yy
@.. τ.i.zz = θs(devi...)*(ε̇.i.zz) + χs(devi...)*τ.i.zz
@.. τ.j.zz = θs(devj...)*(ε̇.j.zz) + χs(devj...)*τ.j.zz
@.. τ.i.xy = θs(devi...)*(ε̇.i.xy) + χs(devi...)*τ.i.xy
@.. τ.j.xy = θs(devj...)*(ε̇.j.xy) + χs(devj...)*τ.j.xy
@.. τ.i.xz = θs(devi...)*(ε̇.i.xz) + χs(devi...)*τ.i.xz
@.. τ.j.xz = θs(devj...)*(ε̇.j.xz) + χs(devj...)*τ.j.xz
@.. τ.i.yz = θs(devi...)*(ε̇.i.yz) + χs(devi...)*τ.i.yz
@.. τ.j.yz = θs(devj...)*(ε̇.j.yz) + χs(devj...)*τ.j.yz
# @.. τ.i.xx = θs(devi...)*(ε̇.i.xx) + χs(devi...)*τ.i.xx
# @.. τ.j.xx = θs(devj...)*(ε̇.j.xx) + χs(devj...)*τ.j.xx
# @.. τ.i.yy = θs(devi...)*(ε̇.i.yy) + χs(devi...)*τ.i.yy
# @.. τ.j.yy = θs(devj...)*(ε̇.j.yy) + χs(devj...)*τ.j.yy
# @.. τ.i.zz = θs(devi...)*(ε̇.i.zz) + χs(devi...)*τ.i.zz
# @.. τ.j.zz = θs(devj...)*(ε̇.j.zz) + χs(devj...)*τ.j.zz
# @.. τ.i.xy = θs(devi...)*(ε̇.i.xy) + χs(devi...)*τ.i.xy
# @.. τ.j.xy = θs(devj...)*(ε̇.j.xy) + χs(devj...)*τ.j.xy
# @.. τ.i.xz = θs(devi...)*(ε̇.i.xz) + χs(devi...)*τ.i.xz
# @.. τ.j.xz = θs(devj...)*(ε̇.j.xz) + χs(devj...)*τ.j.xz
# @.. τ.i.yz = θs(devi...)*(ε̇.i.yz) + χs(devi...)*τ.i.yz
# @.. τ.j.yz = θs(devj...)*(ε̇.j.yz) + χs(devj...)*τ.j.yz


map(τ, ε̇) do τ, ε̇
@.. τ = θs(devi...)*(ε̇) + χs(devi...)*τ
end
# map(τ.i, ε̇.i) do τ, ε̇
# @.. τ = θs(devj...)*(ε̇) + χs(devj...)*τ
# end

# Pressure update
@.. P.i = P.i + θb(voli...)*∇V.i + χb(voli...)*∇V0.i
Expand Down
147 changes: 147 additions & 0 deletions notebooks/FreeSurfaceFlat.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,147 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [],
"source": [
"from sympy import *"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [],
"source": [
"# Symbols\n",
"x, y = symbols('x, y')\n",
"txx, tyy, txy, P = symbols('txx, tyy, txy, P')\n",
"txx0, tyy0, txy0, P0 = symbols('txx0, tyy0, txy0, P0')\n",
"Vx, Vy, Vz = symbols('Vx, Vy, Vz')\n",
"K, G = symbols('K, G')\n",
"etab, etas = symbols('eta_b, eta_s')\n",
"Vx = Function('Vx')(x,y)\n",
"Vy = Function('Vy')(x,y)\n",
"hx, hy = symbols('h_x, h_y')\n",
"Lxx, Lyy, Lxy, Lyx = symbols('Lxx, Lyy, Lxy, Lyx')"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\eta_{b} \\left(\\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}\\right) + 2 \\eta_{s} \\left(\\frac{2 \\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)}}{3} - \\frac{\\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}}{3}\\right) & 2 \\eta_{s} \\left(0.5 \\frac{\\partial}{\\partial y} \\operatorname{Vx}{\\left(x,y \\right)} + 0.5 \\frac{\\partial}{\\partial x} \\operatorname{Vy}{\\left(x,y \\right)}\\right)\\\\2 \\eta_{s} \\left(0.5 \\frac{\\partial}{\\partial y} \\operatorname{Vx}{\\left(x,y \\right)} + 0.5 \\frac{\\partial}{\\partial x} \\operatorname{Vy}{\\left(x,y \\right)}\\right) & \\eta_{b} \\left(\\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}\\right) + 2 \\eta_{s} \\left(- \\frac{\\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)}}{3} + \\frac{2 \\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}}{3}\\right)\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[eta_b*(Derivative(Vx(x, y), x) + Derivative(Vy(x, y), y)) + 2*eta_s*(2*Derivative(Vx(x, y), x)/3 - Derivative(Vy(x, y), y)/3), 2*eta_s*(0.5*Derivative(Vx(x, y), y) + 0.5*Derivative(Vy(x, y), x))],\n",
"[ 2*eta_s*(0.5*Derivative(Vx(x, y), y) + 0.5*Derivative(Vy(x, y), x)), eta_b*(Derivative(Vx(x, y), x) + Derivative(Vy(x, y), y)) + 2*eta_s*(-Derivative(Vx(x, y), x)/3 + 2*Derivative(Vy(x, y), y)/3)]])"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Total stress stensor\n",
"divV = Vx.diff(x) + Vy.diff(y)\n",
"Exx = Vx.diff(x) - Rational(1,3)*divV\n",
"Eyy = Vy.diff(y) - Rational(1,3)*divV\n",
"Exy = 1/2*(Vx.diff(y) + Vy.diff(x))\n",
"P = -etab*divV\n",
"E = Matrix([[Exx, Exy],[Exy,Eyy]])\n",
"D = Matrix([[2*etas, 0],[0, 2*etas]])\n",
"Tau = D*E\n",
"S = -P*eye(2) + Tau\n",
"S"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}2 \\eta_{s} \\left(0.5 \\frac{\\partial}{\\partial y} \\operatorname{Vx}{\\left(x,y \\right)} + 0.5 \\frac{\\partial}{\\partial x} \\operatorname{Vy}{\\left(x,y \\right)}\\right)\\\\\\eta_{b} \\left(\\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)} + \\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}\\right) + 2 \\eta_{s} \\left(- \\frac{\\frac{\\partial}{\\partial x} \\operatorname{Vx}{\\left(x,y \\right)}}{3} + \\frac{2 \\frac{\\partial}{\\partial y} \\operatorname{Vy}{\\left(x,y \\right)}}{3}\\right)\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 2*eta_s*(0.5*Derivative(Vx(x, y), y) + 0.5*Derivative(Vy(x, y), x))],\n",
"[eta_b*(Derivative(Vx(x, y), x) + Derivative(Vy(x, y), y)) + 2*eta_s*(-Derivative(Vx(x, y), x)/3 + 2*Derivative(Vy(x, y), y)/3)]])"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Traction for flat case\n",
"n = Matrix([[0],[1]])\n",
"T = S*n\n",
"T"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Lxy = -Lyx\n",
"Lyy = (-3.0 * Lxx .* eta_b + 2.0 * Lxx .* eta_s) ./ (3.0 * eta_b + 4.0 * eta_s)\n"
]
}
],
"source": [
"# Solve for velocity gradient components to apply as BC\n",
"T = T.subs(Vx.diff(x), Lxx).subs(Vy.diff(y), Lyy)\n",
"T = T.subs(Vx.diff(y), Lxy).subs(Vy.diff(x), Lyx)\n",
"x = solve(T, (Lxy, Lyy))\n",
"print('Lxy = '+ julia_code(x[Lxy]))\n",
"print('Lyy = '+ julia_code(x[Lyy])) "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.9.6"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}

0 comments on commit dd35af1

Please sign in to comment.