Różnice

Różnice między wybraną wersją a wersją aktualną.

--- studia:magisterskie:1sem:modelowanie_procesow_fizycznych [2016/04/20 08:41]
149.156.112.6
+++ studia:magisterskie:1sem:modelowanie_procesow_fizycznych [2016/06/01 08:49] (aktualna)
91.227.212.3 [Drugi Program Maciek]
@@ Linia 42: / Linia 42: @@
 </code>
 ==== Zajęcia 3 ====
-<code matlab>temp=727;
+<code matlab>temp = 727;
+R = 8.3144;
 Q = 140000;
-R = 8.3144;
 temp_k = temp + 273;
-d0 = 0.000041;
+d0= 0.000041;
-D=d0*exp(-Q /(R * temp_k))*1E10;
+D = d0 * exp( -Q /( R * temp_k)) * 1 * 10^10;
 dt = 0.01;
 dx = 0.1;
 pN = 0.1;
-%Kryterium stabilności Neumanna:
-if (D*dt)/dx^2<=0.5
+if (D*dt)/dx^2 <= 0.5
-%(D*dt)/dx^2<=0.5
-%Obszar rozwiązania: 100 mm
+    % Obszar rozwiązania: 100 um
-%pierwszy i ostatni element powielamy
+    %y=[0.67 0.67 0.67 0.67 0.67 0.02 0.02 0.02 0.02];
-    K(6:100) = 0.02;
+    y(1:5) = 0.67;
-    K(1:5) = 0.67;
+    y(6:100) = 0.02;
     ksi = 6;
     for k = 0:dt:10
-        K = [K(1), K, K(length(K))];
+        y = [y(1), y, y(length(y))];
-        K1 = K;
+        y1 = y;
-        for i = 2:1:ksi+1%length(K)-1
+        for i = 2:1:ksi+1
-           K1(i) = D*(dt/dx^2)*(K(i+1)+K(i-1)) + K(i)*(1 - 2*D*(dt/dx^2));
+           y1(i) = D *( dt / dx^2 ) * ( y(i+1) + y(i-1)) + y(i) * (1 - 2*D*(dt/dx^2));
         end
-        K = K1;
+        y = y1;
-        K = K(2:1:length(K)-1);
+        y = y(2:1:length(y)-1);
-        plot(K)
+        plot(y)
         hold on
-        y=-0.004074074074074075*temp_k+3.7155555555555555;
+        % rownanie prostej przechodzacej przez dwa punkty
-        if K(ksi) >= y && ksi < length(K)-1
+        % (x2 - x1)(y - y1) = (y2 - y1)(x - x1)
-            ksi = ksi + 1
+        yLine = -0.00407 * temp_k + 3.715;
+        % zwiekszenie ksi
+        if y(ksi) >= yLine && ksi < length(y)-1
+            ksi = ksi + 1;
         end
+        % zwiekszenie temperatury o zadany prog i kolejna iteracja
         temp_k = temp_k + pN * dt;
-        D=d0*exp(-Q /(R * temp_k))*1E10;
+        D = d0 * exp( -Q /( R * temp_k)) * 1 * 10^10;
     end
-   % disp(K);
+   % disp(y);
 end</code>
+===== Drugi Program =====
+<code matlab>
+temp = 780;
+Q = 140000;
+R = 8.3144;
+temp_k = temp + 273;
+d0 = 0.000041;
+D=d0*exp(-Q /(R * temp_k))*1E10 ;
+dt = 0.01;
+dx = 0.1;
+pN = 0;
+time = 50;
+%Warunek początkowy NEUMANA
+if (D*dt)/dx^2<=0.5
+    K(6:100) = 0.02; %wypełniamy 6-100 elementy stężeniem o wielkości 0.2
+    K(1:5) = 0.67; %wypełniamy 1-5 elementy stężeniem o wielkości 0.67
+    ksi = 6; %ustawiamy na którym elemencie jest granica ksi
+    Clen = int32(time/dt); %obliczamy ilość kroków
+    C(1:Clen) = 0.0; %tworzymy tablicę wypełnioną zerami do której będzie
+    %obliczane stężenie całego materiału w każdym kroku czasowym
+    l = 1; %iterator powyższej tablicy
+    for k = 0:dt:time %zaczynamy proces
+        C(l) = trapz(K); %do tablicy zapisujemy stężenie w materiale
+        l = l + 1;
+        K = [K(1), K, K(length(K))]; %zwiększamy rozmiar tablicy z lewej i z prawej strony
+        K1 = K; %przepisujemy tablicę do tablicy pomocniczej (k+1)
+        K(ksi+1) = K(ksi); %powielamy element graniczny do obliczeń
+        for i = 2:1:ksi+1 %obliczamy stężenie do granicy
+           K1(i) = D*(dt/dx^2)*(K(i+1)+K(i-1)) + K(i)*(1 - 2*D*(dt/dx^2));
+           %obliczanie stężenia dla i-tego elementu
+        end
+        K = K1; %przepisanie tablicy spowrotem
+        K = K(2:1:length(K)-1); %zmniejszenie rozmiaru tablicy do stanu początkowego
+        if (mod(l,100) == 2) %wyświetlamy co 100 element na podstawie iteratora
+            plot(K) %wyświetlamy wykres stężenia w materiale
+            hold on
+        end
+        y=-0.004074074074074075*temp+3.7155555555555555;
+        %na podstawie wyznaczonego wzoru na prostą G-S wyliczamy stężenie
+        %dla danej temperatury
+        if K(ksi) >= y && ksi < length(K)-1 %jeśli stężenie jest wyższe
+            %niż G-S i ksi nie jest ostatnim elementem to zwiększamy granicę
+            ksi = ksi + 1;
+        end
+        temp = temp + pN * dt; %zwiększamy co krok temperaturę
+        temp_k = temp + 273;
+        D=d0*exp(-Q /(R * temp_k))*1E10; %ponownie obliczamy wsp. dyfuzji
+        if (D*dt)/dx^2>0.5 %sprawdzamy czy warunek NEUMANNA nadal jest spełniony
+            disp('warunek neumanna przestal byc spelniony dla temp');
+            disp(temp);
+            break;
+        end
+    end
+    figure;
+    plot(C); %wyświetlamy wykres zmiany stężenia co krok
+end
+</code>
+===== Drugi Program Maciek =====
+<code matlab>
+temp = 780;
+Q = 140000;
+R = 8.3144;
+temp_k = temp + 273;
+d0 = 0.000041;
+D=d0*exp(-Q /(R * temp_k))*1E10 ;
+dt = 0.01;
+dx = 0.1;
+pN = 0;
+time = 50;
+K(6:100) = 0.02;
+K(1:5) = 0.67;
+ksi = 6;
+Clen = int32(time/dt);
+C(1:Clen) = 0.0;
+%Neuman
+if (D*dt)/dx^2<=0.5
+    j = 1;
+    for k = 0:dt:time
+        C(l) = trapz(K); %Calka
+        j = j + 1;
+        K = [K(1), K, K(length(K))];
+        TempK = K;
+        K(ksi+1) = K(ksi);
+        for i = 2:1:ksi+1
+           TempK (i) = D*(dt/dx^2)*(K(i+1)+K(i-1)) + K(i)*(1 - 2*D*(dt/dx^2));
+        end
+        K = TempK;
+        K = K(2:1:length(K)-1);
+        if (mod(j,100) == 2)
+            plot(K)
+            hold on
+        end
+        y=-0.004074074074074075*temp+3.7155555555555555;
+        %zwiekszanie granicy
+        if K(ksi) >= y && ksi < length(K)-1
+            ksi = ksi + 1;
+        end
+        temp = temp + pN * dt;
+        temp_k = temp + 273;
+        D=d0*exp(-Q /(R * temp_k))*1E10;
+        if (D*dt)/dx^2>0.5
+            break;
+        end
+    end
+    figure;
+    plot(C);
+end
+</code>

WoGu Wiki

Narzędzia użytkownika

Narzędzia witryny

Różnice

Narzędzia strony