מצאתם טעות? שלחו הודעה קצרה. גם אם זה רק שגיעת כתיב קטנה. תודה לינאי וגיל ששיכנעו אותי להוסיף את זה...
FEM1_007 משוואת פואסון הדו-ממדית
• זמן קריאה: 38 דק'
צורה חלשה של משוואת פואסון
נפתח כעת הצורה החלשה של משוואת פואסון. בצורתו החזקה, למשל עבור בעיית חום:
לשם גזירת הצורה החלשה, נכפיל את המשוואה הדיפרנציאלית בפונקציית מבחן שרירותית (חלקה מספיק ומתאפסת על שפת דיריכלה , כלומר ), ונבצע אינטגרציה על התחום :
ניישם אינטגרציה בחלקים (לפי משפט הדיברגנץ, או זהות גרין הראשונה) לאיבר הראשון. נשתמש בזהות . נבחר ו-. אז .
לכן, האיבר הראשון הופך ל:
על האיבר השני באגף ימין נפעיל את משפט הדיברגנץ:
כאשר הוא וקטור הנורמל היוצא מהשפה . הביטוי לאיבר הראשון בצורה האינטגרלית המקורית הוא:
נציב זאת חזרה למשוואה האינטגרלית הכוללת:
האינטגרל המשטחי על השפה מתפרק לשני חלקים. על , פונקציית המבחן מתאפסת (), ולכן האינטגרל על הוא אפס. נשארנו עם האינטגרל על . נשתמש בתנאי השפה הנתון על : .
לפיכך, האינטגרל המשטחי הוא:
הצבה חזרה וסידור מחדש מובילים לצורה החלשה של משוואת פואסון:
יש למצוא המקיימת על , כך שלכל פונקציית מבחן קבילה (עם על ):
כעת נציג את קירוב האלמנטים הסופיים עבור הצורה החלשה של משוואת פואסון. נקרב את פונקציית הטמפרטורה ואת פונקציית המבחן באמצעות פונקציות צורה (או פונקציות בסיס) . יהיו הערכים הצומתיים של הטמפרטורה (דרגות החופש של הפתרון) ו- מקדמים שרירותיים עבור פונקציית המבחן.
הקירובים הם:
כאשר פונקציות הצורה נבחרות כך ש- על שפת דיריכלה . הנגזרות (גרדיאנטים) יהיו:
נציב קירובים אלו לצורה החלשה. ניתן להוציא את הסכום על ואת המקדמים מהאינטגרלים המתאימים:
מאחר שמשוואה זו צריכה להתקיים לכל בחירה שרירותית של מקדמי פונקציית המבחן (הכפופים לתנאי השפה ההומוגניים על ), הביטוי שבסוגריים המרובעים הגדולים חייב להתאפס לכל . הדבר מוביל למערכת של משוואות אלגבריות לינאריות מהצורה:
כאשר הוא וקטור המקדמים הצומתיים הנעלמים של הטמפרטורה . איברי המטריצה (מטריצת הקשיחות או המוליכות) והווקטור (וקטור העומס) ניתנים על ידי:
מיפוי אלמנט משולש לינארי כללי
בסעיף זה נפתח את מטריצת הקשיחות האלמנטרית עבור אלמנט משולש לינארי כללי. הגישה המודרנית בשיטת האלמנטים הסופיים מתבססת על מיפוי של כל אלמנט פיזיקלי אל אלמנט אב (master element) סטנדרטי, עליו מוגדרות פונקציות הצורה ומתבצעת האינטגרציה הנומרית.
הביטוי הכללי לאיבר במטריצת הקשיחות הוא:
כאשר הוא מקדם דיפוזיה/מוליכות, הן פונקציות הצורה במונחי קואורדינטות גלובליות של האלמנט הפיזיקלי .
אלמנט המסטר ופונקציות הצורה
נבחר אלמנט אב משולש במערכת קואורדינטות מקומית . קונפיגורציה נפוצה היא משולש ישר זווית עם צמתים בנקודות:
צומת מקומי 1:
צומת מקומי 2:
צומת מקומי 3:
אלמנט מסטר למשולש.
פונקציות הצורה הלינאריות על אלמנט המסטר, המסומנות כ- , הן:
לכל פונקציית צורה יש ערך בצומת המקומי ו- בשני הצמתים האחרים. שטח אלמנט המסטר הוא .
הערה:
שימו לב שבתרגול מגדירים את הצמתים בסדר שונה, ולכן , ו- מוגדרים בסדר שונה.
נוכל להגדיר:
מיפוי איזופרמטרי מאלמנט המסטר לאלמנט הפיזיקלי
אלמנט משולש פיזיקלי במרחב הגלובלי מוגדר על ידי קואורדינטות צמתיו: , , ו-. אנו מניחים שהצומת הגלובלי מתאים לצומת המקומי .
המיפוי (טרנספורמציה) מקואורדינטות מקומיות לקואורדינטות גלובליות ניתן על ידי:
זוהי טרנספורמציה לינארית, הממפה קווים ישרים לקווים ישרים.
אם נגדיר:
נוכל לרשום את המיפוי באופן הבא:
כאשר .
טרנספורמציה של נגזרות ומטריצת היעקוביאן
מטריצת היעקוביאן של הטרנספורמציה, , מוגדרת באופן הבא:
אם נגדיר:
נוכל לרשום את מטריצת היעקוביאן בצורה הבאה:
עבור פונקציות הצורה המשולש הלינארי:
הנגזרות של פונקציות הצורה הגלובליות לפי הקואורדינטות הגלובליות מתקבלות מהנגזרות של פונקציות הצורה המקומיות באמצעות כלל השרשרת:
נוכל לרשום את משוואה זאת באופן מטריצי:
כאשר מכילה את הנגזרות הגלובליות של כל פונקציות הצורה.
חישוב אינטגרלים ומטריצת קשיחות אלמנטרית
תרומת האלמנט למטריצת הקשיחות הגלובלית, לדוגמה עבור משוואת פואסון עם מקדם , היא:
האינטגרל על האלמנט הפיזיקלי מומר לאינטגרל על אלמנט המסטר :
נוכל לרשום את חישובה באופן מטריצי:
כאשר מטריצת היעקוביאן וההופכית שלה מחושבות לפי המיפוי הגיאומטרי בין האלמנט הגלובלי לאלמנט הייחוס, ואינטגרציית גאוס מיושמת כפי שהוסבר באינטגרציית גאוס עבור משולשים.
אינטגרציה נומרית - אינטגרציית גאוס עבור משולשים
עבור אלמנטים משולשיים לינאריים, האינטגרציה מתבצעת באמצעות אינטגרציית גאוס על התחום המשולש. עבור אינטגרנד עד דרגה 2, נקודת גאוס אחת במרכז מספיקה:
נקודה אחת (במרכז המשולש):
שלוש נקודות (במרכז הצלעות):
כאשר נקודות הגאוס הן:
לכן, מטריצת הקשיחות האלמנטרית מחושבת כ:
כאשר הם משקלי הגאוס ו- הוא מספר נקודות הגאוס.
מיפוי אלמנט מרובע ביליניארי כללי
בדומה לאלמנט המשולש, גם עבור אלמנטים מרובעים, נהוג להשתמש באלמנט מסטר סטנדרטי לטובת פיתוח פונקציות הצורה וביצוע אינטגרציה נומרית.
הביטוי הכללי לאיבר במטריצת הקשיחות הוא:
כאשר הוא מקדם דיפוזיה/מוליכות, הן פונקציות הצורה במונחי קואורדינטות גלובליות של האלמנט הפיזיקלי .
אלמנט המסטר ופונקציות הצורה
אלמנט האב הסטנדרטי עבור אלמנט מרובע ביליניארי הוא ריבוע במערכת קואורדינטות מקומית , כאשר ו- . הצמתים המקומיים ממוספרים בדרך כלל נגד כיוון השעון, החל מהפינה השמאלית התחתונה:
צומת מקומי 1:
צומת מקומי 2:
צומת מקומי 3:
צומת מקומי 4:
פונקציות הצורה הביליניאריות עבור אלמנט זה הן:
לכל פונקציית צורה יש ערך בצומת המקומי ואפס בשאר שלושת הצמתים. פונקציות אלו מאפשרות אינטרפולציה ביליניארית של שדה (לדוגמה, תזוזה או טמפרטורה) בתוך האלמנט, מהצורה . חשוב לציין כי סכום פונקציות הצורה בכל נקודה בתוך האלמנט הוא אחד: .
נוכל להגדיר:
מיפוי איזופרמטרי מאלמנט המסטר לאלמנט הפיזיקלי
המיפוי (טרנספורמציה) מקואורדינטות מקומיות לקואורדינטות גלובליות של אלמנט מרובע פיזיקלי כלשהו ניתן על ידי שימוש באותן פונקציות צורה (מיפוי איזופרמטרי):
כאשר הן הקואורדינטות הגלובליות של צומת באלמנט הפיזיקלי.
אם נגדיר:
נוכל לרשום את המיפוי באופן הבא:
כאשר .
טרנספורמציה של נגזרות ומטריצת היעקוביאן
מטריצת היעקוביאן של הטרנספורמציה, , מוגדרת באופן זהה למקרה המשולש:
אם נגדיר:
נוכל לרשום את מטריצת היעקוביאן בצורה הבאה:
עבור פונקציות הצורה הביליניאריות:
אולם, בניגוד לאלמנט משולש לינארי, עבור אלמנט מרובע ביליניארי, רכיבי מטריצת היעקוביאן (ולכן גם הדטרמיננטה שלה, ) הם בדרך כלל פונקציות של ואינם קבועים על פני האלמנט.
הנגזרות של פונקציות הצורה הגלובליות לפי הקואורדינטות הגלובליות מתקבלות מהנגזרות של פונקציות הצורה המקומיות באמצעות כלל השרשרת:
נוכל לרשום את משוואה זאת באופן מטריצי:
הנגזרות המקומיות ו- מחושבות ישירות מפונקציות הצורה. לדוגמה:
וכך הלאה עבור ו- .
חישוב אינטגרלים ומטריצת קשיחות אלמנטרית
תרומת האלמנט למטריצת הקשיחות הגלובלית היא:
האינטגרל על האלמנט הפיזיקלי מומר לאינטגרל על אלמנט המסטר :
נוכל לרשום את חישובה באופן מטריצי:
מכיוון שהנגזרות , (התלויות ב-) והיעקוביאן הם פונקציות של , האינטגרנד בדרך כלל אינו פולינום פשוט או קבוע. לכן, יש לבצע את האינטגרציה באופן נומרי.
אינטגרציה נומרית - אינטגרציית גאוס עבור מרובעים
עבור אלמנטים מרובעים ביליניאריים על התחום , האינטגרציה מתבצעת באמצעות אינטגרציית גאוס דו-ממדי.
נקודה אחת (במרכז):
נקודות גאוס:
כאשר נקודות הגאוס והמשקלים הם:
הנקודות הן:
לכן, מטריצת הקשיחות האלמנטרית מחושבת כ:
מכיוון שהנגזרות , (התלויות ב-) והיעקוביאן הם פונקציות של , יש להעריך אותם בכל נקודת גאוס בנפרד.
תוצאת האינטגרציה הנומרית עבור כל זוג נותנת את איברי מטריצת הקשיחות האלמנטרית , שהיא מטריצה בגודל עבור אלמנט מרובע עם 4 צמתים.
הרכבת מטריצת הקשיחות הגלובלית
לאחר שפיתחנו את הצורה החלשה של המשוואה הדיפרנציאלית ואת הביטויים לאיברי מטריצת הקשיחות והווקטור העומס הגלובליים, וראינו כיצד לחשב מטריצת קשיחות אלמנטרית עבור אלמנט ספציפי, השלב הבא הוא הרכבת המערכת הגלובלית של המשוואות האלגבריות .
האינטגרלים המגדירים את ו- על פני כל התחום (ועל שפת ניומן/רובין ) מחושבים כסכום של תרומות מכל האלמנטים הבודדים ברשת. כלומר:
כאשר הוא מספר האלמנטים הכולל ברשת. ו- מייצגים את התרומה של אלמנט לאיבר במטריצה הגלובלית ולאיבר בווקטור העומס הגלובלי, בהתאמה. תרומות אלו נובעות מהאינטגרלים על תחום האלמנט ועל חלק השפה השייך לאלמנט :
בביטויים אלו, ו- הן פונקציות הבסיס הגלובליות המשויכות לצמתים הגלובליים ו-. כאשר מבצעים אינטגרציה על אלמנט , פונקציות בסיס אלו הופכות לפונקציות הצורה המקומיות של האלמנט אם הצמתים ו- שייכים לאלמנט ; אחרת, תרומתן מאלמנט זה היא אפס.
תהליך ההרכבה מתבצע באופן שיטתי:
אתחול: יוצרים מטריצת קשיחות גלובלית ווקטור עומס גלובלי . מאתחלים את כל איבריהם לאפס.
לולאה על האלמנטים: עוברים על כל אלמנט ברשת, מ- ועד :
חישוב מטריצות אלמנטריות: עבור האלמנט הנוכחי , מחשבים את מטריצת הקשיחות האלמנטרית ואת וקטור העומס האלמנטרי . גודלה של הוא , כאשר הוא מספר הצמתים באלמנט (למשל, לאלמנט משולש לינארי).
זיהוי אינדקסים גלובליים: באמצעות טבלת הקישוריות (כפי שתוארה בסעיף [[FEM1_006 שיטת אלמנטים סופיים בתלת-ממד#|יצירת רשת ופונקציות קישוריות]]), מאתרים את מספרי הצמתים הגלובליים המתאימים לכל צומת מקומי של האלמנט . נסמן את הפונקציה הממפה כ-, המחזירה את האינדקס הגלובלי של הצומת המקומי באלמנט .
פיזור והוספה (Assembly): מוסיפים את איברי המטריצה והווקטור האלמנטריים למיקומים המתאימים במטריצה ובווקטור הגלובליים:
עבור מ- עד , כאשר (אינדקס שורה/רכיב גלובלי):
עבור מ- עד , כאשר (אינדקס עמודה גלובלי):
המערכת הגלובלית: לאחר שעוברים על כל האלמנטים, המטריצה והווקטור מכילים את התרומות מכל חלקי התחום והשפה הרלוונטית. המערכת מייצגת את הקירוב של הבעיה המקורית על פני כל הרשת.
יש לציין כי המטריצה הגלובלית היא לרוב דלילה (sparse), כלומר רוב איבריה הם אפס. זאת מכיוון שכל צומת גלובלי קשור רק למספר קטן של צמתים אחרים (אלו החולקים איתו אלמנטים). ביישומים מעשיים, משתמשים בשיטות אחסון ופתרון יעילות המנצלות דלילות זו.
לאחר הרכבת המערכת הגלובלית, השלב הבא, לפני פתרון המערכת, הוא החלת תנאי שפה מסוג דיריכלה, אשר משנים את המטריצה ואת הווקטור כדי לכפות את הערכים הידועים של בצמתים הרלוונטיים.
חישוב תנאי שפה מסוג רובין
תנאי שפה מסוג רובין הם מהצורה הכללית:
כאשר הוא מקדם המוליכות, , הוא השטף הנכפה, ו- הוא וקטור הנורמל היוצא מהשפה.
תנאי רובין תורמים לשני חלקים בצורה החלשה:
תרומה למטריצת הקשיחות:
תרומה לווקטור העומס:
חישוב אינטגרלים שפתיים עבור אלמנטים לינאריים
עבור קטע שפה לינארי המחבר בין שני צמתים ו- באורך , פונקציות הצורה הלינאריות הן:
כאשר הוא הפרמטר לאורך הקטע ().
מטריצת הקשיחות בשפה
עבור מקדם מעבר קבוע , מטריצת הקשיחות השפתית עבור קטע היא:
חישוב האינטגרלים:
לכן, מטריצת הקשיחות השפתית עבור קטע היא:
וקטור העומס בשפה
עבור שטף נכפה , ווקטור העומס השפתי הוא:
מקרה מיוחד - שטף קבוע: אם קבוע לאורך הקטע:
מקרה כללי - שטף משתנה לינארית: אם משתנה לינארית לאורך הקטע עם ערכים ו- בצמתים:
חישוב האינטגרלים:
לכן:
הרכבה במטריצה הגלובלית
לאחר חישוב המטריצות והווקטורים השפתיים עבור כל קטע שפה, יש להרכיב אותם במטריצת הקשיחות הגלובלית ובווקטור העומס הגלובלי בהתאם לאינדקסים הגלובליים של הצמתים.
דוגמה:
עבור קטע שפה בין צמתים גלובליים ו-, תרומת הקשיחות השפתית תתווסף למיקומים , , , ו- במטריצה הגלובלית.
תרומת תנאי רובין (אינטגרל שפה על ):
נמצא לפי חישוב תנאי שפה מסוג רובין. התרומה היא . הקטע הוא בין צומת לצומת . אורך הקטע . פונקציות הצורה הלינאריות על קטע זה הן ו- . מטריצת התרומה היא:
(שימו לב לסדר השורות והעמודות בהתאם לצמתים הגלובליים 3 ו-4):
מטריצת הקשיחות הגלובלית הכוללת:
וקטור העומס הגלובלי כולל את התרומה מתנאי הרובין על . מכיוון ש-, התרומה היחידה היא מהאינטגרל השפתי .
על (בין צמתים 3 ו-4), . על ו- (תנאי נוימן הומוגניים), .
האינטגרל עבור פונקציית בסיס לינארית על קטע באורך הוא עבור כל אחד משני צמתי הקטע.
הקטע הוא באורך . לכן:
כל שאר איברי וקטור העומס הם אפס ().
לכן, .
החלת תנאי שפה דיריכלה:
נתון , ולכן ו-. נותר למצוא את .
מערכת המשוואות הגלובלית היא . ניקח את המשוואות השלישית והרביעית:
הפתרון עבור הוא ו-. לכן הפתרון הסופי:
תרגיל 2
נתונה הבעיה הבאה:
כאשר:
תחום הבעיה הנתונה.
פתרו בעזרת שני אלמנטים משולשיים לינאריים.
פתרון:
הצורה החלשה הכללית למשוואת פואסון היא:
כאשר איברי המטריצה והווקטור הם:
בבעיה הנתונה, , ולכן מקדם המוליכות , (אין איבר תלוי ב- במשוואה הדיפרנציאלית עצמה), ואין מקור חימום פנימי .
תנאי השפה הם:
על : (תנאי דיריכלה).
על : . זהו תנאי רובין מהצורה . מאחר ו-, אז ו-.
לכן, עבור הבעיה הספציפית שלנו:
הצמתים הגלובליים הם (נסמן אותם בהתאמה לאינדקסים בוקטור הנעלמים ):
נחלק את התחום לשני אלמנטים משולשיים:
אלמנט : צמתים גלובליים .
אלמנט : צמתים גלובליים .
שטח כל אלמנט הוא .
פירוק התחום לשני אלמנטים משולשיים.
עבור אלמנט משולש לינארי כללי, מטריצת הקשיחות האלמנטרית מחושבת באמצעות המיפוי האיזופרמטי כפי שהוסבר במיפוי אלמנט משולש לינארי כללי:
כאשר:
הוא מטריצת הנגזרות של פונקציות הצורה ביחס לקואורדינטות הייחוס
הוא מטריצת היעקוביאן המקשרת בין הקואורדינטות הגלובליות והייחוס
הוא הופכית מטריצת היעקוביאן
הוא מטריצת הנגזרות של פונקציות הצורה ביחס לקואורדינטות הגלובליות
אלמנט 1 (צמתים גלובליים ):
הקואורדינטות של הצמתים הגלובליים המרכיבים את אלמנט הן: . שטח האלמנט .
מיפוי לאלמנט הייחוס:
נמפה את האלמנט לאלמנט משולש ייחוס עם קואורדינטות ברילצנטרי כאשר :
חישוב מטריצת היעקוביאן :
מטריצת הנגזרות של פונקציות הצורה ביחס לקואורדינטות הייחוס:
מטריצת הקואורדינטות:
מטריצת היעקוביאן:
הדטרמיננטה: , לכן .
ההופכית:
חישוב מטריצת :
זה אומר:
חישוב מטריצת הקשיחות :
עבור אינטגרציית גאוס עם נקודה אחת במרכז הכובד , המשקל הוא :