CFD-lösare

[AndersF1]

Jag hittade ett gammalt matlab-program som jag skrev för en massa år sedan när jag skulle implementera s.k. "dual timestepping" i en CFD-lösare. Den är bara i 1D(till skillnad ifrån verkligehetens 3D) och endast för inviskösa problem, men beskriver ändå ganska bra hur en enkel CFD-lösare ser ut och fungerar. Det finns lite "dual timestepping"-pryttlar med, men inte så jättemycket så att det stör.

Programmet: http://www.nada.kth.se/~ayt/programs/eulsolv.m

Nu har ju inte alla matlab hemma, men det borde även fungera med nån av gratisvarianterna som finns. Kan möjligen bli problem med grafiken. Jag har dock ingen länk till något sånt program. Någon?


Den viktigaste delen där det mesta av lösningen räknas ut:

Kod: Markera allt

 for nr = 1:nrmax,
%
% Inviscid, convective fluxes
%
         f(nfirst:nlast,1) = (w(nfirst+1:nlast+1,2)-w(nfirst-1:nlast-1,2))*0.5;
         f(nfirst:nlast,2) = ((w(nfirst+1:nlast+1,2).*u(nfirst+1:nlast+1) + ...
                               p(nfirst+1:nlast+1)) - ... 
                              (w(nfirst-1:nlast-1,2).*u(nfirst-1:nlast-1) + ...
                               p(nfirst-1:nlast-1)))*0.5;
         f(nfirst:nlast,3) = ((w(nfirst+1:nlast+1,3) + p(nfirst+1:nlast+1)).*...
                               u(nfirst+1:nlast+1) - ...
                              (w(nfirst-1:nlast-1,3) + p(nfirst-1:nlast-1)).*...
                               u(nfirst-1:nlast-1))*0.5;
%
% Update artificial dissipation for the first 2 R-K stages
%
         if (nr < 3),
            mu = abs((p(nfirst:nlast+2) - 2*p(nfirst-1:nlast+1) + ...
                      p(nfirst-2:nlast))./(p(nfirst:nlast+2) + ...
                      2*p(nfirst-1:nlast+1) + p(nfirst-2:nlast)));
%
            eps2 = art2*max(mu(1:n+1),mu(2:n+2));
            eps4 = max(zeros(size(eps2)),art4 - eps2);
%
            fd(nfirst:nlast,1) = ...
              0.5*(lambda(nfirst+1:nlast+1) + lambda(nfirst:nlast)).* ...
                (eps2(2:n+1).*(w(nfirst+1:nlast+1,1)-w(nfirst:nlast,    1)) -...
                 eps4(2:n+1).*(w(nfirst+2:nlast+2,1)-3*w(nfirst+1:nlast+1,1)+...
                 3*       w(nfirst:nlast    ,1) -   w(nfirst-1:nlast-1,1))) -...
              0.5*(lambda(nfirst:nlast) + lambda(nfirst-1:nlast-1)).* ...
                (eps2(1:n).*  (w(nfirst:nlast,    1)-w(nfirst-1:nlast-1,1)) -...
                 eps4(1:n).*  (w(nfirst+1:nlast+1,1) - 3*w(nfirst:nlast,1)  +...
                 3*           w(nfirst-1:nlast-1,1) -   w(nfirst-2:nlast-2,1)));
            fd(nfirst:nlast,2) = ...
              0.5*(lambda(nfirst+1:nlast+1) + lambda(nfirst:nlast)).* ...
                (eps2(2:n+1).*(w(nfirst+1:nlast+1,2) - w(nfirst:nlast,  2)) -...
                 eps4(2:n+1).*(w(nfirst+2:nlast+2,2)-3*w(nfirst+1:nlast+1,2)+...
                 3*       w(nfirst:nlast    ,2) -   w(nfirst-1:nlast-1,2))) -...
              0.5*(lambda(nfirst:nlast) + lambda(nfirst-1:nlast-1)).* ...
                (eps2(1:n).* (w(nfirst:nlast, 2) -   w(nfirst-1:nlast-1,2)) -...
                 eps4(1:n).*(w(nfirst+1:nlast+1,2)-3*w(nfirst:nlast,    2)  +...
                 3*          w(nfirst-1:nlast-1,2) -   w(nfirst-2:nlast-2,2)));
            fd(nfirst:nlast,3) = ...
              0.5*(lambda(nfirst+1:nlast+1) + lambda(nfirst:nlast)).* ...
                (eps2(2:n+1).*(w(nfirst+1:nlast+1,3) - w(nfirst:nlast, 3))  -...
                 eps4(2:n+1).*(w(nfirst+2:nlast+2,3)-3*w(nfirst+1:nlast+1,3)+...
                 3*       w(nfirst:nlast    ,3) -   w(nfirst-1:nlast-1,3))) -...
               0.5*(lambda(nfirst:nlast) + lambda(nfirst-1:nlast-1)).* ...
                 (eps2(1:n).*  (w(nfirst:nlast,3) -  w(nfirst-1:nlast-1,3)) -...
                  eps4(1:n).* (w(nfirst+1:nlast+1,3) - 3*w(nfirst:nlast,3)  +...
                  3*         w(nfirst-1:nlast-1,3) -   w(nfirst-2:nlast-2,3)));
        end
%
        if (timeaccurate & dualt),
%
% Add the source term and flux for the dual time stepping to the
% convective flux
%
           if (outerstep > 1),
% 2nd order backward difference
              fac = 1.5/dtout;
           else
% 1st order backward Euler
              fac = 1/dtout;
           end
%
           f(nfirst:nlast,:) = f(nfirst:nlast,:) + fac*w(nfirst:nlast,:).*...
                               vol(nfirst:nlast,:) + duals(nfirst:nlast,:);
        end
%
% Update the state vector
%
        w(nfirst:nlast,1) = wt(nfirst:nlast,1)-nrcoef(nr)*dt(nfirst:nlast).*(...
                            f(nfirst:nlast,1) - fd(nfirst:nlast,1));
        w(nfirst:nlast,2) = wt(nfirst:nlast,2)-nrcoef(nr)*dt(nfirst:nlast).*(...
                            f(nfirst:nlast,2) - fd(nfirst:nlast,2));
        w(nfirst:nlast,3) = wt(nfirst:nlast,3)-nrcoef(nr)*dt(nfirst:nlast).*(...
                            f(nfirst:nlast,3) - fd(nfirst:nlast,3));
%
% Update the boundary condition
%
        w(nlast+1,1) = w(nlast,1);
        w(nlast+2,1) = w(nlast-1,1);
        w(nlast+1,2) = -w(nlast,2);
        w(nlast+2,2) = -w(nlast-1,2);
        w(nlast+1,3) = w(nlast,3);
        w(nlast+2,3) = w(nlast-1,3);
%
% Update pressure, velocity and square of sound velocity
%
        w(:,1) = abs(w(:,1));
        u = w(:,2)./w(:,1);
        p = gammafac*(w(:,3) - 0.5*u.*u.*w(:,1));
        c2 = gamma*p./w(:,1);
%
      end

[AndersF1]

Ja just det, en liten beskrivning för den intresserade:
Rumsdiskretisering med 2:a ordningens centralt schema samt 2:a och 4:e ordningens artificiell viskositet av Jameson-typ. Tidsdiskretisering med explicit Runge-Kutta i variabelt antal steg, men högst andra ordningens nogrannhet.

Fallet som körs är ett stötrör, där man får en plötslig höjning av trycket i ena änden, den till höger, p.g.a. stoppat flöde.


_________________
"Ingen är perfekt, själv är jag känslig för drag"
-Oscar Wilde

<font size=-1>[ Detta Inlägg ändrades av: AndersF1 den 2002-03-08 18:09 ]</font>

[Granis]

Intressant Anders, den måste jag testa.

Jag har en länk till ett Matlabderivat som heter MatrixCalc och körs faktiskt på PalmOS. http://www.mtrxcal.com/mtrxcal/index.htm

Hmmm, jag har Nummetenta imorgon lördag. Kan man försöka använda Palm'en som miniräknare på tentan kanske.

[Hood]

On 2002-03-08 18:08, AndersF1 wrote:
Rumsdiskretisering med 2:a ordningens centralt schema samt 2:a och 4:e ordningens artificiell viskositet av Jameson-typ. Tidsdiskretisering med explicit Runge-Kutta i variabelt antal steg, men högst andra ordningens nogrannhet.

Ja, det är inte många svenska ord i den meningen som en vanlig nasare som jag fattar...

// J.

[Gäst]

"Runge-Kutta"? Men AndersF1 då, vet du inte att det finns barn som läser här?

[Robert.]

Jari: Förnatt inte tala om: "Fallet som körs är ett stötrör, där man får en plötslig höjning av trycket i ena änden, den till höger, p.g.a. stoppat flöde. "

Illa om man skulle få stoppat flöde i stötröret förresten

[Gäst]

Absolut!
Jag körde folk i natt som skulle ut och supa med Benny. Var det vänner till dig också och varför var du inte med ute?
Du kanske ville avvakta till ikväll när Forza går åstad?