{\rtf1\ansi\ansicpg936\cocoartf1404\cocoasubrtf470 {\fonttbl\f0\fmodern\fcharset0 Courier;} {\colortbl;\red255\green255\blue255;\red0\green0\blue0;} \paperw11900\paperh16840\margl1440\margr1440\vieww10800\viewh8400\viewkind0 \deftab720 \pard\pardeftab720\sl280\partightenfactor0 \f0\fs24 \cf2 \expnd0\expndtw0\kerning0 \outl0\strokewidth0 \strokec2 function T=tubes(h,z,r)\ \ l=0.3;\ k1=0.58;\ k2=0.023;\ Q=2e4;\ %z=0;\ a=1;\ H=0.6;\ r0=2.5e-3;\ %r0=2.5e-3;\ mu=8e-4;\ alpha=800;\ rho=1000;\ g=9.8;\ bess=[3.83171,7.01559,10.1735,13.3237,16.4706]/r0;\ \ G=rho*g*(h-H/3)/l;\ A=-k2*Q/(k2*(alpha*h+2*k1)+k1*alpha*(H-h));\ B=(k1*A+Q)/alpha;\ T0=A*H/3+B;\ Th=A*h+B;\ a(1)=(T0+(G^2)*(r0^4)/(192*mu*k1))/l;\ T=a(1)*(l-z);\ for i=1:length(bess)\ a(i+1)=r0*besselj(3,r0*bess(i))*(8-r0^2*(bess(i)^2))*(G^2)/(32*besselj(0,r0*bess(i))^2*bess(i)^3*mu*k1*sinh(bess(i)*l));\ if abs(a(i+1))>1e-4*abs(a(i))\ T=T+a(i+1)*sinh(bess(i)*(l-z))*besselj(0,r*bess(i));\ end\ end\ T=T-(G^2)*(r^4)/(64*mu*k1);\ %Th\ %T0\ \ }