Różnice

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

--- studia:magisterskie:1sem:modelowanie_procesow_fizycznych [2016/04/20 09:20]
149.156.112.6 [Zajęcia 3]
+++ studia:magisterskie:1sem:modelowanie_procesow_fizycznych [2016/06/01 08:49] (aktualna)
91.227.212.3 [Drugi Program Maciek]
@@ Linia 93: / Linia 93: @@
 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