| 1 |
# |
| 2 |
# AXI-SYMMETRIC NEWTONIAN MODEL ; UPDATED LAGRANGIAN FORMULATION |
| 3 |
# |
| 4 |
# |
| 5 |
# step 1 rho*(v_star-v) = dt * (sigma'_ij,j-teta3*p,i+f_i) |
| 6 |
# step 2 dp=-dt*B*(v_j,j+teta1*v_star_j,j-dt*teta1*((1-teta3)*p_,jj+teta2*dp_,jj)) |
| 7 |
# step 3 rho*(v+-v) = -dt*((1-teta3)*p_,jj+teta2*dp_,jj) |
| 8 |
# step 4 p+=1/2(p+dp+abs(p+dp)) |
| 9 |
# step 4 sigma'i+_ij,j=f(v+,p+,...) |
| 10 |
# |
| 11 |
# |
| 12 |
from esys.escript import * |
| 13 |
from esys.escript.linearPDEs import LinearSinglePDE, LinearPDESystem |
| 14 |
from esys.finley import Rectangle |
| 15 |
iter = 0 |
| 16 |
nstep = 3000 |
| 17 |
w_step=max(int(nstep/50),1)*0+1 |
| 18 |
mstep = 5 |
| 19 |
H = 0.5 |
| 20 |
eta = 1.0 # shear viscosity |
| 21 |
ro = 1.0 |
| 22 |
g = 1.00 |
| 23 |
tauy0 = 0.9*ro*g*H/sqrt(3) #0.11/(3*sqrt(3)) |
| 24 |
Pen = 100 |
| 25 |
vtop = 0.01 |
| 26 |
dt = 10**(-10) |
| 27 |
small = EPSILON |
| 28 |
alpha = 1.00 |
| 29 |
nel = 20 |
| 30 |
# len = sqrt(0.5/nel) |
| 31 |
toler = 0.001 |
| 32 |
alphaw = 1.00 |
| 33 |
L = 1.0 |
| 34 |
teta1 = 0.5 |
| 35 |
teta2 = 0.5 |
| 36 |
teta3 = 0 # =0 split A; =1 split B |
| 37 |
Etau=1000000000. |
| 38 |
|
| 39 |
# create domain: |
| 40 |
dom=Rectangle(int(nel*L/min(L,H)),int(nel*H/min(L,H)),order=1, l0=L, l1=H) |
| 41 |
x=dom.getX() |
| 42 |
|
| 43 |
|
| 44 |
momentumStep1=LinearPDESystem(dom) # A momentumStep1 |
| 45 |
momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.]) |
| 46 |
face_mask=whereZero(FunctionOnBoundary(dom).getX()[1]) |
| 47 |
|
| 48 |
# b 1 U1_2={0.0} |
| 49 |
# b 1 U1_1={0,0} |
| 50 |
# b 4 U1_1={0.0} |
| 51 |
|
| 52 |
pressureStep2=LinearSinglePDE(dom) # A [reduced] pressureStep2 |
| 53 |
pressureStep2.setReducedOrderOn() # |
| 54 |
pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H)) |
| 55 |
# b free U1={0.0} ? |
| 56 |
|
| 57 |
|
| 58 |
|
| 59 |
momentumStep3=LinearPDESystem(dom) # A momentumStep3 |
| 60 |
momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]) # +whereZero(x[1])*[1.,1.]) |
| 61 |
# b 1 U1_2={0.0} |
| 62 |
# b 1 U1_1={0,0} |
| 63 |
# b 4 U1_1={0.0} |
| 64 |
|
| 65 |
|
| 66 |
# A deformation |
| 67 |
# e V700= [grad] V100 |
| 68 |
# e V901= V701+V704 +V101/(X1+small) |
| 69 |
# e V601=sqrt(2*((V701)^2+(V704)^2+(V101/(X1+small))^2+\ |
| 70 |
# (V702+V703)^2/2)) |
| 71 |
# e U1={abs(V901)/(V601+small)} |
| 72 |
|
| 73 |
|
| 74 |
|
| 75 |
|
| 76 |
# A stress |
| 77 |
# e V700= [grad] V100 |
| 78 |
# e V901= V701+V704 +V101/(X1+small) |
| 79 |
# e V601=sqrt(2*((V701-V901)^2+(V704-V901)^2+(V101/(X1+small)-V901)^2+\ |
| 80 |
# (V702+V703)^2/2)) |
| 81 |
# e V602=((eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))) |
| 82 |
# e V704=2*V602*V704-V401 |
| 83 |
# b 3 [integrated] {V704} |
| 84 |
# |
| 85 |
# A surface |
| 86 |
# b 3 [integrated] {1} |
| 87 |
|
| 88 |
|
| 89 |
U=Vector(0.,Solution(dom)) # V100=0 |
| 90 |
p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 # V401=ro*g*(L-X1)*(H-X2)/3 |
| 91 |
t=dt |
| 92 |
istep=0 |
| 93 |
En=0 |
| 94 |
|
| 95 |
|
| 96 |
while istep < nstep: |
| 97 |
istep=istep+1 |
| 98 |
print "time step :",istep," t = ",t |
| 99 |
|
| 100 |
n_d=dom.getNormal() |
| 101 |
t_d=matrixmult(numarray.array([[0.,1.],[-1.,0]]),n_d) |
| 102 |
|
| 103 |
s=-En*inner(U,n_d)*face_mask |
| 104 |
nStressn=s*wherePositive(s) |
| 105 |
u_boundary=Etau*inner(U,t_d)*face_mask |
| 106 |
m=whereNegative(u_boundary-nStressn) |
| 107 |
nStressTau=nStressn*(1.-m)+u_boundary*m |
| 108 |
print nStressn |
| 109 |
# print nStressn*n_d+nStressTau*t_d |
| 110 |
if istep == 20: |
| 111 |
# print nStressTau |
| 112 |
saveVTK("stress.vtu",s=(nStressn*n_d+nStressTau*t_d)) |
| 113 |
1/0 |
| 114 |
# print nStressTau |
| 115 |
|
| 116 |
r=dom.getX()[0] |
| 117 |
# A viscosity |
| 118 |
gg=grad(U) # V700= [grad] V100 |
| 119 |
vol=gg[0,0]+gg[1,1]+U[0]/(r+small) # V901= V701+V704 +V101/(X1+small) |
| 120 |
tau=sqrt(2*((gg[0,0]-vol/3)**2+(gg[1,1]-vol/3)**2+(U[0]/(r+small)-vol/3)**2+(gg[1,0]+gg[0,1])**2/2)) |
| 121 |
# V601=sqrt(2*((V701-V901/3)^2+(V704-V901/3)^2+(V101/(X1+small)-V901/3)^2+(V702+V703)^2/2)) |
| 122 |
m=whereNegative(eta*tau-alpha*p) # eta*V601<alpha*V401 |
| 123 |
eta_d=m*eta+(1.-m)*alpha*p/(tau+small) # |
| 124 |
# U1={(eta*V601<alpha*V401)*eta+(eta*V601>(alpha*V401-small))*(alpha*V401/(V601+small))} |
| 125 |
# viscosity 1 100 V100 200 V200 300 V300 400 V400 |
| 126 |
# solve |
| 127 |
print " viscosity =",inf(eta_d),sup(eta_d) # show V801 |
| 128 |
stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom)) |
| 129 |
# solve momentum equationStep1 |
| 130 |
# momentumStep1 1 100 V100 200 V200 300 V300 400 V400 |
| 131 |
momentumStep1.setValue(D=ro/dt*kronecker(dom), # {ro/dt}U1_i |
| 132 |
Y=stress[:,0]/(r+small)+[0.,-ro*g], # -{0.0,ro*g}_i |
| 133 |
X=-(stress+p*kronecker(dom)), # -D_j{V602}D_j{V100}_i-D_i{V602}D_j{V100}_j+D_i{(2/3)*V602}D_j{V100}_j-D_i{V401} |
| 134 |
y=-(nStressn*n_d+nStressTau*t_d)) |
| 135 |
dU2=momentumStep1.getSolution() |
| 136 |
# solve |
| 137 |
# dU2-> V100 |
| 138 |
# U -> V200 |
| 139 |
# p -> V500 |
| 140 |
|
| 141 |
# pressureStep2 1 100 V100 200 V200 300 V300 400 V500 |
| 142 |
|
| 143 |
gg2=grad(dU2) |
| 144 |
div_dU2=gg2[0,0]+gg2[1,1]+dU2[0]/(r+small) |
| 145 |
pressureStep2.setValue(A=r*teta1*teta2*dt**2*kronecker(dom), # -D_j{teta1*teta2*dt^2}D_jU1 |
| 146 |
D=r*ro, # {ro/Pen}U1 |
| 147 |
Y=-r*(ro*dt*vol+teta1*ro*dt*div_dU2)) # -{ro*dt}D_j{V200}_j-{teta1*ro*dt}D_j{V100}_j |
| 148 |
# solve |
| 149 |
dp=pressureStep2.getSolution() |
| 150 |
# dp -> V100 |
| 151 |
# dU2 -> V200 |
| 152 |
# U -> V300 |
| 153 |
# p -> V600 |
| 154 |
p+=dp # V601=V601+V101 |
| 155 |
# V701=V101 : dp -> V700 |
| 156 |
p=(p+abs(p))/2+small # V601=cut V601 |
| 157 |
# V601=V601+small |
| 158 |
print " pressure increment :",inf(dp),sup(dp) # show V601 |
| 159 |
print " pressure range :",inf(p),sup(p) # show V601 |
| 160 |
# popp |
| 161 |
# dU2 -> V100 |
| 162 |
# U -> V200 |
| 163 |
# p -> V500 |
| 164 |
# dp -> V600 |
| 165 |
|
| 166 |
# momentumStep3 1 100 V100 200 V200 300 V300 400 V600 |
| 167 |
A2=Tensor4(0.0,Function(dom)) |
| 168 |
for i in range(dom.getDim()): |
| 169 |
for j in range(dom.getDim()): |
| 170 |
A2[i,i,j,j]+=1 |
| 171 |
momentumStep3.setValue(A=Pen*A2, # -D_i{Pen}D_jU1_j |
| 172 |
D=ro/dt*kronecker(dom), # {ro/dt}U1_i |
| 173 |
Y=(ro/dt*(dU2+U)-teta2*grad(dp))) |
| 174 |
# Y=r*(ro/dt*(dU2+U)-teta2*dp*[1.,0]), # {ro/dt}{V200}_i+{ro/dt}{V100}_i |
| 175 |
# X=r*teta2*dp**kronecker(dom)) # -D_i{teta2*V401} |
| 176 |
U, U_old=momentumStep3.getSolution(), U |
| 177 |
# U -> V100 |
| 178 |
# dU2 -> V200 |
| 179 |
# U_old -> V300 |
| 180 |
# p -> V600 |
| 181 |
# dp -> V700 |
| 182 |
# V400=V300 |
| 183 |
# V300=V100 |
| 184 |
# popp |
| 185 |
# popp |
| 186 |
# U -> V100 |
| 187 |
# U_old -> V200 |
| 188 |
# p -> V400 |
| 189 |
# dp -> V500 |
| 190 |
# deformation |
| 191 |
# solve |
| 192 |
# show ratio |
| 193 |
# show V101 |
| 194 |
# popp |
| 195 |
|
| 196 |
print " new U:",inf(U),sup(U) # show v100 |
| 197 |
print " old U",inf(U_old),sup(U_old) # show v200 |
| 198 |
print " relative change:",Lsup(U_old)/Lsup(U) # test=L1 V200 |
| 199 |
# test0=L1 V100 |
| 200 |
# show test/test0 |
| 201 |
# black |
| 202 |
# arrow |
| 203 |
|
| 204 |
vmax=Lsup(U) # vmax1=abs(max V100) |
| 205 |
# vmax2=abs(min V100) |
| 206 |
# show vmax1 |
| 207 |
# show vmax2 |
| 208 |
|
| 209 |
# vmax=vmax1 |
| 210 |
# if vmax < vmax2 |
| 211 |
# vmax=vmax2 |
| 212 |
# endif |
| 213 |
print " max velocity = ",vmax # show vmax |
| 214 |
len=inf(dom.getSize()) |
| 215 |
dt1=len/vmax # Courant condition |
| 216 |
dt2=0.5*ro*(len**2)/eta |
| 217 |
dt=dt1*dt2/(dt1+dt2) |
| 218 |
En=1/dt |
| 219 |
t=t+dt |
| 220 |
print " time , time step ",t,dt # show t |
| 221 |
# show dt |
| 222 |
|
| 223 |
|
| 224 |
dom.setX(dom.getX()+U*dt) # mapm 500 |
| 225 |
# V500=V500+V100*dt |
| 226 |
# V600=V500 |
| 227 |
# mapm 500 |
| 228 |
|
| 229 |
print " volume : ",integrate(r)# vol=Out |
| 230 |
# show vol |
| 231 |
# dU -> V100 |
| 232 |
# U -> V200 |
| 233 |
# p -> V400 |
| 234 |
# dp -> V600 |
| 235 |
# clear |
| 236 |
# prim |
| 237 |
# saveVTK("u.%d.xml"%(istep),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp) |
| 238 |
if (istep-1)%w_step==0:saveVTK("u.%d.xml"%((istep-1)/w_step),p=p,eta=eta_d,U=U,dU2=dU2,tau=tau,dp=dp) |
Parent Directory
| 
Revision Log