Advanced Mathematics and Mechanics Applications Using MATLAB phần 9 pps

86: print(’Input data are incorrect. The ’);

87: print(’following r values lie outside the ’);

88: print(’unit circle:’); disp(rvec(kout)’); 89: return

90: end 91:

92: if bvtyp==1 % Solve a Dirichlet problem

93: % Check for points on the boundary where 94: % function values are known. Interpolate

95: % these directly

96: konbd=find(r==1); onbndry=length(konbd); 97: if onbndry > 0

98: u(konbd)=lintrp(thbv,fbv,th(konbd));

99: end 100:

101: % Evaluate the series solution 102: kinsid=find(r<1); inside=length(kinsid);

103:

104: if inside > 0 105: a=fft(lintrp(thbv,fbv,thfft));

106: a=a(1:nsum)/(nft/2);

107: a(1)=a(1)/2; Z=r(kinsid).*exp(i*th(kinsid)); 108: u(kinsid)=real(polyval(flipud(a(:)),Z));

109: end 110:

111: titl= ...

112: ’Dirichlet Problem Inside the Unit Circle’; 113:

114: else % Solve a Neumann problem

115: gbv=lintrp(thbv,fbv,thfft); 116: a=fft(gbv)/(nft/2);

117: erchek=abs(a(1))/sum(abs(gbv)); 118: if erchek>1e-3

119: disp(’ ’);

120: disp(’ERROR DUE TO NONZERO AVERAGE VALUE’); 121: disp(’OF NORMAL GRADIENT ON THE BOUNDARY.’);

122: disp(’CORRECT THE INPUT DATA AND RERUN.’);

123: disp(’ ’); u=[]; r=[]; th=[]; return; 124: end

125: a=a(2:nsum)./(1:nsum-1)’; z=r.*exp(i*th); 126: u=real(polyval(flipud([0;a(:)]),z));

127: titl=’Neumann Problem Inside the Unit Circle’;

128: end 129:

130: u=reshape(u,nth,nr); r=R; th=Th;

© 2003 by CRC Press LLC

131: surf(r.*cos(th),r.*sin(th),u);

132: xlabel(’x axis’); ylabel(’y axis’);

133: zlabel(’function u’); title(titl); 134: colormap(’default’);

135: grid on; figure(gcf); 136: % print -deps dirich

137:

138: %============================================= 139:

140: % function y=lintrp(xd,yd,x)

141: % See Appendix B

Function cauchtst

1: function u=cauchtst(z,nquad)

2: %

3: % u=cauchtst(z,nquad)

4: % ~~~~~~~~~~~~~~~~~~~

5: %

6: % This function solves a mixed boundary

7: % value problem for the interior of a circle

8: % by numerically evaluating a Cauchy integral.

9: % 10: % z - matrix of complex coordinates where

11: % function values are computed

12: % nquad - order of Gauss quadrature used to 13: % perform numerical integration

14: % 15: % u - computed values of the approximate

16: % solution

17: % 18: % User m functions called: cauchint, mbvtest,

19: % gcquad, splined

20:

21: if nargin<2, nquad=50; end; nbdat=61;

22: if nargin==0 23: z=linspace(0,.99,10)’* ...

24: exp(i*linspace(0,2*pi,91));

25: end 26: th=linspace(-pi/2,pi/2,nbdat); zb=exp(i*th);

27: fb=sqrt(zb-i).*sqrt(zb+i); fb(1)=1; fb(nbdat)=1;

© 2003 by CRC Press LLC

28: fb=cos(th)./fb; fb(1)=0; fb(end)=0;

29: F=cauchint(fb,zb,z,nquad);

30: F=F.*sqrt(z-i).*sqrt(z+i); u=2*real(F); 31:

32: surf(real(z),imag(z),u); xlabel(’x axis’); 33: ylabel(’y axis’); zlabel(’Solution Value’)

34: title([’Approximate Solution to ’, ...

35: ’a Mixed Boundary Value Problem’]); 36: grid on; figure(gcf); %gra(.4);

37: fprintf(’\nPress [Enter] to solution error\n’);

38: pause 39: %print -deps caucher1

40: uexact=mbvtest(z,1); udif=u-uexact;

41: clf; surf(real(z),imag(z),udif); 42: title([’Difference Between Exact and ’, ...

43: ’Approximate Solutions’]); 44: xlabel(’x axis’); ylabel(’y axis’);

45: zlabel(’Solution Error’)

46: grid on; figure(gcf); %gra(.4) 47: %print -deps caucher2

48:

49: %============================================= 50:

51: function u=mbvtest(z,noplot) 52: %

53: % u=mbvtest(z,noplot)

54: % ~~~~~~~~~~~~~~~~~~~ 55: %

56: % This function determines a function which is

57: % harmonic for abs(z)<1 and satisfies at r=1, 58: % u=cos(theta), -pi/2<theta<pi/2

59: % du/dr=0, pi/2<theta<3*pi/2 60: % The solution only applies for points inside

61: % or on the unit circle.

62: % 63: % z - matrix of complex values where the

64: % solution is computed.

65: % noplot - option set to one if no plot is 66: % requested, otherwise option is not

67: % required. 68: %

69: % u - values of the harmonic function

70: % defined inside the unit circle 71: %

72: % User m functions called: none

© 2003 by CRC Press LLC

73: %----------------------------------------------

74:

75: if nargin==0 76: noplot=0;

77: z=linspace(0,1,10)’* ... 78: exp(i*linspace(0,2*pi,81));

79: end

80: [n,m]=size(z); z=z(:); u=1/2*ones(size(z)); 81: k=find(abs(z)>0); Z=z(k);

82: U=(Z+1./Z+(1-1./Z).*sqrt(Z-i).*sqrt(Z+i))/2;

83: u(k)=real(U); u=reshape(u,n,m); 84: if nargin==1 | noplot==0

85: z=reshape(z,n,m);

86: surf(real(z),imag(z),u); xlabel(’x axis’); 87: ylabel(’y axis’);

88: title([’Mixed Boundary Value Problem ’, ... 89: ’for a Circular Disk’]);

90: grid; figure(gcf); %gra(.4), pause

91: %print -deps mbvtest 92: end

93:

94: %============================================= 95:

96: function F=cauchint(fb,zb,z,nquad) 97: %

98: % F=cauchint(fb,zb,z,nquad)

99: % ~~~~~~~~~~~~~~~~~~~~~~~~~ 100: %

101: % This function numerically evaluates a Cauchy

102: % integral of the form: 103: %

104: % F(z)=1/(2*pi*i)*Integral(f(t)/(t-z)*dt) 105: %

106: % where t denotes points on a curve in the

107: % complex plane. The boundary curve is defined 108: % by spline interpolation through data points

109: % zb lying on the curve. The values of f(t)

110: % are also specified by spline interpolation 111: % through values fb corresponding to the

112: % points zb. Numerical evaluation of the 113: % integral is performed using a composite

114: % Gauss formula of arbitrary order.

115: % 116: % fb - values of density function f

117: % at point on the curve

© 2003 by CRC Press LLC

118: % zb - points where fb is given. The

119: % number of values of zb must be

120: % adequate to define the curve 121: % accurately.

122: % z - a matrix of values at which the 123: % Cauchy integral is to be evaluated.

124: % If any of the z-values lie on path

125: % of integration or too close to the 126: % path of integration, incorrect

127: % results will be obtained.

128: % nquad - the order of Gauss quadrature 129: % formula used to perform numerical

130: % integration

131: % 132: % F - The value of the Cauchy integral

133: % corresponding to matrix argument z 134: %

135: % User m functions called: gcquad splined

136: %---------------------------------------------- 137:

138: n=length(fb); [nr,nc]=size(z); z=z(:).’;

139: nz=length(z); t=1:n; 140: [dummy,bp,wf]=gcquad(’’,1,n,nquad,n-1);

141: fq=spline(t,fb,bp); zq=spline(t,zb,bp); 142: zqd=splined(t,zb,bp); nq=length(fq);

143: fq=fq(:).*zqd(:);

144:

145: bdrylen=sum(abs(zq(2:nq)-zq(1:nq-1)));

146:

147: closnes=1e100; bigz=max(abs(z)); 148: for j=1:nq

149: closnes=min([closnes,abs(zq(j)-z)]); 150: end

151: if closnes/bdrylen<.01 | closnes/bigz<.01

152: disp(’ ’) 153: disp([’WARNING! SOME DATA VALUES ARE ’, ...

154: ’EITHER NEAR OR ON’]);

155: disp([’THE BOUNDARY. COMPUTED RESULTS ’, ... 156: ’MAY BE INACCURATE’]); disp(’ ’)

157: end 158: F=wf(:)’*(fq(:,ones(1,nz))./(zq(:,ones(1,nz))...

159: -z(ones(nq,1),:)));

160: F=reshape(F,nr,nc)/(2*pi*i); 161:

162: %=============================================

© 2003 by CRC Press LLC

163:

164: % function [val,bp,wf]=gcquad(func,xlow,...

165: % xhigh,nquad,mparts,varargin) 166: % See Appendix B

167:

168: %=============================================

169:

170: % function val=splined(xd,yd,x,if2) 171: % See Appendix B

12.12 Inviscid Fluid Flow around an Elliptic Cylinder

This section analyzes inviscid flow around an elliptic cylinder in an infinite field.

Flow around a circular cylinder is treated first. Then the function conformally map￾ping the exterior of a circle onto the exterior of an ellipse is used in conjunction with

the invariance of harmonic functions under a conformal transformation. Results de￾scribing the elliptic cylinder flow field for uniform velocity components at infinity

are presented.

Let us solve for the flow around a circular cylinder in the region |ζ| ≥ 1, ζ = ξ+iη

with the requirement that the velocity components at infinity have constant values

u = U,v = V

where (u, v) are the horizontal and vertical components of velocity. These compo￾nents are derivable from a potential function φ such that

u = ∂φ

∂ξ , v = ∂φ

∂η

where φ is a harmonic function. The velocity normal to the cylinder boundary must

be zero. This requires that the function ψ, the harmonic conjugate of φ, must be con￾stant on the boundary. The constant can be taken as zero without loss of generality.

In terms of the complex velocity potential

f(ζ) = φ + iψ

we need

f(ζ) − f(ζ)=0 on |ζ| = 1

The velocity field is related to the complex velocity potential by

u − iv = f

(ζ)

so the flow condition at infinity is satisfied by

f(ζ) = pζ + O(1) where p = U − iV

© 2003 by CRC Press LLC