# Annotation of /trunk/finley/test/python/axisymm-splitB.py

Revision 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (12 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 3952 byte(s)

 1 gross 1562 # 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 gross 1639 p=ro*g*(H-ReducedSolution(dom).getX()[1]) 60 dev_stress=Tensor(0.,Function(dom)) 61 gross 1562 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 gross 1639 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 gross 1562 U_star=momentumStep1.getSolution() 84 gross 1639 saveVTK("u.xml",u=U_star,u0=U) 85 1/0 86 gross 1562 # 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 gross 1639 dev_stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom)) 126 gross 1562 # 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

Name Value
svn:executable *