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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
| %%%%%%%
%Mold dimensions and mesh elements
clear all, close all; clc,
Lx = 0.600; Nx = 60; hx = Lx / (Nx);
Ly = 0.300; Ny = 30; hy = Ly / (Ny);
%Loading LIMS/LB output files:
%%A text file containing 8 columns: number of the node, node coordinates
%%(x,y and z respectively),pressure, flow rate, time and fill factor in all
%%the nodes
B1 = load('case_0_8.txt'); Pb1= B1(:,5);
B2 = load('case_1_8.txt'); Pb2= B2(:,5);
B3 = load('case_2_8.txt'); Pb3= B3(:,5);
B4 = load('case_3_8.txt'); Pb4= B4(:,5);
B5 = load('case_4_8.txt'); Pb5= B5(:,5);
B6 = load('case_5_8.txt'); Pb6= B6(:,5);
B7 = load('case_6_8.txt'); Pb7= B7(:,5);
B8 = load('case_7_8.txt'); Pb8= B8(:,5);
B9 = load('case_8_8.txt'); Pb9= B9(:,5);
B10 = load('case_9_8.txt'); Pb10= B10(:,5);
B11 = load('case_10_8.txt');Pb11= B11(:,5);
B12 = load('case_11_8.txt');Pb12= B12(:,5);
B13 = load('case_12_8.txt');Pb13= B13(:,5);
B14 = load('case_13_8.txt');Pb14= B14(:,5);
B15 = load('case_14_8.txt');Pb15= B15(:,5);
for j = 1:Ny+1
for i = 1:Nx+1
k = (j-1)*(Nx+1) + i;
PRESS1(i,j)= Pb1(k);
PRESS2(i,j)= Pb2(k);
PRESS3(i,j)= Pb3(k);
PRESS4(i,j)= Pb4(k);
PRESS5(i,j)= Pb5(k);
PRESS6(i,j)= Pb6(k);
PRESS7(i,j)= Pb7(k);
PRESS8(i,j)= Pb8(k);
PRESS9(i,j)= Pb9(k);
PRESS10(i,j)= Pb10(k);
PRESS11(i,j)= Pb11(k);
PRESS12(i,j)= Pb12(k);
PRESS13(i,j)= Pb13(k);
PRESS14(i,j)= Pb14(k);
PRESS15(i,j)= Pb15(k);
end
end
PRESS1 = PRESS1'; PRESS1=PRESS1(end:-1:1,:);
PRESS2 = PRESS2'; PRESS2=PRESS2(end:-1:1,:);
PRESS3 = PRESS3'; PRESS3=PRESS3(end:-1:1,:);
PRESS4 = PRESS4'; PRESS4=PRESS4(end:-1:1,:);
PRESS5 = PRESS5'; PRESS5=PRESS5(end:-1:1,:);
PRESS6 = PRESS6'; PRESS6=PRESS6(end:-1:1,:);
PRESS7 = PRESS7'; PRESS7=PRESS7(end:-1:1,:);
PRESS8 = PRESS8'; PRESS8=PRESS8(end:-1:1,:);
PRESS9 = PRESS9'; PRESS9=PRESS9(end:-1:1,:);
PRESS10 = PRESS10'; PRESS10=PRESS10(end:-1:1,:);
PRESS11 = PRESS11'; PRESS11=PRESS11(end:-1:1,:);
PRESS12 = PRESS12'; PRESS12=PRESS12(end:-1:1,:);
PRESS13 = PRESS13'; PRESS13=PRESS13(end:-1:1,:);
PRESS14 = PRESS14'; PRESS14=PRESS14(end:-1:1,:);
PRESS15 = PRESS15'; PRESS15=PRESS15(end:-1:1,:);
PRESST=PRESS1+PRESS2+PRESS3+PRESS4+PRESS5+PRESS6+PRESS7+PRESS8+PRESS9+PRESS10+PRESS11+PRESS12+PRESS13+PRESS14+PRESS15;
[X,Y] = meshgrid(0:hx:Lx, Ly:-hy:0);
% Using the built-in function gradient directly on the pressure values (map1)
figure(1)
colormap('jet');
[px,py] = gradient(PRESST,hx,hy);
pxpy =sqrt(px.^2+py.^2);
contour(pxpy), quiver(px,py)
contour(X,Y,pxpy), shading flat
[aa,bb] = contourf(X,Y,pxpy); clabel(aa,bb), colorbar
title('map 1')
% Calculating pressure gradients in x and y directions and taking the norm
PRESS1_dx = [zeros(Ny+1,1) (diff(PRESS1,1,2))]/hx; PRESS1_dy = [zeros(1,Nx+1); (diff(PRESS1,1,1))]/hy; grad1=sqrt(PRESS1_dx.^2+PRESS1_dy.^2);
PRESS2_dx = [zeros(Ny+1,1) (diff(PRESS2,1,2))]/hx; PRESS2_dy = [zeros(1,Nx+1); (diff(PRESS2,1,1))]/hy; grad2=sqrt(PRESS2_dx.^2+PRESS2_dy.^2);
PRESS3_dx = [zeros(Ny+1,1) (diff(PRESS3,1,2))]/hx; PRESS3_dy = [zeros(1,Nx+1); (diff(PRESS3,1,1))]/hy; grad3=sqrt(PRESS3_dx.^2+PRESS3_dy.^2);
PRESS4_dx = [zeros(Ny+1,1) (diff(PRESS4,1,2))]/hx; PRESS4_dy = [zeros(1,Nx+1); (diff(PRESS4,1,1))]/hy; grad4=sqrt(PRESS4_dx.^2+PRESS4_dy.^2);
PRESS5_dx = [zeros(Ny+1,1) (diff(PRESS5,1,2))]/hx; PRESS5_dy = [zeros(1,Nx+1); (diff(PRESS5,1,1))]/hy; grad5=sqrt(PRESS5_dx.^2+PRESS5_dy.^2);
PRESS6_dx = [zeros(Ny+1,1) (diff(PRESS6,1,2))]/hx; PRESS6_dy = [zeros(1,Nx+1); (diff(PRESS6,1,1))]/hy; grad6=sqrt(PRESS6_dx.^2+PRESS6_dy.^2);
PRESS7_dx = [zeros(Ny+1,1) (diff(PRESS7,1,2))]/hx; PRESS7_dy = [zeros(1,Nx+1); (diff(PRESS7,1,1))]/hy; grad7=sqrt(PRESS7_dx.^2+PRESS7_dy.^2);
PRESS8_dx = [zeros(Ny+1,1) (diff(PRESS8,1,2))]/hx; PRESS8_dy = [zeros(1,Nx+1); (diff(PRESS8,1,1))]/hy; grad8=sqrt(PRESS8_dx.^2+PRESS8_dy.^2);
PRESS9_dx = [zeros(Ny+1,1) (diff(PRESS9,1,2))]/hx; PRESS9_dy = [zeros(1,Nx+1); (diff(PRESS9,1,1))]/hy; grad9=sqrt(PRESS9_dx.^2+PRESS9_dy.^2);
PRESS10_dx = [zeros(Ny+1,1) (diff(PRESS10,1,2))]/hx; PRESS10_dy = [zeros(1,Nx+1); (diff(PRESS10,1,1))]/hy; grad10=sqrt(PRESS10_dx.^2+PRESS10_dy.^2);
PRESS11_dx = [zeros(Ny+1,1) (diff(PRESS11,1,2))]/hx; PRESS11_dy = [zeros(1,Nx+1); (diff(PRESS11,1,1))]/hy; grad11=sqrt(PRESS11_dx.^2+PRESS11_dy.^2);
PRESS12_dx = [zeros(Ny+1,1) (diff(PRESS12,1,2))]/hx; PRESS12_dy = [zeros(1,Nx+1); (diff(PRESS12,1,1))]/hy; grad12=sqrt(PRESS12_dx.^2+PRESS12_dy.^2);
PRESS13_dx = [zeros(Ny+1,1) (diff(PRESS13,1,2))]/hx; PRESS13_dy = [zeros(1,Nx+1); (diff(PRESS13,1,1))]/hy; grad13=sqrt(PRESS13_dx.^2+PRESS13_dy.^2);
PRESS14_dx = [zeros(Ny+1,1) (diff(PRESS14,1,2))]/hx; PRESS14_dy = [zeros(1,Nx+1); (diff(PRESS14,1,1))]/hy; grad14=sqrt(PRESS14_dx.^2+PRESS14_dy.^2);
PRESS15_dx = [zeros(Ny+1,1) (diff(PRESS15,1,2))]/hx; PRESS15_dy = [zeros(1,Nx+1); (diff(PRESS15,1,1))]/hy; grad15=sqrt(PRESS15_dx.^2+PRESS15_dy.^2);
% Superimposed pressure gradients
GRADT=grad1+grad2+grad3+grad4+grad5+grad6+grad7+grad8+grad9+grad10+grad11+grad12+grad13+grad14+grad15;
% Finding the 3 highest pressure gradients zones
Obj=reshape(GRADT,1,numel(GRADT));
Xf=reshape(X,1,numel(X));
Yf=reshape(Y,1,numel(Y));
[B,I] = sort(Obj);
X_op = Xf(I(numel(I):-1:1));
Y_op = Yf(I(numel(I):-1:1));
i=1;
j=2;
X1=X_op(1)
Y1=Y_op(1)
p=0.2;
while j<= length(Obj)
if i<2
if sqrt((X_op(j)-X1)^2+(Y_op(j)-Y1)^2)<= p
j=j+1;
else
X2=X_op(j)
Y2=Y_op(j)
i=i+1;
end
else
if sqrt((X_op(j)-X1)^2+(Y_op(j)-Y1)^2)<= 0.2 || sqrt((X_op(j)-X2)^2+(Y_op(j)-Y2)^2)<= p
j=j+1;
else
X3=X_op(j)
Y3=Y_op(j)
i=i+1
break
end
end
end
% plotting pressure gradients map2 and sensors locations
figure(2);
colormap('jet');
contour(X,Y,GRADT), hold on, shading flat
[aa,bb] = contourf(X,Y,GRADT); colorbar, clabel(aa,bb)
title('map 2')
hold on
plot(X1,Y1,'k*',X2,Y2,'k*',X3,Y3,'k*','linewidth', 10)
txt1 = ' 1';txt2 = ' 2';txt3 = ' 3';
text(X1,Y1,txt1,'FontSize', 20), text(X2,Y2,txt2,'FontSize', 20), text(X3,Y3,txt3,'FontSize', 20) |
Calculs des gradients et trouver les coordonnées du N le plus élevé
Répondre avec citation
Partager