| 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 3b 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 |
|
| 16 |
|
| 17 |
nel = 20 |
| 18 |
H = 0.5 |
| 19 |
L = 1.0 |
| 20 |
|
| 21 |
eta = 1.0 # shear viscosity |
| 22 |
ro = 1.0 |
| 23 |
g = 1.00 |
| 24 |
|
| 25 |
alpha_w = 1.00 |
| 26 |
alpha = 1.00*1000000. |
| 27 |
Pen=0. |
| 28 |
B=100. |
| 29 |
|
| 30 |
nstep = 3000 |
| 31 |
dt = 1. |
| 32 |
small = EPSILON |
| 33 |
w_step=max(int(nstep/50),1)*0+1 |
| 34 |
toler = 0.001 |
| 35 |
teta1 = 0.5 |
| 36 |
teta2 = 0.5 |
| 37 |
teta3 = 1 # =0 split A; =1 split B |
| 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) |
| 45 |
momentumStep1.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) # fix x0=0 and x1=0 |
| 46 |
face_mask=whereZero(FunctionOnBoundary(dom).getX()[1]) |
| 47 |
|
| 48 |
pressureStep2=LinearSinglePDE(dom) |
| 49 |
pressureStep2.setReducedOrderOn() |
| 50 |
pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H)) |
| 51 |
|
| 52 |
momentumStep3=LinearPDESystem(dom) |
| 53 |
momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.]) |
| 54 |
# |
| 55 |
# initial values: |
| 56 |
# |
| 57 |
U=Vector(0.,Solution(dom)) |
| 58 |
p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3 |
| 59 |
p=ro*g*(H-ReducedSolution(dom).getX()[1]) |
| 60 |
dev_stress=Tensor(0.,Function(dom)) |
| 61 |
|
| 62 |
t=dt |
| 63 |
istep=0 |
| 64 |
while istep < nstep: |
| 65 |
istep=istep+1 |
| 66 |
print "time step :",istep," t = ",t |
| 67 |
r=Function(dom).getX()[0] |
| 68 |
r_b=FunctionOnBoundary(dom).getX()[0] |
| 69 |
print " volume : ",integrate(r) |
| 70 |
# |
| 71 |
# step 1: |
| 72 |
# |
| 73 |
# calculate normal |
| 74 |
n_d=dom.getNormal() |
| 75 |
t_d=matrixmult(numarray.array([[0.,-1.],[1.,0]]),n_d) |
| 76 |
sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.) |
| 77 |
print " sigma_d =",inf(sigma_d),sup(sigma_d) |
| 78 |
|
| 79 |
momentumStep1.setValue(D=r*ro*kronecker(dom), |
| 80 |
Y=r*ro*U+dt*r*[0.,-ro*g], |
| 81 |
X=-dt*r*(dev_stress-teta3*p*kronecker(dom)), |
| 82 |
y=sigma_d*face_mask*r_b) |
| 83 |
U_star=momentumStep1.getSolution() |
| 84 |
saveVTK("u.xml",u=U_star,u0=U) |
| 85 |
1/0 |
| 86 |
# |
| 87 |
# step 2: |
| 88 |
# |
| 89 |
# U2=U+teta1*(U_star-U) |
| 90 |
U2=U+teta1*U_star |
| 91 |
gg2=grad(U2) |
| 92 |
div_U2=gg2[0,0]+gg2[1,1]+U2[0]/r |
| 93 |
|
| 94 |
grad_p=grad(p) |
| 95 |
|
| 96 |
pressureStep2.setValue(A=r*dt*B*teta1*teta2/ro*dt*kronecker(dom), |
| 97 |
D=r, |
| 98 |
Y=-dt*B*r*div_U2, |
| 99 |
X=-r*B*dt**2/ro*teta1*(1-teta3)*grad_p) |
| 100 |
dp=pressureStep2.getSolution() |
| 101 |
# |
| 102 |
# step 3: |
| 103 |
# |
| 104 |
p2=(1-teta3)*p+teta2*dp |
| 105 |
grad_p2=grad(p2) |
| 106 |
momentumStep3.setValue(D=r*ro*kronecker(dom), |
| 107 |
Y=r*(ro*U_star-dt*teta2*grad_p2)) |
| 108 |
U_new=momentumStep3.getSolution() |
| 109 |
# |
| 110 |
# update: |
| 111 |
# |
| 112 |
p+=dp |
| 113 |
U=U_new |
| 114 |
print " U:",inf(U),sup(U) |
| 115 |
print " P:",inf(p),sup(p) |
| 116 |
|
| 117 |
|
| 118 |
p_pos=clip(p,small) |
| 119 |
gg=grad(U) |
| 120 |
vol=gg[0,0]+gg[1,1]+U[0]/r |
| 121 |
gamma=sqrt(2*((gg[0,0]-vol/3)**2+(gg[1,1]-vol/3)**2+(U[0]/r-vol/3)**2+(gg[1,0]+gg[0,1])**2/2)) |
| 122 |
m=whereNegative(eta*gamma-alpha*p_pos) |
| 123 |
eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small) |
| 124 |
print " viscosity =",inf(eta_d),sup(eta_d) |
| 125 |
dev_stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom)) |
| 126 |
# |
| 127 |
# step size control: |
| 128 |
# |
| 129 |
len=inf(dom.getSize()) |
| 130 |
dt1=inf(dom.getSize()/(length(U)+small)) |
| 131 |
dt2=inf(0.5*ro*(len**2)/eta_d) |
| 132 |
dt=dt1*dt2/(dt1+dt2) |
| 133 |
print " new step size = ",dt |
| 134 |
# |
| 135 |
# update geometry |
| 136 |
# |
| 137 |
dom.setX(dom.getX()+U*dt) |
| 138 |
t=t+dt |
| 139 |
if (istep-1)%w_step==0:saveVTK("u.%d.xml"%((istep-1)/w_step),p=p,eta=eta_d,U=U_star,U_star=U_star,gamma=gamma) |
| 140 |
if istep == 3: 1/0 |
Parent Directory
| 
Revision Log