Bonjour;
Merci de m'expliquer pourquoi j'obtiens des valeurs en nombres complexes pour F et f.
Code:clc;clear all;close all; %%%%Geometry of domain%%%%%%%%%%% W=1.95; % width of domain L=1.22; % High of domain Hs=.85; %Nr=55;%51; %No of grids in x diretcion Nr=55; Nz=55; %101;%No of grids in z diretcion dr=W/(Nr-1); % width of element dz=L/(Nz-1); % height of element nt = 10; dt=1; ks=0.0045949074; nV=4; mV=1-1/nV; alphaV=9; thetas=.43; thetar=0; %x=(0:Nr)*dr; y=(0:Nz)*dy; % x , y values on the grid %%%% Boundary conditions %%% F=ones(Nr,Nz);%zeros(Nr,Nz);% Ss=ones(Nr,Nz); f=zeros(Nr,Nz); wt=ones(Nr,Nz); C=zeros(Nr,Nz); k=zeros(Nr,Nz); T=20; Z=[0:dz:L]; F(Nr,:)=1.22; Ndt=round(T/dt); % Nombre de pas de temps t=0; for nz=1:Nz f(Nr,:)=.305 - nz*dz; end %k(Nr,:)=1000*ks; k(Nr,:)=ks; F(round(.122/dr):round(1.95/dr),1)=F(round(.122/dr):round(1.95/dr),2); %Bottom condition for F f(round(.122/dr):round(1.95/dr),1)=F(round(.122/dr):round(1.95/dr),1); %Bottom condition for f k(round(.122/dr):round(1.95/dr),1)=0; %Bottom condition for k F(round(.122/dr):round(1.95/dr),Nz)=F(round(.122/dr):round(1.95/dr),Nz-1); %Top condition for F %f(round(.122/dr):round(1.95/dr),Nz)=0; %k(round(.122/dr):round(1.95/dr),Nz)=1000*ks; %Top condition for k F(round(.122/dr),1:round(0.305/dz))=0.305; F(1,1:round(0.305/dz))=0.305; for nz=1:round(1.22/dz) f(round(1.95/dr),nz)=0.305; end for nz=1:round(0.305/dz) f(round(.122/dr),nz)=0.305-nz*dz; end k(round(.122/dr),1:round(.305/dz))=1000*ks; k(1,1:round(.305/dz))=1000*ks; for nz=round(.305/dz):round(Hs/dz) F(round(.122/dr),nz)=nz*dz; f(round(.122/dr),nz)=0; k(round(.122/dr),nz)=1000*ks; F(1,nz)=nz*dz; %f(1,nz)=0; k(1,nz)=1000*ks; end for nz=round(Hs/dz):Nz F(round(.122/dr),nz)=F(round(.122/dr)-1,nz); %f(round(.122/dr),nz)=F(round(.122/dr)+1,nz); k(round(.122/dr),nz)=1000*ks; end %%%%% Computation %%%%%%%% f_old=f; %F_intial=F; f_intial=f; %k_initial=k; Gauss_Seidel_iteration = 1; for nx=1:Ndt % Incrément de temps t=t+dt; tolerance=1e-5; Iter=0; for m = 1:nt %error = 9e-9; error=.5; while(error > tolerance) %Iter=Iter+1; %disp(Iter); %Hold=F; for nz=2:Nz-1 for nx=round(.122/dr)+1:Nr-1 a=(k(nx,nz)+k(nx,nz-1))/(2*dz^2); b=((dr*nx).*k(nx,nz)+(dr*(nx-1)).*k(nx-1,nz))./(2*(dr*nx)*dr^2); d=((dr*nx).*k(nx,nz)+(dr*(nx+1)).*k(nx+1,nz))./(2*(dr*nx)*dr^2); e=(k(nx,nz)+k(nx,nz+1))/(2*dz^2); f1=C(nx,nz)/dt; g=-(k(nx,nz+1)+k(nx,nz-1))/(2*dz); c=-(a+b+d+e+f1); f(nx,nz)=-[a.*f(nx,nz-1)+b.*f(nx-1,nz)+d.*f(nx+1,nz)+e.*f(nx,nz+1)+f1.*f_old(nx,nz)+g]./c; Ss(nx,nz)=(1+(alphaV*abs(min(f(nx,nz),0))).^nV).^(-mV); %Ss(nx,nz)=(1+(-alphaV*f(nx,nz))^nV)^(-mV); C(nx,nz)=alphaV*(thetas-thetar)*(nV-1)*(alphaV*abs(min(f(nx,nz),0))).^(nV-1)./... ((1+(alphaV*abs(min(f(nx,nz),0))).^nV).^((2*nV-1)/nV)); % C(nx,nz)=alphaV*(thetas-thetar)*(nV-1)*(Ss(nx,nz).^(nV/(nV-1))-1).^((nV-1)/nV)./... %((Ss(nx,nz)).^((1-2*nV)/(nV-1))); k(nx,nz)=ks*(1-(1-Ss(nx,nz).^(-mV)).^mV).^2.*Ss(nx,nz).^.5; F(nx,nz)=f(nx,nz)+nz*dz; % F(nx,nz)=[k(nx,nz).*((F(nx+1,nz)+F(nx-1,nz-1))/dr^2+(F(nx,nz+1)+F(nx,nz-1))/dz^2)+(k(nx+1,nz)-k(nx-1,nz)).*F(nx+1,nz)/(2*dr^2)+... %(k(nx,nz+1)-k(nx,nz-1)).*F(nx,nz+1)/(2*dz^2)]./[2*k(nx,nz)*(1/dr^2+1/dz^2)(k(nx+1,nz)-k(nx-1,nz))/(2*dr^2)+... %(k(nx,nz+1)-k(nx,nz-1))/(2*dz^2)]; if f(nx,nz)==0; wt(nx,nz)=nz*dz; %F(nx,nz); endif endfor endfor %Error=sqrt(sumsq(F-Hold)); %Error=sqrt(F-Fold); %Error=max(abs(F-Fold)); error = max(abs(f_old - f)); %max(max(abs(f_old - f))); f_old = f; Gauss_Seidel_iteration = Gauss_Seidel_iteration + 1; end f_intial = f; end end %%% Plotting the results %%% x=0:dr:W; % width array (in m) z=0:dz:L; % Height array (in m) %[X,Y] =meshgrid(x,z); %v=[0.8 0.6 0.4 0.2 0.1 0.05 0.01]; %colormap(jet); %contourf(x, z, f'); %surf(f) %shading flat %plot(x,F) %plot(x,wt) plot(x,f) title("Charge distribution"); xlabel("Radial coordinate (in m)"); ylabel("Position (in m)"); title(sprintf('No. of Unsteady Gauss-Seidel Iterations(IMPLICIT) = %d', Gauss_Seidel_iteration));
-----