11. Προβλήματα οριακών συνθηκών

Τα προβλήματα οριακών συνθηκών ορίζονται σε έναν πεδίο υπολογισμού (domain) στο οποίο ισχύει μία γνωστή διαφορική εξίσωση ή σύστημα εξισώσεων. Απαραίτητη είναι και η γνώση των συνθηκών στα όρια του πεδίου υπολογισμού. Στόχος μας είναι να επιλύσουμε την διαφορική εξίσωση και να προσδιορίσουμε την διακύμανση των εξαρτημένων μεταβλητών στο εσωτερικό του πεδίου. Στην περίπτωση που είναι δύσκολο να εξάγουμε αναλυτική λύση της διαφορικής εξίσωσης, μπορούμε να χρησιμοποιήσουμε αριθμητικές μεθόδους για τον υπολογισμό της εξαρτημένης μεταβλητής σε συγκεκριμένα σημεία. Τέτοιες μέθοδοι είναι οι:

• Μέθοδος πεπερασμένων διαφορών

• Μέθοδος πεπερασμένων όγκων (FVM)

• Μέθοδος πεπερασμένων στοιχείων (FEM)

Με την εφαρμογή τους η διαφορική εξίσωση μετατρέπεται σε ένα σύστημα γραμμικών εξισώσεων, το οποίο μπορεί να λυθεί με τις μεθόδους που συζητήθηκαν ήδη. Εδώ θα συζητηθούν οι πρώτες δύο. Η μέθοδος των πεπερασμένων διαφορών είναι μια καθαρά μαθηματική προσέγγιση με την χρήση σειρών Taylor. Η μέθοδος πεπερασμένων όγκων αξιοποιεί τις πεπερασμένες διαφορές, αλλά δίνει έμφαση στην ορθή τήρηση ισοζυγίων σε στοιχειώδεις όγκους ελέγχου και στην ακριβή περιγραφή των οριακών συνθηκών. Το κλασσικό βιβλίο του Patankar [11] παρουσιάζει μια πολύ καλή εισαγωγή στην μέθοδο αυτή. Η μέθοδος πεπερασμένων όγκων είναι πολύ δημοφιλής σε ερευνητικά και εμπορικά πακέτα λογισμικού για μηχανικούς τα οποία αφορούν την επίλυση ισοζυγίων (ενέργειας, μάζας, ορμής, ηλεκτρικού φορτίου). Μάλιστα κάποια από αυτά είναι ανοιχτού κώδικα [12, 13].

11.1. Διακριτοποίηση

Αρχικά το πεδίο υπολογισμού χωρίζεται αυθαίρετα σε \(n+1\) μικρότερα κελιά με αρίθμηση \(i=0,n\). Αναλόγως την μέθοδο τα κελιά ονομάζονται και όγκοι ελέγχου (FVM) ή στοιχεία (FEM). Τα κελιά αυτά μπορούν να έχουν ομοιόμορφες ή ανομοιόμορφες διαστάσεις. Συνήθως επιλέγεται ένα ομοιόμορφο βήμα διακριτοποίησης. Για κάθε κελί πρέπει να οριστεί ένα κεντρικό (αντιπροσωπευτικό) σημείο, το οποίο συνήθως ταυτίζεται με το γεωμετρικό κέντρο, καθώς και τα όριά του.

Exercise 11.1

Διακριτοποιήστε το πεδίο υπολογισμού για μονοδιάστατο πρόβλημα με όρια στο \(x=0\) και \(x=L\).

Solution to Exercise 11.1

Επιλέγουμε σταθερό βήμα \(Δx=\frac{L}{n}\). Τα κέντρα των κελιών ορίζονται στα σημεία \(0,\Delta x,2\Delta x,\dots,L\) στα οποία αντιστοιχούν οι τιμές της εξαρτημένης μεταβλητής \(y_0,y_1,\dots,y_n\). Τα όρια των κελιών ορίζονται στις ενδιάμεσες αποστάσεις \(0,\frac{1}{2}\Delta x,\frac{3}{2}\Delta x,\dots,L-\frac{1}{2}\Delta x,L\).

• το κελί \(1\) θα έχει κεντρικό σημείο στο \(\Delta x\) και θα καλύπτει την περιοχή \([\frac{1}2\Delta x,\frac{3}2\Delta x]\)

• το οριακό κελί \(0\) θα έχει κεντρικό σημείο στο \(x=0\) (δεν είναι υποχρεωτικό να ταυτίζεται με το γεωμετρικό κέντρο) και θα καλύπτει την περιοχή \([0,\frac{1}2\Delta x]\) (μισό \(\Delta x\)). Αντιστοίχως και το τελευταίο κελί \(n\) στο άλλο όριο.

• Καθώς τα δύο οριακά κελιά είναι μισά, \(L=(n-1)\Delta x +2 \frac{\Delta x}2= n\Delta x\).

11.2. Πεπερασμένες διαφορές

Στην συνέχεια εφαρμόζουμε απλοποιητικές παραδοχές, έτσι ώστε να μετατρέψουμε την διαφορική εξίσωση σε αλγεβρική. Στην μέθοδο πεπερασμένων διαφορών, οι παράγωγοι αντικαθίστανται άμεσα με τις προσεγγιστικές σχέσεις.

Exercise 11.2

Να διατυπωθεί με πεπερασμένες διαφορές η διαφορική εξίσωση:

\[ \frac{d^2y}{dx^2}=0 \]

Η παραπάνω εξίσωση ισοδυναμεί με την εξίσωση Laplace σε μια διάσταση και βρίσκει εφαρμογή σε πλήθος φυσικών φαινομένων.

Solution to Exercise 11.2

Αρκεί να προσεγγιστεί η δεύτερη παράγωγος σύμφωνα με την Ενότητα 7.4 ως:

\[\left. \frac{d^2y}{dx^2}\right|_{x_i} \approx \frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x^2}\]

\[\frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x^2}=0\]

11.3. Πεπερασμένοι όγκοι

Στην μέθοδο των πεπερασμένων όγκων, εκφράζουμε πρώτα την διαφορική μας εξίσωση ως ισοζύγιο του φυσικού μεγέθους που διατηρείται. Κάθε όρος στην διαφορική εξίσωση πρέπει να εκφράζει την ροή του διατηρούμενου μεγέθους διαμέσου μίας πλευράς του κελιού ή την παραγωγή/κατανάλωση εντός κελιού. Έτσι μπορεί να χρειαστεί να ολοκληρώσουμε την αρχική εξίσωση για το κελί \(i\):

Exercise 11.3

Να διατυπωθεί με πεπερασμένους όγκους η διαφορική εξίσωση της Άσκησης 11.2:

Solution to Exercise 11.3

\[ \int_{x_i-\frac{\Delta x}2}^{x_i+\frac{\Delta x}2}\frac{d^2y}{dx}dx=0 \]

ή

\[ \left. \frac{dy}{dx}\right|_{x_i-\frac{\Delta x}2}- \left. \frac{dy}{dx}\right|_{x_i+\frac{\Delta x}2}=0 \]

Κάθε παράγωγος μπορεί να προσεγγιστεί με κεντρικές διαφορές:

\[\left. \frac{dy}{dx}\right|_{x_i+\frac{\Delta x}2} \approx \frac{y_{i+1} - y_{i}}{\Delta x}\]

Τελικά:

\[\frac{y_{i+1} - 2y_i + y_{i-1}}{\Delta x}=0\]

Αν απαλειφθεί ο παρονομαστής, η αλγεβρική εξίσωση ταυτίζεται με αυτή που προέκυψε με την μέθοδο πεπερασμένων διαφορών.

Για κάθε κελί \(i=0,n\) ορίζεται μία τιμή εξαρτημένης μεταβλητής \(y_i\) , επομένως προκύπτουν \(n+1\) άγνωστοι. Εφαρμόζοντας τις προσεγγιστικές σχέσεις στην διαφορική εξίσωση προκύπτουν \(n-1\) εξισώσεις για τα εσωτερικά κελιά.

11.4. Οριακές συνθήκες

Η διαδικασία της διακριτοποίησης εφαρμόζεται και στα δύο οριακά κελιά λαμβάνοντας υπόψη τις οριακές συνθήκες.
Έτσι προκύπτουν ακόμα δύο εξισώσεις (σύνολο \(n+1\)) και μπορεί να επιλυθεί το αλγεβρικό σύστημα για όλα τα κεντρικά σημεία.

• Σταθερής τιμής ή Dirichlet

  \[ y=f \]

• Neumann

  \[ \frac{dy}{dx}=f \]

• Συνδυασμοί των δύο παραπάνω

11.5. Εφαρμογή

Exercise 11.4

Η μετάδοση θερμότητας εντός ενός μονοδιάστατου τοιχώματος σε συνθήκες μόνιμης κατάστασης δίνεται από την διαφορική εξίσωση:

\[ -\frac{d}{dx}\left(k\frac{dT}{dx}\right)=0 \]

με τις παρακάτω οριακές συνθήκες σε πρακτικές εφαρμογές:

\[\begin{split} \left. {-k\frac{dT}{dx}}\right|_{0}=h_{in}(T_{in}-T_{x=0})\\ \left. -k\frac{dT}{dx}\right|_{L}=h_{out}(T_{x=L}-T_{out}) \end{split}\]

με

• \(k\) τον συντελεστή θερμικής αγωγιμότητας του υλικού του τοιχώματος

• \(h_{in}\) τον συντελεστή συναγωγής στο εσωτερικό

• \(h_{out}\) τον συντελεστή συναγωγής στο εξωτερικό

• \(L\) το πάχος του τοιχώματος

Solution to Exercise 11.4

Στην περίπτωση των πεπερασμένων διαφορών αντικαθιστούμε απευθείας την δεύτερη παράγωγο ως:

\[\left. \frac{d^2T}{dx^2}\right|_{x_i} \approx \frac{T_{i+1} - 2T_i + T_{i-1}}{\Delta x^2}\]

Στην περίπτωση των πεπερασμένων όγκων ολοκληρώνουμε την εξίσωσή μας ως προς \(x\), οπότε:

\[ -\left.k\frac{dT}{dx}\right|_{x_i+\frac{\Delta x}2}+\left.k\frac{dT}{dx}\right|_{x_i-\frac{\Delta x}2}=0 \]

και στην συνέχεια αντικαθιστούμε τις πρώτες παραγώγους:

\[\left. \frac{dT}{dx}\right|_{x_i-\frac{\Delta x}2} \approx \frac{T_{i} - T_{i-1}}{\Delta x}\]

\[\left. \frac{dT}{dx}\right|_{x_i+\frac{\Delta x}2} \approx \frac{T_{i+1} - T_i}{\Delta x}\]

Και στις δύο περιπτώσεις προκύπτει η ίδια αλγεβρική εξίσωση που εφαρμόζεται σε όλα τα εσωτερικά σημεία (\(n-1\)):

\[-\frac{k}{\Delta x}T_{i+1} + 2\frac{k}{\Delta x}T_i - \frac{k}{\Delta x}T_{i-1}=0\]

Άλλες δύο εξισώσεις (σύνολο \(n+1\)) προκύπτουν από τις οριακές συνθήκες, εφαρμόζοντας εμπρός και πίσω διαφορές σύμφωνα με τη Ενότητα 7.2:

\[-k\frac{T_1-T_0}{\Delta x} =h_{in}(T_{in}-T_0)\]

\[-k\frac{T_{n}-T_{n-1}}{\Delta x} =h_{out}(T_{n}-T_{out})\]

ή

\[\left(\frac{k}{\Delta x}+h_{in}\right)T_0 - \frac{k}{\Delta x}T_1 =h_{in}T_{in}\]

\[- \frac{k}{\Delta x}T_{n-1} +\left(\frac{k}{\Delta x}+h_{out}\right)T_{n} =h_{out}T_{out}\]

Έτσι προκύπτει το τριδιαγώνιο γραμμικό σύστημα:

\[\begin{split} \begin{bmatrix} \frac{k}{Δx}+h_{in} & -\frac{k}{Δx} & 0 & {\ldots} & 0 \\ -\frac{k}{Δx} & 2\frac{k}{Δx} & -\frac{k}{Δx} & 0 & {\ldots} \\ 0 & -\frac{k}{Δx} & 2\frac{k}{Δx} & -\frac{k}{Δx} &{\ldots} \\ {\ldots} &{\ldots} &{\ldots} & {\ldots}&{\ldots} \\ 0 & {\ldots}& 0 & -\frac{k}{Δx} &\frac{k}{Δx}+h_{out} \\ \end{bmatrix} \begin{bmatrix} T_0\\ T_1\\ T_2\\ {\ldots}\\ T_{n}\\ \end{bmatrix} = \begin{bmatrix} h_{in}T_{in}\\ 0\\ 0\\ {\ldots}\\ h_{out}T_{out}\\ \end{bmatrix} \end{split}\]

Για την επίλυση του γραμμικού συστήματος χρησιμοποιούνται οι μέθοδοι του Κεφαλαίου 4.