-
Notifications
You must be signed in to change notification settings - Fork 3
/
Stokes2D_ve_pureshear.jl
116 lines (114 loc) · 5.29 KB
/
Stokes2D_ve_pureshear.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
# Initialisation
using Plots, Printf, Statistics, LinearAlgebra
Dat = Float64 # Precision (double=Float64 or single=Float32)
# Macros
@views av(A) = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
@views av_xa(A) = 0.5*(A[1:end-1,:].+A[2:end,:])
@views av_ya(A) = 0.5*(A[:,1:end-1].+A[:,2:end])
# 2D Stokes routine
@views function Stokes2D_ve()
# Physics
Lx, Ly = 1.0, 1.0 # domain size
radi = 0.01 # inclusion radius
ξ = 2.0 # Maxwell relaxation time
μ0 = 1.0 # viscous viscosity
G0 = 1.0 # elastic shear modulus
Gi = 1.0/2.0 # elastic shear modulus perturbation
εbg = 1.0 # background strain-rate
# Numerics
nt = 10 # number of time steps
nx, ny = 31, 31 # numerical grid resolution
Vdmp = 4.0 # convergence acceleration (damping)
Ptsc = 8.0 # iterative time step limiter
ε = 1e-6 # nonlinear tolerence
iterMax = 1e5 # max number of iters
nout = 200 # check frequency
# Preprocessing
dx, dy = Lx/nx, Ly/ny
dt = μ0/(G0*ξ + 1e-15)
# Array initialisation
Pt = zeros(Dat, nx ,ny )
∇V = zeros(Dat, nx ,ny )
Vx = zeros(Dat, nx+1,ny )
Vy = zeros(Dat, nx ,ny+1)
Exx = zeros(Dat, nx ,ny )
Eyy = zeros(Dat, nx ,ny )
Exy = zeros(Dat, nx-1,ny-1)
Txx = zeros(Dat, nx ,ny )
Tyy = zeros(Dat, nx ,ny )
Txy = zeros(Dat, nx+1,ny+1)
Txx_o = zeros(Dat, nx ,ny )
Tyy_o = zeros(Dat, nx ,ny )
Txy_o = zeros(Dat, nx+1,ny+1)
Tii = zeros(Dat, nx ,ny )
Rx = zeros(Dat, nx-1,ny )
Ry = zeros(Dat, nx ,ny-1)
dVxdt = zeros(Dat, nx-1,ny )
dVydt = zeros(Dat, nx ,ny-1)
Rog = zeros(Dat, nx ,ny )
Xsi = ξ*ones(Dat, nx, ny)
Mus = μ0*ones(Dat, nx, ny)
G = G0*ones(Dat, nx, ny)
# Initialisation
xc, yc = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny)
xc, yc = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny)
xv, yv = LinRange(0.0, Lx, nx+1), LinRange(0.0, Ly, ny+1)
(Xvx,Yvx) = ([x for x=xv,y=yc], [y for x=xv,y=yc])
(Xvy,Yvy) = ([x for x=xc,y=yv], [y for x=xc,y=yv])
rad = (xc.-Lx./2).^2 .+ (yc'.-Ly./2).^2
G[rad.<radi].= Gi
Vx .= εbg.*Xvx
Vy .= .-εbg.*Yvy
dtVx = min(dx,dy)^2.0./av_xa(Mus)./4.1./2.0
dtVy = min(dx,dy)^2.0./av_ya(Mus)./4.1./2.0
dtPt = 4.1*Mus/max(nx,ny)/Ptsc
# Time loop
t=0.0; evo_t=[]; evo_Txx=[]
for it = 1:nt
iter=1; err=2*ε; err_evo1=[]; err_evo2=[]
Txx_o.=Txx; Tyy_o.=Tyy; Txy_o.=Txy
local itg
while (err>ε && iter<=iterMax)
# divergence - pressure
∇V .= diff(Vx, dims=1)./dx .+ diff(Vy, dims=2)./dy
Pt .= Pt .- dtPt.*∇V
# strain rates
Exx .= diff(Vx, dims=1)./dx .- 1.0/3.0*∇V
Eyy .= diff(Vy, dims=2)./dy .- 1.0/3.0*∇V
Exy .= 0.5.*(diff(Vx[2:end-1,:], dims=2)./dy .+ diff(Vy[:,2:end-1], dims=1)./dx)
# stresses
Xsi .= Mus./(G.*dt)
Txx .= Txx_o.*Xsi./(Xsi.+1.0) .+ 2.0.*Mus.*Exx./(Xsi.+1.0)
Tyy .= Tyy_o.*Xsi./(Xsi.+1.0) .+ 2.0.*Mus.*Eyy./(Xsi.+1.0)
Txy[2:end-1,2:end-1] .= Txy_o[2:end-1,2:end-1].*av(Xsi)./(av(Xsi).+1.0) .+ 2.0.*av(Mus).*Exy./(av(Xsi).+1.0)
Tii .= sqrt.(0.5*(Txx.^2 .+ Tyy.^2) .+ av(Txy).^2)
# velocities
Rx .= .-diff(Pt, dims=1)./dx .+ diff(Txx, dims=1)./dx .+ diff(Txy[2:end-1,:], dims=2)./dy
Ry .= .-diff(Pt, dims=2)./dy .+ diff(Tyy, dims=2)./dy .+ diff(Txy[:,2:end-1], dims=1)./dx .+ av_ya(Rog)
dVxdt .= dVxdt.*(1-Vdmp/nx) .+ Rx
dVydt .= dVydt.*(1-Vdmp/ny) .+ Ry
Vx[2:end-1,:] .= Vx[2:end-1,:] .+ dVxdt.*dtVx
Vy[:,2:end-1] .= Vy[:,2:end-1] .+ dVydt.*dtVy
# convergence check
if mod(iter, nout)==0
norm_Rx = norm(Rx)/length(Rx); norm_Ry = norm(Ry)/length(Ry); norm_∇V = norm(∇V)/length(∇V)
err = maximum([norm_Rx, norm_Ry, norm_∇V])
push!(err_evo1, err); push!(err_evo2, itg)
@printf("it = %d, iter = %d, err = %1.3e norm[Rx=%1.3e, Ry=%1.3e, ∇V=%1.3e] \n", it, itg, err, norm_Rx, norm_Ry, norm_∇V)
end
iter+=1; itg=iter
end
t = t + dt
push!(evo_t, t); push!(evo_Txx, maximum(Txx))
# Plotting
p1 = heatmap(xv, yc, Vx' , aspect_ratio=1, xlims=(0, Lx), ylims=(dy/2, Ly-dy/2), c=:inferno, title="Vx")
# p2 = heatmap(xc, yv, Vy' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="Vy")
p2 = heatmap(xc, yc, G' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="G")
p3 = heatmap(xc, yc, Tii' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="τii")
p4 = plot(evo_t, evo_Txx , legend=false, xlabel="time", ylabel="max(τxx)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
plot!(evo_t, 2.0.*εbg.*μ0.*(1.0.-exp.(.-evo_t.*G0./μ0)), linewidth=2.0) # analytica solution
display(plot(p1, p2, p3, p4))
end
return
end
Stokes2D_ve()