Ghostscript wrapper for D:\Digitalizacja\MTS88_t26_z1_4_PDF_artykuly\01mts88_t26_zeszyt_1.pdf M ECH AN IKA TEORETYCZNA STOSOWANA 1, 26,  1988 WYZNACZANIE  POWIERZCHNI  S WOBODNEJ  CIECZY  W  TRÓJKĄ TNYM ROWKU M AR I U SZ  KACZ M AREK IPPT   PAN , Poznań JAN .  A.  KOŁ OD Z I E J Politechnika Poznań ska G R Z E G OR Z  M U SI ELAK IPPT   PAN , Poznań W  pracy  wyznacza  się  kształ t  powierzchni  swobodnej  cieczy  znajdują cej  się  w  trój- ką tnym  rowku.  W tym  celu  korzysta  się  z  równania  Laplace'a- Younga.  Parametrami ustalonymi  są : ką t rozwarcia  rowka,  ką t zwilż enia,  liczba  Bonda.  Opisano  stosowaną numeryczną   metodę   optymalizacyjną ,  podano  przykł adowe  wyniki  oraz  program  w  ję - zyku  Basic  mikrokomputera  H P86B,  przy  pomocy  którego  wykonywano  obliczenia. 1.  Wprowadzenie W  niektórych  zagadnieniach  technicznych dą ży  się  do zwię kszenia  powierzchni swo- bodnej cieczy.  Celem  może być intensyfikacja  procesów  wymiany  ciepł y i masy  pomię dzy cieczą   i  gazem  na powierzchni  swobodnej  (parowanie  lub kondensacja).  Jedną   z  dróg zwię kszenia  powierzchni  swobodnej  jest  stosowanie  powierzchni  ż ł obionych  trójką tnymi rowkami lub powierzchni pofał dowanych, na których  zachodzi  spływ  cieczy  pod wpł ywem sił   grawitacji.  Przy  przepł ywie  wzdł uż rowków,  przy  niewielkiej  iloś ci  cieczy w  rowkach, wpływ  sił  kapilarnych  jest  istotny  i  nastę puje  znaczny  wzrost  efektywnej  powierzchni wymiany  w porównaniu  do spł ywu  po powierzchni  pł askiej. Problem  laminarnego  przepł ywu  cieczy  lepkiej,  nieś ciś liwej  w trójką tnym  rowku pod wpływem  sił  grawitacji  był   rozważ any  w  pracy  [1], U wzglę dniono  tam efekt  napię cia powierzchniowego  na powierzchni  swobodnej  cieczy, jednak  przyję to  zał oż enie upraszcza- ją ce, że powierzchnia  swobodna  ma stał y promień  krzywizny.  Zał oż enie to jest  korzystne z  punktu  widzenia  stosowanej  metody  rozwią zywania  zagadnienia  przepł ywu,  zapewnia bowiem  analityczną   postać  funkcji  okreś lają cej  kształ t powierzchni  swobodnej,  może  być jednak  przybliż eniem  niewystarczają co  dokł adnym  z  punktu  widzenia  modelowania rzeczywistej  powierzchni  swobodnej.  Ponadto  zał oż enie takie  ogranicza  budowę   modelu teoretycznego  optymalizacji  procesu  wymiany. 84  M .  KACZMARUK  I  IN N I Celem  niniejszej  pracy  jest  podanie  metody  numerycznego  wyznaczania  kształ tu po- wierzchni  swobodnej  dla  dowolnych  iloś ci  cieczy  w  trójką tnym  rowku,  przy  dowolnym, n apię ciu  powierzchniowym  i  dowolnej  gę stoś ci  cieczy,  to jest  dla  dowolnych  liczb  Bonda Z ał oż enie  stał ego  promienia  krzywizny,  jakie  przyję to  w  pracy  [1] jest  uzasadnione przy bardzo  mał ych  wartoś ciach  liczby  Bonda.  Podstawą   rozważ ań  jest  zwią zek  Laplace'a- Youn ga,  który  w  rozważ anym  przypadku  prowadzi  do  dwupunktowego  zagadnienia brzegowego  i  nieliniowym  równaniem  róż niczkowym  drugiego  rzę du,  w  którym  niewia- dom ą   wielkoś cią   jest  współ rzę dna  powierzchni  swobodnej. W  przyję tej  metodzie rozwią zywania  zagadnienie brzegowe  sprowadza  się   do zagadnie- nia  począ tkowego.  W  takim  uję ciu  istotna  róż nica w  stosunku  do  znanych  autorom roz- wią zań  tego  typu,  n p.  z  pracy  [2],  w  której  okreś lano  kształ t  powierzchni  swobodnej cieczy,  polega  n a  tym,  iż  w  niniejszej  pracy  nieznany jest  jeden  z  warunków  począ tko- wych.  Z am iast tego  zadany jest  warunek  okreś lonego  pola pod  powierzchnią   swobodną . Z  uwagi  na  duże  rozpowszechnienie  mikrokomputerów  do  pracy  doł ą czony jest  prog- ram  obliczeń  w  ję zyku  Basic. 2.  Sformułowanie  problemu Weź my  pod  uwagę   rowek  trójką tny  o  ką cie  rozwarcia  20  (rys.  1).  Przyjmijmy  bie- gunowy  ukł ad  współ rzę dnych  (r,    =  0. (6) D rugi  warun ek  brzegowy  wyn ika  z ką ta  zwilż enia  n a gran icy  trzech oś rodków  i m a  p o st a ć : _ tge> (7) D odatkowy  warun ek,  jaki  musi  speł niać rozwią zanie  r s ( - y r j  [(Y,)2  +  ( 7 2 ) 2 ] 3 ' 2 ,  (15) dcp  -   '  Y, gdzie:  T t   =   R.Qp),  Y 2   =   - ^ ^ - . Rozwią zując  problem  począ tkowy  sformuł owany  równaniami  (14) i  (15) oraz warunkami (10)  i  (13)  dla  dowolnie  dobranych  parametrów  A  i  S  otrzymamy  postać  powierzchni swobodnej,  która  z  okreś lonymi  bł ę dami  speł nia  warunki  (11)  i  (12).  Zakł adają c,  że istnieje  para  parametrów  A  i  S,  dla  której  warunki  (11)  i  (12)  są   speł nione dokł adnie, naszym  celem jest  znalezienie  rozwią zania  z  okreś lonym  bł ę dem  wzglę dnym  DT   dla  ką ta WYZN ACZAN IE  POWIERZCH N I  SWOBOD N EJ...  87 zwilż enia i  bł ę dem  DS  dla  pola  powierzchni  pod  krzywą gdzie  T  i P OL  są  obliczonymi  wartoś ciami  ką ta  zwilż enia  i pola  powierzchni  pod  krzywą dla  zadanych  parametrów  A  i  S.  Realizacja  postawionego  celu  polega  na  minimalizacji, przy  pomocy  omówionej  niż ej  procedury  numerycznej,  funkcji  bł ę du  w  postaci: FBL   =  DT 2 +DS Z .  (18) Jako  warunek  zakoń czenia  obliczeń  numerycznych  przyjmuje  się  odpowiednio  mał e wartoś ci  dla  bł ę dów  wzglę dnych  DT   i  DS. W  pracy  proponuje  się  minimalizowanie  wprowadzonej  funkcji  bł ę du  (18)  przy  po- mocy  metody  bę dą cej  poś rednią  pomię dzy  metodą  G aussa- Seidela  (bezgradientową) {[4],  str.  200}  a  metodą  najwię kszego  spadku  (gradientową)  {[4],  str.  203}.  Zasadniczo w tej  procedurze  powtarzają  się  na  przemian  dwa  etapy  obliczeń:  poszukiwanie  kierunku najwię kszego  spadku  oraz przeszukiwanie  tego  kierunku.  Z  uwagi  na  fakt,  że  w  badanym zagadnieniu  nie  dysponuje  się  informacją  o  gradiencie  funkcji  bł ę du  FBL   w  pierwszym etapie  poszukuje  się  minimum FBL  w  skoń czonej  iloś ci  punktów  na elipsie  w  pł aszczyź nie (A, S).  Ilość  tych  punktów,  osie  i  ś rodek  elipsy  okreś la  się  arbitralnie  na  począ tku  pro- cedury  poszukiwania  kierunku  najwię kszego  spadku.  Punkt  na  elipsie,  w  którym  FBL osią ga  minimum  i  ś rodek  elipsy  wyznaczają  kierunek  przybliż ony  do  kierunku  naj- wię kszego  spadku,  i jest  on  dalej  nazywany  kierunkiem  najwię kszego  spadku.  O  ile  FBL na  elipsie  nie  osią ga  wartoś ci  mniejszej  aniż eli  w  jej  ś rodku  nastę puje  zmniejszenie  osi elipsy.  W  kierunku  najwię kszego  spadku  minimum funkcji  bł ę du  poszukuje  się  przy  po- mocy  metody  przeszukiwania  ze  zmiennym  krokiem. D la przyspieszenia  obliczeń pierwsze  przybliż enie  rozwią zania  poszukiwano  dla  mał ej iloś ci  kroków  cał kowania  (oznaczonej  w  algorytmie  N P)  w  metodzie  Rungego- Kutty. N astę pnie  ilość  tych  kroków  zwię kszano  tak  dł ugo,  dopóki  wyniki  z  dwóch  kolejnych kroków  cał kowania  róż niły  się  mniej  niż  zał oż one  kryterium  (w  pracy  1%  uzyskanego wyniku).  N ależy  zwrócić  uwagę  na  fakt,  że  czas  obliczeń  jest  w  duż ym  stopniu  zależ ny od  trafnoś ci  przyję cia  pierwszej  pary  parametrów  A  i  £   oraz  kroków  procedury  przeszu- kiwania  (w programie  oznaczonych AA  dla A  i SW dla S).  D latego też w programie  wpro- wadzono  moż liwość  „ strzelania"  punktem  startowym  (danymi  A,  S,  AA,  SS,  N F)  tak, aby  uzyskane  pierwsze  przybliż enie  zapewniał o  moż liwie  krótki  czas  obliczeń. N a  rys.  2  podano  algorytm  obliczeń  programu  wykorzystanego  w  pracy,  który  jest zamieszczony  w dodatku. Oddzielnie na rys.  3 i 4 rozpisane został y fragmenty  tej procedury dotyczą ce  poszukiwania  kierunku  najwię kszego  spadku  oraz  poszukiwania  minimum na  tym  kierunk- u.  Algorytm  obliczeń  z  rys.  3,  2  i  4  ma  n a  celu  uł atwienie  ewentualnemu uż ytkownikowi  korzystania  z  zał ą czonego  programu  (obliczenia  wykonywano  na  H P S6B  w  ję zyku  Basic).. I C zyt a ć  :  8  ,  Hi .  Bp  ( | Czyt ać :  A  ,  S  ,  AA  ,  SS  ,  NP Ro zwią zan ie  z a g a d n i e n i a  p o c z ą t k o w e go jako  w y n i k :  FBL  ,  DT  ,  DS  ;  FB- FBL ± < | Wyś wietlenie  wyników  na  monitorze Czy  punkt  sta r t owy  obl i czeń  nas  zadawał  o ? \ n i e / p yt an ie  do  uży tkowni k a /   ^  ' \   tak Wydrukowanie  i  zapamię t anie  a k t u a l n e g o wyn i k u  j a k o  p unkt  bazowy / Czy  d l a  d aneg o  NP  uzyskano  ż a,danq  d o k ł a d n o ś ć ? > . X  | 0T|  < 1 0 ' 4  | DS| <   lO- 4  / We, We, Po szuk iwanie  k i e r u n k u  najwię kszego  s p a d k u f u n k c j i  b ł ę d u.  K i e r u  n e k  w y zn a c za j ą ,:  p u n k t b azowy  i  p u n k t  na  e l i p s i e  o  ś r o d ku  w  p u n - k c i e  b a zo w y m  i  o s i a c h  ZxSS  i  2 xAA 5 S= SS/ 5 AA=AA/ B CN^A , t ak Sp r awd zenie  zm iany  wyniku  sp owod owanej  zmianaX. lic zb y  kroków  c a ł k o w a n i a .  Czy  X ' c | (RN- S) / S |   > 10" 3  lub  | (CN- A)/ A|   > 1 0 " 3 1nie | Drukowanie  wyników- .  Bo  ,  NP ,  FB ,  DS ,  DT ,  S ,  A} ( "S T O P ) R ys.  2.  Algorytm  obliczeń  wyznaczania  kształ tu  powierzchni  swobodnej [88] r* L We; .We, Rozwią zanie  równania w dwóch  kierunkach  prostopa- dłych  do  ostatniego, najlep- szego  kierunku ;  wybór lepszego  z  nic h. Rozwią zanie  równania  w  czte- rech  kierunkach  przestrzeni (S,A):  ( S± SS, A) ,  ( S,A i AA) ; wybór  k ie r u n k u , dla  którego FBL=min;  BL^FBLmin /   C z y  i s t n i e j e  s ą s i e d ni  k i e r u n e k  ,  d a j ą cy w y n i k  ? t ak P r z y j ę c ie  t e g o  o s t a t n i e g o  za  n a j l e p s z y k i e r u n e K  ;  BL=  FB L s ą s i e d ni  k i e r u n e k  d a j e  l e p s z y  w y n i k ? m   n i e N t 7 Wy Rys.  3.  Algorytm  procedury  wyznaczania  kierunku  najwię kszego  spadku Pr zyję c ie  nowego  p unkt u  bazowego  ; nowe  S,  A ,  FB=BL,  DT ,  D5 PARIsPARl+1 W  zn a l e zi o n y m  u p r ze d n i o  k i e r u n k u  r o z- w i ą zu j e my  zag ad n i e n i e  z  k r o k i e m SSxPARl  .  AA»PARI;  o k r e ś l e n ie  FBL t a k FBL  < FB ? PARI  >   3 ? \   t ak Ro zwią zan ie  zag ad n i e n i a  od  p u n k t u  bazowe- l go  z  kr okiem  SS, AA  w  b ad anym  kier unku < Xzy  n a s t ą p i ła  p o p r awa  r o zw i ą za n ia  ? t a k 1 ,, nie Ro zwią zan ie  od  p u n k t u  b azo weg o  z  k r o - kiem  S S ,  AA  w  k i e r u n k u  p r ze c i w n y m n a s t ą p i ła  p o p r a w a  r o zwi  ą za n ia  ? t a k [ Zmiana  k ier unk u  poszukiwania  na  przeciwny  | Rys.  4.  Algorytm  procedury  metody  przeszukiwania [89] 90  M.  KACZMAREK  I  IN N I 4.  Przykł adowe  wyniki Przedstawiona  metoda  wyznaczania  powierzchni  swobodnej  może  być  wykorzystana dla  powierzchni  wypukł ych  zarówno  w  dodatnim  (A  >  S) jak  i  ujemnym  (A  <  S) kie- runku  osi  ukł adu  biegunowego.  Przykł adowe  przebiegi  powierzchni  swobodnych  po- kazano  na  rys.  5 i 6. N a rys.  5 przedstawiono  powierzchnie  swobodne  dla równych ką tów Rys.  5. Profile  powierzchni  swobodnej  dla  róż nych  Rys. 6. Profile powierzchni swobodnych dlaióż nych liczb  Bonda  ką tów  zwilż enia zwilż enia  i  róż nych  liczb  Bonda.  N a  rys.  6  przedstawiono  sytuację   odwrotną .  W  celu oszacowania  bł ę du  wynikają cego  z  zał oż enia  stał ego  promienia krzywizny  [1]  na  rys. 7 przedstawiono  zależ noś ci  promienia  krzywizny  od  współ rzę dnej biegunowej  ę   dla  róż- 3   2 - Rys.  7.  Bezwymiarowe  promienie  krzywizny  dla  róż nych  liczb  Bonda WYZN ACZAN IE  POWIERZCH N I SWOBODNEJ 91 nych  liczb  Bonda.  Z  przebiegów  moż na  stwierdzić,  że  popeł niony  bł ą d  jest  niewielki jeż eli liczba Bonda jest  mniejsza  od  1. Dla wię kszych liczb Bonda  bł ą d ten znacznie roś n ie. Uż ywany  przez  autorów  program  wymaga  okreś lenia  wyjś ciowej  pary  parametrów A  i  S.  Trafność  ich  doboru  decyduje  o  szybkoś ci  osią gnię cia  zadowalają cych  wyników a tę  pierwszą   uł atwia znajomość  charakteru zmian tych  wielkoś ci  w funkcji  liczby  Bonda. Przykł adowe  zależ noś ci  pokazują ce  jak  zmieniają   się   parametry  A  i  S  w  funkcji  liczby Bonda  przedstawiono  na rys.  8 i  9.  Znajomość jakoś ciowych  zależ noś ci  tych parametrów może  być  wykorzystana  w  trakcie  obliczeń  dotyczą cych  innych  ką tów  rozwarcia 

FBL  THEN  KATA=KATB  ®  BL«FBL    MOD  (2#F'l ) 580 C=A+AA*COS  FB THEN AA=AA*,2 @ SB=SS*.2 @ DISP "zmiana kroku AA,SS="; AA;3S @ 0OTQ 220 650 KAT- KATA @ BOTO 360 660 KATA- KAT1 e KAT1=(KATA+ZKAT) MOD  (2*PI ) 670 C=A+AA*CQ5  (KAT1) @ Y(1)=S+SS*SIN  (KAT1) 6S0 GOBUB RQZW  @ IF FBLFB THEN AA=AA*.2 i SS=SS*.2 © DISP "zmiana kroku AA,SS="; AA;S3 @ GOTO 220 7 0 0   KAT- KATA  S  GOTO  3 6 0 710  i  Sp r awd zen ie  d o k ł ad n o ś ci  uzyskanego  wyn ik u 720  SPR:  IF  ABS  ( DT) > . 0 0 0 1  DR  ABS  (DS)> .0001  THEN RETURN 730  IF  ABS  ( CRN- S) / S)  >.OO1  OR  ABS  ( . 001  THEN  DISP  " zmiana k r o k u  c a ł k o w a n i a  NP";NP* 2 740  IF  ABS  ( (RN- S)  / S)  >.OO1   OR  ABS  (   > . 0 0 1   THEN  PARN=1  > 920 D(I)=2*U*H»F  *QCI) 930 NEXT I 940 U=.2928932 950 GOTO F WYZN ACZAN IE  POWIERZCH N I  SWOBOD N EJ... 95 O O - C / Y U) )*SQR 960 E: FOR 1=1 TO 2 970  Y(I)=Y(I)+H*F(I)/6- Q(I)/3 980 NEXT 1 990 M=0 @ K=2 1000 GDTO (3 1010 F: K=l 1020 B; IF K=2 THEN K 1030 F(l)- Y(2) 1040  F<2)=Y(1)+2*Y(2)*- Y(2)/Y(1>+B0- *(COS Y(2) *Y(2) )'- '3 1050  GOTO LD 1060 K: AB=Y(1>*COS  (X) @ AC=YU)*SIN  (X) 1070  POLB=PQLB+(AB+G1>*(AC- G2) 1080 G1=AB © G2=AC @ KR=KR+1 1090 IF KR/TA