ЛБА дахь Умбуур хязгаарын арга /Immersed Boundary method in LBM/

Умбуур хязгаарын нөхцөл

Тогтоол усанд унаж доошоо живж буй бөмбөлөг

Судалж буй шингэнээр бүрэн хүрээлэгдсэн хатуу биеийн гадаргууг умбуур хязгаар гэж нэрлэнэ. Жишээ нь тэнгэрт нисэж буй онгоц агаарт умбасан байна, голын усанд зөөгдөж буй хагшаас мөн л умбасан байна, сэлж буй хүн мөн л усанд умбасан байна. Умбасан бие шингэний урсгалын улмаар хөдлөх ба хөдлөхдөө эргээд шингэний хөдөлгөөндөө нөлөөлдөг. Өөрөөр хэлбэл шингэн дотор хөдөлж байу хатуу биет болон шингэн хоорондоо харилцан үйлчлэлтэй байна. Иймд хатуу-шингэний систем гэж нэрлэнэ. Хөдөлж буй хязгаарын нөхцлийн тухай стандарт ЛБА-ын тухай  Шингэн дахь хатуу биеийн хөдөлгөөн гэсэн нийтлэлээс уншаарай. 


Умбуур хязгаарын нөхцлийн олон хувилбар тооцооны өөр өөр аргууд байгаатай адил ЛБА мөн л өөр өөр хувилбараар боловсруулж байна. Түүний нэг болох хялбар бөгөөд нарийвчлал сайтай аргыг авч үзье. 1998 онд Нобле болон Торкзински нар умбасан биетэд ЛБА-ын мөргөлдөх алхамыг хэрхэн тооцож болохыг

гэж тодорхойлсон ба энд байх betta нь нүдний хатуу хэсгийн харьцаа болон хэмжээсгүй тэнцвэрт орох хугацаанаас хамаарах жинлэх функц бөгөөд

гэж илэрхийлэгдэнэ. ЛБА-ын тэгшитгэлд байх f_i^m нь нэмэлт мөргөлдөөний гишүүн бөгөөд түгэлтийн функцын тэнвэрт бус хэсгийг буцаж ойлгох үүргийг гүйцэтгэх ба

гэж өгөгдөнө. Энд u_s нь t хугацаан дахь хатуу хэсгийн хурд ба –i нь i–ынхаа эсрэг чиглэлийг тэмдэглэхэд хэрэглэгдсэн юм.

Тэгшитгэл 3.14-ын зохион байгуулалтаас В нь хатуу хэсгийн харьцаа 0 ээс 1 хүртэл өөрчлөгдөхтэй адилаар 0 – ээс 1 гэсэн утгыг авахаар байгаа нь тодорхой харагдаж байна. Тэгшитгэл 3.13-т орж буй өөрчлөлт нь хатуу хэсгээр бүрхэгдсэн нүднүүдийн хувьд шингэний урсгалыг байхгүй болгох, хатуу хэсгийн хөдөлгөөнийг шингэний хурдтай тохируулах гэх мэт замаар түгэлтийн функцуудыг засварлах явдал юм. Энэ хязгаарын нөхцлийн чухал тал нь хатуу хэсгийн харьцаа гэдэг утгын дагуу хатуу хэсгийн хөдөлгөөний дэд торны шийдлийг бий болгодог ба энэ талаар дараа авч үзье.

Салангид элементийн аргатай хослох нь

Товчхондоо хатуу-шингэний систем нь хатуу хэсгийн нөлөөллөөр шингэний урсгал өөрчлөгдөх боловч шингэний урсгалын нөлөөллөөр хатуу хэсэг шилжилтэнд орох ёстой юм. Хөдлөх хязгаарын нөхцөл нь дээрх дурьдсан нөлөөллүүдийг шингэн доторх хөдлөх хатуу хэсгийн хувьд тооцож өгнө. Шингэний зүгээс хатуу хэсэгт үйлчлэх хүч нь хатуу хэсгээр бүрхэгдсэн n зангилааны дагуу бий болох моментын шилжилтийг нийлбэрчилсэнээр

гэж тодорхойлогдоно. Энэ хүчний тооцоолол нь тоон аргынхаа хувьд үнэн зөв. Ямарч шингэний зүгээс бий болох ойролцоолсон хүчдэл байхгүй. Харин хүчнээс үүдэлтэй мушгих хүчийг

ба энд xs нь t хугацаан дахь хатуу биетийн төв болно. Хатуу-шингэний систем нь хугацааны алхам бүрт эхлээд шингэний урсгалыг тооцож дараа нь хатуу биетийн шилжилтийг тооцох маягаар хийгдэнэ. Хугацааны агшин бүрд дараах алгоритм хүчинтэй.

Шингэний урсгалыг шийдэх

1.       t хугацаан дахь шингэний сүлжээнд давхцаж байгаа салангид элементүүдийг тодорхойлно. Х ба Ү тэнхлэг мэдэгдэх учир сүлжээн дотор байгаа элементүүдийг тодорхойлж тэдгээрийн байрлалыг мэдэхэд хэрэглэгдэнэ.

2.       Давхцан байрлаж байгаа бүх элементийн хувьд бүх сүлжээний хувьд хагас болон бүтэн бүрхэгдсэн эсэхийг тогтоож хатуу хэсгийн харьцааг олно. Мөн, зангилаанууд дээрх хатуу биеийн хөдөлгөөн дээр суурилж эдгээр хагас ба бүтэн бүрхэгдсэн зангилаанд хатуу хэсгийн шилжилтийн хурднуудыг өгнө.

3.       ЛБА-ыг тэгшитгэл 3.13-ын дагуу тооцоолно, ингэхдээ 2-т бичигдсэн мэдээллийг ашиглана. тэг бус хатуу хэсгийн харьцааны утгатай зангилаа бүрийн хувьд хатуу хэсгийн хүчийг тэгшитгэл 3.16-ийн дагуу дотоод нийлбэр авч тооцно. Эцэст нь урсах алхам хийгдэнэ. Дотоод сумма

4.       салангид элемент дээрх шингэний хүчийг салангид элемент бүрд зангилааны хүчнүүдийг нэмсэнээр тодорхойлно. Ингэдхээ 3.16 ба 3.17 -ын хувьд гадаад нийлбэрийг авна. гадаад сумма

Салангид элементийн байрлалыг тогтоох

1.       салангид элементийн аргын харилцан үйлчлэлийг тодорхойлох алгоритмыг хэрэглэж ойролцоо салангад элементүүдийг тодорхойлох

2.       шүргэлцэлүүдийг тодорхойлж зохих хүчнүүдийг өгөх

3.       харилцан үйлчлэлийн хүч, шингэний хүч, болон бүх биеийн хүчнүүдийн нийлбэрчлэлийг хэрэглэн хөдөлгөөний тэгшитгэлийг нэгтгэх. Шингэний хөвүүлэх хүч нь хүндийн хүчний хурдатгалыг живсэн хэсгийн эзлэхүүнээр үржүүлж тодорхойлогдоно. Хатуу хэсгүүдийн байрлал тодорхойлогдох ба гарсан хурд нь шингэний урсгалыг шийдэх хэсэгт хэрэглэгдэх болно.

Эхлээд зангилаа ба хатуу хэсэг хоёрын холбоог авч үзье. Биет буюу элемент нь олон зангилааг хамрах ба зангилааны нүд нь олон элементээр бүрхэгдэнэ. Хоёрт сүлжээний нүдний хатуу хэсгийн харьцааг маш нарийн тодорхойлох хэрэгтэй, дээрхийн 2 дахь нь. Ялангуяа зангилааны алхам элементийн хэмжээнээс бага үед нарийн тооцох хэрэгтэй. хэрвээ бид салангид элементийн характеристик уртыг О(mh) гэж авч үзвэл \m нь зангилааны тоо\ хоёр хэмжээстэд элемент бүр O(m2) зангилааг эзлэнэ. Тиймээс тэнд O(m2) гэсэн хатуу хэсгийн харьцаа элемент бүрт тодорхойлогдох ба n салангад элементын хувьд хугацааны алхам бүрт O(nm2) тооцоолол хийгдэнэ. Тооцоо их хугацаа, багтаамж эзлэхээр байна.

Хатуу хэсгийн харьцааны тооцоо

Зангилаа бүр дээр, хатуу хэсгийн харьцааны тооцоо нь зангилааны нүдний хоёр хэмжээстэд хатуу хэсгээр бүрхэгдсэн хувь, гурван хэмжээстэд эзлэхүүний тооцооллыг агуулна. Үр ашигтай байхын тулд зураг 3.3-т квадрат нүднүүд, элементийг харуулсан байна. Зураг 3.3-т харуулсан шиг дотроо дугуй элемент агуулах гадаад нүднүүдийг авч үзье. Нүднүүдэд багтсан байх дискийн хувьд нүд болон дискний хооронд хоёроос багагүй огтлолцол байна. Eрөнхийдөө, энэ хоёр цэг нь нүдний дөрвөн талтай дискний хүрээтэй хэрхэн огтлолцож байгааг хайснаар тодорхойлогдох болно. Нэгэнт эдгээр цэгүүд олдвол, бүрхэгдсэн талбай нь 3-аас 5 хэсэгтэй полигон болон дугуй сегмэнтэд хуваагдах ба эдгээрт үл мэдэгдэх талбайнууд харгалзана. 

Жишээ бодлого

Шингэн дотор суух хатуу бөмбөлөгийг бодсон кодны хэсгээс оруулъя. Хамгийн гол хэсэг мөргөлдөөнийг бодох хэсэгт хийгдэнэ. Энд дээр өгүүлсэн тэгшитгэл 3.13-3.15 зэргийг кодлоно. 

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

! *********************************************************************! ! ***** Шингэний урсгалд зориулсан мөргөлдөөний дэд програм *****! ! *********************************************************************! subroutine collesion include 'paramc.h' ! force on fluid! 2016.01.24 feqib=0. !initialization of temporal Dfs bb=0. ! B parameter in IB do i=0,n do j=0,m ubx=0. vbx=0. if(flag(i,j).eq.1.or.flag(i,j).eq.2) then fox(i,j)=rho(i,j)*gx foy(i,j)=0.!rho(i,j)*gy !-gbeta*((th(i,j)-tref)**2)*rho(i,j) if(lf(i,j).ne.1) then ubx=ub(1)+abs(wang(1))*abs(real(i)-xb(1)) vbx=vb(1)+abs(wang(1))*abs(real(j)-yb(1)) else ubx=0. vbx=0. endif t1=u(i,j)*u(i,j)+v(i,j)*v(i,j) ti=ubx*ubx+vbx*vbx ! moving velocity of solid t1 do k=0,8 t2=u(i,j)*cx(k)+v(i,j)*cy(k) tb=ubx*cx(k)+vbx*cy(k) ! moving velocity of solid t2 feqib(k,i,j)=w(k)*rho(i,j)*(1.+3.0*tb+4.50*tb*tb-1.50*ti) - нэмэлт мөргөлдөөний гишүүн feq(k,i,j)=w(k)*rho(i,j)*(1.+3.0*t2+4.50*t2*t2-1.50*t1) enddo endif enddo enddo DO i=0,n DO j=0,m if(flag(i,j).eq.1.or.flag(i,j).eq.2) then ! 2015.12.5 afl=lf(i,j) bb(i,j)=(1.-afl)*(tau(i,j)-0.5)/(afl+tau(i,j)-0.5) - жинлэх функц t1=u(i,j)*u(i,j)+v(i,j)*v(i,j) tf2=u(i,j)*fox(i,j)+v(i,j)*foy(i,j) DO k=0,8 op=opp(k) omega=1./tau(i,j) tf1=cx(k)*fox(i,j)+cy(k)*foy(i,j) t2=u(i,j)*cx(k)+v(i,j)*cy(k) ! force=3.*w(k)*gg*cy(k)*rho(i,j) ! force scheme-1 ! feq(k,i,j)=w(k)*rho(i,j)*(1.+3.0*t2+4.50*t2*t2-1.50*t1) force=w(k)*(1.-0.5*omega) & *(3.*tf1+9.*tf2*t2-3.*tf2) ! force scheme-2 if(i.eq.0.or.i.eq.n) force =0.0 if(j.eq.0.or.j.eq.m) force =0.0 omef(k,i,j)=f(op,i,j)-f(k,i,j)+feqib(k,i,j)-feq(op,i,j) f(k,i,j)=omega*feq(k,i,j)*(1-bb(i,j))+(1.-omega+omega*bb(i,j)) & *f(k,i,j)+force+omef(k,i,j)*bb(i,j) ! f(k,i,j)=f(k,i,j)-omega*(1.-bb(i,j))*(f(k,i,j)-feq(k,i,j))+force ! & +omef(k,i,j)*bb(i,j) ! f(k,i,j)=omega*feq(k,i,j)+(1.-omega)*f(k,i,j) if(isnan(f(k,i,j))) then print*, 'Distribution function is NaN' print*, 'failed time step is', kk, i, j stop endif END DO end if END DO END DO ! print*, 'collision has been performed' return end

Үүний дараа хатуу биеийн хөдөлгөөнийг тодорхойлох хэрэгтэй. Шингэний зүгээс хатуу биед үйлчлэх гидродинамикийн хүч, мушгих хүчийг тооцох хэрэгтэй. 

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 117 118 119 120 121 122 123 124 125 126

! *********************************************************************! ! ***** умбах биеийн хүчийг тооцож, хурд, хөдөлгөөнийг *****! ! *********************************************************************! subroutine forces include 'paramc.h' real forx_old(1), fory_old(1), torz_old(1) real forx(1),fory(1),torz(1) real uu(1), vv(1), ww(1), xx(1), yy(1), rr(1) save uu, vv, ww, xx, yy, rr, forx_old,fory_old,torz_old if(kk.eq.0) then uu(1)=0. vv(1)=0. ww(1)=0. xx(1)=xb(1)*dx yy(1)=yb(1)*dx rr(1)=rb(1)*dx forx_old=0. fory_old=0. torz_old=0. print*, 'initialized' else endif ! fluid force on object scale_f=1000.*dx**3/dt**2 !dx**3/dt**2 scale_t=1000.*dx**4/dt**2 !dx**4/dt**2 scale_u=dx/dt sum_outx=0.0 sum_outy=0.0 sum_torg=0.0 do i=0,n do j=0,m sum_innx=0.0 sum_inny=0.0 if(lf(i,j).ne.1.) then do k=1,8 sum_innx=sum_innx+omef(k,i,j)*cx(k) sum_inny=sum_inny+omef(k,i,j)*cy(k) enddo sum_outx=sum_outx+bb(i,j)*sum_innx sum_outy=sum_outy+bb(i,j)*sum_inny sum_torg=sum_torg & +(float(i)-xb(1))*bb(i,j)*sum_innx & +(float(j)-yb(1))*bb(i,j)*sum_inny endif enddo enddo forx(1)=-sum_outx*scale_f fory(1)=-sum_outy*scale_f-pi*rr(1)**2*9.81*(srho-1000.)!pi*rr(1)**2*9.81*1.65 torz(1)=-sum_torg*scale_t forx(1)=0.5*(forx_old(1)+forx(1)) fory(1)=0.5*(fory_old(1)+fory(1)) torz(1)=0.5*(torz_old(1)+torz(1)) ! force checking print open(111,file='forces.dat') write(111,*) kk, sum_outy*scale_f, fory(1), torz(1) ! old force copy forx_old(1)=forx(1) fory_old(1)=fory(1) torz_old(1)=torz(1) ! object velocity updates uu(1)=uu(1)+dt*forx(1)/(srho*pi*rr(1)**2) ! uu real velocity vv(1)=vv(1)+dt*fory(1)/(srho*pi*rr(1)**2) ! vv real velocity ww(1)=ww(1)+dt*torz(1)/(srho *rr(1)**2) ! print*, 'uu', uu(1), 'vv', vv(1) ! object position updates xx(1)=dt*uu(1)+xx(1) yy(1)=dt*vv(1)+yy(1) theta(1)=dt*ww(1)+theta(1) !think rotation need old value ! print*, 'xx', xx(1), 'yy', yy(1), 'theta', theta(1) ! scaling to lattice variables ub(1)=uu(1)/scale_u vb(1)=vv(1)/scale_u wang(1)=ww(1)*dt xb(1)=xx(1)/dx yb(1)=yy(1)/dx ! stop positions if(yb(1).le.0) then yb(1)=20 vb(1)=0. else if(xb(1).le.0.or.xb(1).ge.n) then xb(1)=n/20 ub(1)=0. endif ! restimation of object position do i=xb(1)-15,xb(1)+15!0,n do j=yb(1)-15,yb(1)+15!0,m if(flag(i,j).eq.1) then lf(i,j)=1. ! for all if ((real(i)-xb(1))**2+(real(j)-yb(1))**2 & .le.(rb(1)+rb(1)*0.095)**2) then ! flag(i,j)=2 lf(i,j)=0.5 e(i,j)=1.0 ma(i,j)=rho(i,j)*e(i,j) !0.5 endif if((real(i)-xb(1))**2+(real(j)-yb(1))**2 & .le.(rb(1)-rb(1)*0.019)**2) then ! flag(i,j)=1 ! flag 0- changed ! movflag(i,j)=1 lf(i,j)=0.0 e(i,j)=1.0 ma(i,j)=rho(i,j)*e(i,j) !1.0 endif endif enddo enddo ! solid velocity updates do i=0,n do j=0,m if(lf(i,j).ne.1.) then u(i,j)=ub(1)!0. v(i,j)=vb(1)!0. endif enddo enddo return end

Энэ кодны боловсруулалт нь зөвхөн нэг хатуу бөмбөлөгт зориулсан гэдгийг анхааруулъя. Хэрэв маш олон хатуу обьектыг тооцох бол Салангид элементийн аргатай хослох шаардлагатай болно. 

ЛБА-ын бусад үйлдлүүд хэвээр байх ба дээрх кодоор хэдэн туршилт явуулсаныг доор хуваалцая. 

Хөдөлгөөнгүй шингэн дотор чөлөөтэй тунаж буй хатуу бөмбөлөг.

Энэхүү умбуур хязгаарын аргыг шингэний чөлөөт гадаргууг тодорхойлох алгоритмтой хослуулж члөөт гадаргуу дээр хөвөх хатуу биетийг загварчлах юм. Мөн голын мөс ханзрах үед бий болох мөсөн хахааны үзэгдлийг загварчлахад нэг алхам ойртоно. Энэ арга маш амархан бөгөөд шингэн ба хатуу хэсгийн харилцан хамаарлыг төвөггүй тодорхойлж байгаа нь бусад тооцонгийн аргуудаас хамаагүй давуу харагдаж байгаа юм. Арга нь хөдөлгөөний тэгшитгэлээс бусад хэсэгтээ бүрэн параллел чанартай учир ЛБА-ын параллель тооцоололд төвөггүй орно.