fluid.web. ?? ...

  • Published on
    15-Apr-2018

  • View
    228

  • Download
    15

Transcript

  • CAE ( 6 ) 2

    k LES (Large Eddy Simulation)DNS (Direct Numerical Simulation)

    1.

    2

    1.1 1 1 PrPr 1 1 (Pr = 0.7) Pr 1Pr

    1

  • Pr < 1

    Pr > 1

    Pr ~ 1

    Tw

    uT

    1: Pr

    (Pr > 1)

    xy

    Tmi Tmo

    x

    TFDHFD

    0

    2:

    1.2 2 (Hydrodynamically FullyDeveloped = HFD)Thermally Fully Developed = TFD 3 3 (qw = (T/y)|y=0 = const.) (Tw = const.) 3

    2

  • 0

    Tw Tw

    T (x, y)

    (TFD)

    (HFD)

    u (y)

    Tmi

    3:

    2.

    2.1 2 (1) (2) x y u, v

    u

    x+

    v

    y= 0 (1)

    uu

    x+ v

    u

    y= 1

    p

    x+

    (2u

    x2+

    2u

    y2

    )

    uv

    x+ v

    v

    y= 1

    p

    y+

    (2v

    x2+

    2v

    y2

    )

    (2)

    , [kg/m3] [m2/s] (3)

    [4]

    uT

    x+ v

    T

    y=

    (2T

    x2+

    2T

    y2

    )(3)

    [m2/s] u(x, y), v(x, y) T (x, y) T (1),(2)

    3

  • y

    x2Hu

    4:

    2.2 4 y = 0 2H y = H

    u

    x= 0, v = 0 (4)

    ux v (4) (4) (2)

    0 = 1

    p

    x+

    2u

    y2

    0 = 1

    p

    y

    (5)

    u(x, y) u(y), p(x, y) p(x)

    d2u

    dy2=

    1

    dp

    dx(6)

    (6) y 2

    u =1

    2

    (dp

    dx

    )y2 + C1y + C2 (7)

    y = H u = 0 (7)

    0 =1

    2

    (dp

    dx

    )H2 C1H + C2 (8)

    C1, C2

    C1 = 0

    C2 = 12

    (dp

    dx

    )H2

    (9)

    4

  • (7) u(y)

    u =1

    2

    (dp

    dx

    )(y2 H2) = 1

    2

    (dp

    dx

    )(H2 y2) (10)

    y = 0(dp/dx) > 0 y = 0 umax

    umax = u|y=0 = 12

    (dp

    dx

    )H2 (11)

    umax (10)

    u = umax

    (1 y

    2

    H2

    )(12)

    um

    um 12H

    HH

    udy

    =umax2H

    [y y

    3

    3H2

    ]H

    H=

    2

    3umax (13)

    umax um 1.5 2

    2.3 Tm T (y) y

    Tm =

    HH

    u(y)T (y)dy HH

    u(y)dy=

    HH

    u(y)T (y)dy

    2umH(14)

    um (14) cp (14) u(y) Tav =

    HH T (y)dy

    2.4

    (x, y) =Tw(x) T (x, y)Tw(x) Tm(x) (15)

    5

  • Tw, Tm(x, y) (y) (15)y

    (x, y)

    x= 0 (16)

    (15) x

    x=

    1

    (Tw Tm)2[(

    dTwdx

    Tx

    )(Tw Tm) (Tw T )

    (dTwdx

    dTmdx

    )]

    = 0 (17)

    (Tw Tm)Tx

    = (Tw Tm)dTwdx

    (Tw T )dTwdx

    + (Tw T )dTmdx

    = (T Tm)dTwdx

    + (Tw T )dTmdx

    (18)

    x T/x

    T

    x=

    T TmTw Tm

    (dTwdx

    )+

    Tw TTw Tm

    (dTmdx

    )(19)

    (19)2.6

    2.5 h qw h(Tw T) T Tw Tm h

    qw h(Tw Tm) (20)

    qw

    qw = Ty

    y=H

    = T

    y

    y=H

    (21)

    2.6 2 T (x, y)

    6

  • dqwdx

    = 0 (22)

    h (22) (20)dTwdx

    dTmdx

    = 0 (23)

    dTwdx

    =dTmdx

    (24)

    (19) T/x x

    T

    x=

    dTwdx

    =dTmdx

    = Function(x) (25)

    T/x Tw Tm x (25) 5(a)TmW x x + dxx

    mcpTm(x + dx) = mcpTm(x) + 2qwWdx (26)

    m, cp [kg/s] [J/(kgK)] (26) m =um W 2H

    dTmdx

    =qw

    cpumH(27)

    5(a)qw = const. dTm/dx = const.

    dTwdx

    = 0 (28)

    (19)T

    x=

    Tw TTw Tm

    dTmdx

    = Function(x, y) (29)

    T/xT/x x y (29) 5(b)

    7

  • 0

    x

    x

    0

    Tm

    Tw

    Tm

    Tw =const.

    TwTm = const.

    Tw

    Tmi

    Tmi

    TFD(Thermally Fully Developed)

    TFD

    (a)

    (b)

    5:

    x 6 6 Tw T Tw

    2.72.6 (25) (3)

    2T

    x2

    2T

    y2(30)

    v = 0 (3)

    uT

    x=

    2T

    y2(31)

    (25)T/x = dTm/dx (31)2.6T/x

    8

  • qw

    qw

    qw

    qw

    Tw Tw

    Tw Tw

    x

    y

    x

    y

    (a) (qw = const.)

    (b) (Tw = const.)

    6:

    2T/x2 = 0 (31) (31) (12), (13)

    2T

    y2=

    3um2

    (dTmdx

    ) (1 y

    2

    H2

    )(32)

    (32)Ax

    A 3um2

    (dTmdx

    )(33)

    y

    T

    y= A

    (y y

    3

    3H2

    )+ C1 (34)

    T = A

    (y2

    2 y

    4

    12H2

    )+ C1y + C2 (35)

    C1, C2

    y = 0 Ty

    = 0

    y = H T = Tw

    (36)

    1 y = 0 2y = H T Tw (36) (35)

    9

  • C1 = 0

    Tw = A

    (H2

    2 H

    4

    12H2

    )+ C2 = A

    5H2

    12+ C2

    (37)

    T (x, y)

    T (x, y) = A

    (y2

    2 y

    4

    12H2 5H

    2

    12

    )+ Tw(x)

    = Tw(x) AH2

    12

    [(y

    H

    )4 6

    (y

    H

    )2+ 5

    ](38)

    qw

    qw = T

    y

    y=H

    =2AH

    3(39)

    Tm(x) (14)

    Tm(x) =1

    2umH

    HH

    u Tdy

    =1

    2umH

    HH

    3um2

    (1 y

    2

    H2

    ){

    Tw(x) AH2

    12

    [(y

    H

    )4 6

    (y

    H

    )2+ 5

    ]}dy (40)

    (40) y/H dy = Hd (40)

    Tm(x) =3

    4

    11

    (1 2

    ) [Tw(x) AH

    2

    12

    (4 62 + 5

    )]d

    =3Tw(x)

    4

    11

    (1 2

    )d +

    AH2

    16

    11

    (2 1)(4 62 + 5

    )d

    = Tw(x) 34AH2

    105(41)

    (20), (39), (41) h

    h =qw

    Tw Tm =2AH/3

    34AH2/105=

    35

    17

    H(42)

    Nu 2H

    Nu h (2H)

    =35

    17H 2H

    =

    70

    17= 4.118 (43)

    10

  • 102

    110 10 10 1

    10

    (HFD, TFD)

    4.36

    3.66

    x / (D Re Pr)

    Nu

    3 2 1

    7:

    (W )Dh

    Dh = 4 = limW

    2WH 42(W + 2H)

    = 4H (44)

    Nu = 8.24NuNu = 4.12

    2.8 7 7 xDGz (Re Pr D)/xGz 7Gz1 > 0.05 4.36 3.66

    11

  • yx

    2H

    0

    FlowInlet

    8:

    3.

    Nu = 4.12

    3.12.1u, v, p,

    T

    u =u

    Uin, v =

    v

    Uin, p =

    p

    U2in, =

    (T Tin)Tin

    (45)

    Uin, Tin 8 4 8 x, y t

    x =x

    2H, y =

    y

    2H, t =

    t

    (2H/Uin)(46)

    u

    t+ u

    u

    x+ v

    u

    y= p

    x+

    1

    Re

    (2u

    x2+

    2u

    y2

    )

    v

    t+ u

    v

    x+ v

    v

    y= p

    y+

    1

    Re

    (2v

    x2+

    2v

    y2

    )

    (47)

    t+ u

    x+ v

    y=

    1

    RePr

    (2

    x2+

    2

    y2

    )(48)

    Re Uin 2H/, Pr = / (48)Re Pr Pe

    12

  • 3.2 qw

    qw = Ty

    y=0

    = (

    Tin2H

    )

    y

    y=0

    (49)

    [W/(mK)] (49) qw

    qw = qw(

    2H

    Tin

    )(50)

    qw

    qw =

    y

    y=0

    = i, j=1 i, j=0y

    (51)

    (51) (i, 0)

    (i, 0) = (i, 1) + qwy (52)

    [] (i, Ny + 1)

    qw h

    qw = h(Tw Tm) (53)

    qw = Nu(w m) (54)

    Nu = h 2H/ xxNuNum

    [] (54)

    (54)

    Nu =qw

    w m (55)

    x w m

    w =i, j=0 + i, j=1

    2(56)

    13

  • (14) u (i, j) u y m

    [] m

    m =

    10

    udy

    10

    udy=

    10

    udy (57)

    3.3Nu = 4.1248 2.7 4.118 0.16%yx 9 12 7

    c SMAC (original version by Dr. T. Ushijima)c Simplified Marker and Cell methodc Solving Heat Transfer in flow between parallel walls.c

    implicit double precision(a-h,o-z)c n

    parameter(n=20, nx=n*10, ny=n)c p: c u, v: c phi: c divup: c up, vp: c psi: c the: (e.g. )

    dimension p(0:nx+1,0:ny+1)+ ,u(0:nx,0:ny+1),v(0:nx+1,0:ny)+ ,phi(0:nx+1,0:ny+1),divup(1:nx,1:ny)+ ,up(0:nx,0:ny+1),vp(0:nx+1,0:ny)+ ,psi(0:nx,0:ny+1)

    + ,the(0:nx+1,0:ny+1),the0(0:nx+1,0:ny+1)

    c looploop=20000

    c rere=50.d0

    c pr ()pr=0.7d0

    c q_w=1.d0

    14

  • c dx(=dy)dy=1.d0/dble(n)dx=dy

    c c dx=dy*2.d0c dt

    dt=0.002d0c

    dt=min(dt,0.25*dx)c ()

    dt=min(dt,0.2*re*dx*dx)write(6,*) dt = ,dtddx=1.d0/dxddy=1.d0/dyddx2=ddx*ddxddy2=ddy*ddyddt=1.d0/dt

    c initial conditionc icont=1 c

    icont=0if (icont.eq.1) then

    open(unit=9,file=fort.21+ ,form=unformatted, status=unknown)

    read(9) u,v,pclose(9)

    elsedo 131 j=0,ny

    do 132 i=0,nxu(i,j)=1.d0v(i,j)=0.d0p(i,j)=0.d0

    the(i,j)=0.d0

    132 continue131 continue

    end ifdo 133 i=0,nx

    p(i,ny+1)=0.d0u(i,ny+1)=0.d0

    the(i,ny+1)=0.d0

    133 continuedo 134 j=0,ny

    p(nx+1,j)=0.d0v(nx+1,j)=0.d0

    the(nx+1,j)=0.d0

    134 continuec

    un=0.d0uw=1.d0us=0.d0ue=0.d0vn=0.d0vw=0.d0vs=0.d0

    15

  • ve=0.d0

    c do 1000 it=1,loop

    c c boundary condition

    do 135 j=0,nyc right wall (east) or outlet

    v(nx+1,j)=v(nx,j)u(nx,j)=u(nx-1,j)p(nx+1,j)=0.0

    the(nx+1,j)=2.*the(nx,j)-the(nx-1,j)

    c left wall (west) or inletv(0,j)=vwu(0,j)=uwp(0,j)=p(1,j)

    the(0,j)=0.d0

    135 continuec u(0,ny/2+1)=uw

    do 136 i=0,nxc lower wall (south)

    u(i,0)=2.*us-u(i,1)v(i,0)=vsp(i,0)=p(i,1)

    the(i,0)=the(i,1)+q_w*dy !heat flux constant

    c upper wall (north)u(i,ny+1)=2.*un-u(i,ny)v(i,ny)=vnp(i,ny+1)=p(i,ny)

    the(i,ny+1)=the(i,ny)+q_w*dy !heat flux constant

    136 continue

    c do 801 j=1,ny

    do 802 i=1,nxcnvt=(ddx*((the(i,j)+the(i-1,j))*u(i-1,j)

    + -(the(i,j)+the(i+1,j))*u(i,j))+ +ddy*((the(i,j)+the(i,j-1))*v(i,j-1)+ -(the(i,j)+the(i,j+1))*v(i,j)))/2.d0

    dift=(ddx2*(the(i+1,j)-2.*the(i,j)+the(i-1,j))+ +ddy2*(the(i,j+1)-2.*the(i,j)+the(i,j-1)))/(pr*re)

    the0(i,j)=the(i,j)+dt*(cnvt+dift)802 continue801 continue

    do 803 j=1,nydo 804 i=1,nx

    the(i,j)=the0(i,j)804 continue803 continue

    16

  • c c predictor stepc for u_ijc (up-u)/dt=-dp/dx-duu/dx-duv/dy+(nabla)2 u

    write(6,*)updo 125 j=1,ny

    do 126 i=1,nx-1c vij=(v(i,j)+v(i,j-1)+v(i+1,j)+v(i+1,j+1))*0.25c cnvu=0.5*ddx*(u(i,j)*(u(i+1,j)-u(i-1,j))c + -abs(u(i,j))*(u(i-1,j)-2.*u(i,j)+u(i+1,j)))c + +0.5*ddy*(vij*(u(i,j+1)-u(i,j-1))c + -abs(vij)*(u(i,j-1)-2.*u(i,j)+u(i,j+1)))c cnvu:

    cnvu=ddx*((u(i+1,j)+u(i,j))**2+ -(u(i-1,j)+u(i,j))**2)/4.d0+ +ddy*((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j))+ -(u(i,j)+u(i,j-1))*(v(i,j-1)+v(i+1,j-1)))/4.d0

    fij=-ddx*(p(i+1,j)-p(i,j))-cnvu+ +ddx2*(u(i+1,j)-2.d0*u(i,j)+u(i-1,j))/re+ +ddy2*(u(i,j+1)-2.d0*u(i,j)+u(i,j-1))/re

    up(i,j)=u(i,j)+dt*fij126 continue

    c write(6,603) (up(i,j),i=1,nx-1)125 continue

    cdo 124 j=1,ny

    up(0,j)=uwup(nx,j)=up(nx-1,j)

    124 continuec up(0,ny/2+1)=uw*1.01

    do 224 i=1,nx-1up(i,0)=2.*us-up(i,1)up(i,ny+1)=2.*un-up(i,ny)

    224 continuec for v_ijc (vp-v)/dt=-dp/dy-duv/dx-dvv/dy+(nabla)2 v

    write(6,*)vpdo 122 j=1,ny-1

    do 123 i=1,nxc uij=0.25*(u(i,j)+u(i+1,j)+u(i,j+1)+u(i+1,j+1))c cnvv=0.5*ddx*(uij*(v(i+1,j)-v(i-1,j))c + -abs(uij)*(v(i-1,j)-2.*v(i,j)+v(i+1,j)))c + +0.5*ddy*(v(i,j)*(v(i,j+1)-v(i,j-1))c + -abs(v(i,j))*(v(i,j-1)-2.*v(i,j)+v(i,j+1)))c cnvv:

    cnvv=ddx*((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j))+ -(u(i-1,j+1)+u(i-1,j))*(v(i-1,j)+v(i,j)))/4.d0+ +ddy*((v(i,j+1)+v(i,j))**2+ -(v(i,j)+v(i,j-1))**2)/4.d0

    gij=-ddy*(p(i,j+1)-p(i,j))-cnvv+ +ddx2*(v(i+1,j)-2.d0*v(i,j)+v(i-1,j))/re+ +ddy2*(v(i,j+1)-2.d0*v(i,j)+v(i,j-1))/re

    vp(i,j)=v(i,j)+dt*gij123 continue

    c write(6,603) (vp(i,j),i=1,nx)122 continue

    17

  • c do 121 i=1,nx

    vp(i,0)=vsvp(i,ny)=vn

    121 continuedo 221 j=1,ny-1

    vp(0,j)=vwvp(nx+1,j)=vp(nx,j)

    221 continuec evaluate continuity

    write(6,*) evaluate continuityic=0div=0.0do 112 j=1,ny

    do 111 i=1,nxdivup(i,j)=ddx*(up(i,j)-up(i-1,j))

    + +ddy*(vp(i,j)-vp(i,j-1))div=div+divup(i,j)**2ic=ic+1

    111 continuec write(6,603) (divup(i,j),i=1,nx)112 continue

    write(6,*) sqrt(div/dble(ic))c c solve the poisson equation (nabla)2 p=(nabla)up/dt by SOR

    write(6,*) solve the poisson equation for pressurec initialisation

    do 107 i=0,nx+1do 108 j=0,ny+1

    phi(i,j)=0.d0108 continue107 continue

    eps=1.D-6c maxitrc maxitr=nx*nyC C

    maxitr=nx*ny/10c alpha=1.51.7

    alpha=1.7do 100 iter=1,maxitr

    error=0.d0do 101 j=1,ny

    do 102 i=1,nxrhs=ddt*divup(i,j)resid=ddx2*(phi(i-1,j)-2.d0*phi(i,j)+phi(i+1,j))

    + +ddy2*(phi(i,j-1)-2.d0*phi(i,j)+phi(i,j+1))+ -rhs

    den=2.d0*(ddx2+ddy2)dphi=alpha*resid/denerror=max(abs(dphi),error)phi(i,j)=phi(i,j)+dphi

    102 continue101 continue

    do 103 j=1,ny

    18

  • phi(0,j)=phi(1,j)phi(nx+1,j)=0.0

    103 continuedo 104 i=1,nx

    phi(i,0)=phi(i,1)phi(i,ny+1)=phi(i,ny)

    104 continuec

    if (error.lt.eps) goto 998100 continue998 continue

    write(6,*) iter =, iter, it, errorc pause

    if (iter.ge.maxitr) write(6,*)maximum iteration exceeded!c c corrector step

    do 150 j=1,nydo 151 i=1,nx-1

    u(i,j)=up(i,j)-dt*ddx*(phi(i+1,j)-phi(i,j))151 continue150 continue

    do 152 j=1,ny-1do 153 i=1,nx

    v(i,j)=vp(i,j)-dt*ddy*(phi(i,j+1)-phi(i,j))153 continue152 continue

    do 160 j=1,nydo 161 i=1,nx

    p(i,j)=p(i,j)+phi(i,j)161 continue160 continue

    c c check the continuity for n+1 th step

    write(6,*) check the continuty for n+1 th stepic=0div=0.0do 155 j=1,ny

    do 156 i=1,nxdivup(i,j)=ddx*(u(i,j)-u(i-1,j))+ddy*(v(i,j)-v(i,j-1))ic=ic+1div=div+divup(i,j)**2

    156 continuec write(6,603)(divup(i,j),i=1,nx)603 format(20(1X,E10.2))155 continue

    write(6,*) sqrt(div/dble(ic))

    c *********************************************c i=nxc *********************************************

    iN=nxtm=0.d0

    c um=0.d0do 170 j=1,nytm=tm+0.5d0*(u(iN,j)+u(iN-1,j))*the(iN,j)*dy

    19

  • c um=um+0.5d0*(u(iN,j)+u(iN-1,j))*dy170 continue

    c tm=tm/umtw=0.5d0*(the(iN,0)+the(iN,1))s_Nu=q_w/(tw-tm)write(6,*) s_Nu

    1000 continuec output c

    do 491 i=0,nxpsi(i,0)=0.d0do 492 j=1,ny+1

    psi(i,j)=psi(i,j-1)+0.5*dy*(u(i,j-1)+u(i,j))492 continue491 continue

    cdo 501 j=1,ny

    do 502 i=1,nxx0=dx*(dble(i)-0.5)y0=dy*(dble(j)-0.5)u0=0.5*(u(i,j)+u(i-1,j))v0=0.5*(v(i,j)+v(i,j-1))p0=p(i,j)divu0=divup(i,j)psi0=0.5*(psi(i,j)+psi(i-1,j))

    c u,v Pwrite(10,699) x0,y0,u0,v0,p0,divu0,psi0,the(i,j)

    699 format (8(1X,E12.5))

    502 continuewrite(10,*)

    501 continuedo 503 j=1,ny-1

    do 504 i=1,nx-1x0=dx*dble(i)y0=dy*dble(j)u0=0.5*(u(i,j)+u(i,j+1))v0=0.5*(v(i,j)+v(i+1,j))omega=dx*(v(i+1,j)-v(i,j))-dy*(u(i,j+1)-u(i,j))psi0=0.5*(psi(i,j)+psi(i,j+1))

    C xy u,vwrite(11,699) x0,y0,u0,v0,omega,psi0

    504 continuewrite(11,699)

    503 continuec

    write(21) u,v,pend

    20

  • 2.0

    1.5

    1.0

    0.5

    0

    u / U

    in

    1086420

    x / (2H)

    y/(2H)=0.05

    y/(2H)=0.5

    y/(2H)=0.10

    9: 1.5

    1.0

    0.8

    0.6

    0.4

    0.2

    0

    (T -

    Tin)

    /Tin

    1086420

    x / (2H)

    y/(2H)=0.05

    y/(2H)=0.5

    10: 5(a)

    21

  • 1.0

    0.8

    0.6

    0.4

    0.2

    0

    y / (

    2H)

    1086420x / (2H)

    1.4 1.4

    1.4

    1.4 1.4 1.4

    1.3 1.3

    1.3 1.3 1.2

    1.2

    1.1

    1.1 1.1

    1.1 1.1

    1 1

    1

    0.9 0.9

    0.9

    0.8 0.8

    0.8

    0.7 0.7

    0.7

    0.6 0.6 0.6

    0.6

    0.5 0.5

    0.5

    0.4 0.4 0.4

    0.4

    0.3 0.3 0.3

    0.3

    11: u

    1.0

    0.8

    0.6

    0.4

    0.2

    0

    y / (

    2H)

    1086420x / (2H)

    0.75

    0.75

    0.7

    0.7

    0.65

    0.65

    0.6

    0.6

    0.55

    0.55

    0.55

    0.5

    0.5

    0.45

    0.45

    0.4

    0.4

    0.35

    0.35

    0.3 0.25

    0.25

    0.25

    0.2

    0.2

    0.15 0.1

    0.1

    0.0

    5

    12:

    1. , (1995),2. JSME, (2005),3. White, F. M. Heat and Mass Transfer, (1988), Addison-Wesley.4. ,43 178, (2004.1), pp. 2631.

    (http://www.htsj.or.jp/dennetsu/denpdf/2004_01.pdf)

    22