Bonjour, j'ai fais plusieurs codes sur les "CARACTÉRISATION DE SURFACE PAR LA MÉTHODE DU PONT CAPILLAIRE" en MATLAB,
j'aimerais savoir si quelqu'un pourrait traduire cela en SCILAB, svp !
La partie la plus dure est celle de la résolution d'équation différentielle de Laplace (du 2nd ordre). Savez vous résoudre celle ci sur SCILAB, sans passer par des solvers svp ?
Pour plus d'information vous pouvez m'envoyer un mail : janot92@hotmail.fr
Voici ci dessous les codes :
Fonction Dérivé de Y par rapport à z :
Fonction Point capillaire (en 2D) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5 function [ v ] = DYdz( z,Y ) %v : paramètre de sortie ; z et y : paramètres d'entrée global lc; %lc = kappa v=[Y(2) (1+Y(2)^2)/Y(1)+(1+Y(2)^2)^(3/2)*z/(lc^2)]; end
Fonction variation de Thêta (pour obtenir la courbe A0=f(), pour h=0) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 function [z,Y] = ptcap2D (h,theta_deg) global lc; warning off; lc=2; R=100; z_lentille=linspace(0,10,101); r_lentille=sqrt(2*R.*(z_lentille-h)-((z_lentille-h).^2)); zpas=1; zini=h+0.01; for iii=1:1:4 zpas=zpas/(10); for z0=zini:zpas:R/10 r0=sqrt(2*R.*(z0-h)-((z0-h)^2)); alpha=asin(r0/R); theta=theta_deg/180*pi; rp0=cot(alpha+theta); %solveur [z,Y]=ode45('DYdz',[z0 0],[r0 rp0]); taille=size(z); zfinal=z(taille(1)); taille2=size(Y); y2final=Y(2,taille2(2)); if (zfinal~=0 || abs(y2final)>100) zini=z0-zpas; [z,Y]=ode45('DYdz',[zini 0],[r0 rp0]); break end end end %tracés (Attention : il faut mettre ces tracés en commentaire pour obtenir les autres courbes) hold on; plot(r_lentille,z_lentille); plot(-r_lentille,z_lentille); plot(Y(:,1),z); plot(-Y(:,1),z); axis equal; axis([-60 60 0 10]) end Fonction Pont capillaire (en 3D) : function [X,Y,Zlentille] = ptcap3D ( ) R=100; h=0; for zi=0:1:100 r_lentille=sqrt(2.*R.*(zi/10-h)-(zi/10-h).^2); for the=0:100 X(zi+1,the+1)=r_lentille.*cos(the*3.6*pi/180); Y(zi+1,the+1)=r_lentille.*sin(the*3.6*pi/180); end end for xi=1:101 for yi=1:101 Zlentille(xi,yi)=((2*R+2*h)... -sqrt(((-2*R-2*h)^2)-4.*(2*R*h+h^2+X(xi).^2+Y(yi).^2)))./2; end end [Zpt,Ytmp]=ptcap2D(h,0); % on se sert de ptcap2D tailleY=size(Ytmp); for zi=1:tailleY(1,1) r_pt=Ytmp(zi,1); for the=1:tailleY(1,1) Xpt(zi,the)=r_pt.*cos((the-1)/tailleY(1,1)*360*pi/180); Ypt(zi,the)=r_pt.*sin((the-1)/tailleY(1,1)*360*pi/180); end end for xi=1:tailleY(1,1) for yi=1:tailleY(1,1) Zpt_tab(xi,yi)=Zpt(xi); end end %tracés du graphe 3D figure( ) hold on mesh(X,Y,Zlent); % la function mess équivaut à un plot, mais en 3D mesh(Xpt,Ypt,Zpt_tab); colormap cool; %colormap est une fonction permettant d'avoir plusieurs coloris du graphe axis equal; hold off; end
Fonction variation de h (pour obtenir la courbe A=f(h) pour différents fixé) :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14 function variation_theta () h = 0; kkk = 0; for theta_deg=0:1:178 kkk = kkk +1; [z,Y] =ptcap2D(h,theta_deg); tab_theta(kkk,1)=theta_deg; % on stock les valeurs de theta_deg dans un tableau Aire(kkk,1)=pi*Y(1,1)^2; %pi*Y(1,1)^2 end hold on plot(tab_theta,Aire,,'+r'); hold off end
Code : Sélectionner tout - Visualiser dans une fenêtre à part
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 function variation_h() num_theta=0 ; for theta_deg=0:20:180 %faire varier h avec un thêta fixe qui se trouve dans argument grâce à un for % ici on varie thêta entre 0° et 180° num_theta= num_theta+1 ; hpas=0.5; %Aire=zeros(100,1); %tab_h=zeros(100,1); jjj=0; for jjj=0:1:1 kkk=0; for h=0:hpas:100 kkk=kkk+1; [z,Y] =ptcap2D(h,theta_deg); %h tab_h(kkk, num_theta)=h; Aire(kkk, num_theta)=pi*Y(1,1)^2; %pi*Y(1,1)^2 taille=size(z); zfinal=z(taille(1)); if(zfinal~=0) if(kkk<=30) hpas=hpas/10; break; else break; end end %utiliser ptcap2D pour avoir dernière valeur du tableau z qui correspond au pont capillaire %si cette valeur est différente de 0 on arrête avec un break %On doit prendre la première valeur de Y pour calculer l'aire de contact %pi*r2=pi*Y2 => aire de contact %tracer a(h) if(kkk>30) break; end end end hold on plot(tab_h,Aire,'+r'); hold off end end
Partager