FEM1_003 מימוש אלמנט סופי בחד-ממד

מבוא

שיטות קלאסיות לפתרון בשיטת גלרקין בונות קירובים מפונקציות קבילות קינמטית (kinematically admissible functions), שהן פונקציות שמקיימות את תנאי השפה על הזזה (דהיינו תנאי דיריכלה) מלכתחילה. שתי בעיות מרכזיות עולות:

1. לפעמים זה די קשה למצוא פונקציות קבילות קינמטית בכל המרחב הרלוונטי.
2. אם כן מצאנו פונקציות כאלה, הן מובילות למערכת משוואות מצומדות ומסובכות.

ניתן להתעלות על בעיות מסוג זה בעזרת העובדה שקירובים מקומיים יכולים לספק פתרונות בדיוק גבוה, ובאותו הזמן להוביל למערכת משוואות שיש להן מבנים מתמטיים נוחים לחישוב, גם בממדים גבוהים. קירובים חלקיים אלו, או “אלמנטיים” אלו, כבר הוכרו לפני לפחות 60 שנה ע”י קורנט (Courant). היה מגוון של שיטות קירוב דומות לפתרון משוואות פיזיקליות. הכי נפוצה ביניהן היא שיטת האלמנטים הסופיים, FEM. המאפיין המרכזי של שיטה זו היא חילוק המרחב באופן שיטתי להרכבה של תת-תחומים בדידים, או, “אלמנטים”, ואז לקרב את הפתרון של כל אחד מחלקים אלו באופן שמצמד אותם כדי ליצור פתרון גלובלי תקף לאורך כל התחום. התהליך פותח באופן כזה שהמשוואות האלגבריות יהיו ניתנות לחישוב (במחשב), ויעילות מבחינת הזיכרון שהן משתמשות בו.

קירוב אלמנט סופי

נביט בצורה החלשה של הבעיה שראינו בפרק קודם (משוואה ):

עלינו למצוא (הפתרון שלנו), כאשר (תנאי שפה כלשהו מסוג דיריכלה), כך שלכל שמקיים וגם , מתקיים:

אנו מקרבים את ע”י

אם נבחר עם אותם הפונקציות קירוב, אבל עם קומבינציה לינארית שונה

אז אנו יכולים לרשום:

מאחר וה--ים שרירותיים, גם שרירותיים (חוץ מהקצוות, שם תנאי השפה קובעים את ), כך ש:

כאשר:

את אנו יכולים גם לרשום בצורה מטריציונית:

כאשר היא מטריצה , ש-, באופן לא מפתיע, קוראים לה מטריצת הקשיחות. ל- אנו קוראים וקטור העומס (load vector). ככל שאנו מגדילים את אנו מגדילים את סיבוכיות הפתרון, אבל גם את דיוק השיעור שלנו. למעשה, גדול הוא לא בעייתי כל כך בבעיות חד-ממדיות, אבל הוא יהווה מכשול די משמעותי כאשר נידון בבעיות תלת-ממדיות.

הערה:

בתרגול אנו מפרקים את וקטור העומס למטריצת המסה (mass matrix) ווקטור הכוח (force vector). אנו עושים זאת על קירוב בעזרת אותם הפונקציות בסיס:

לכן במקרה שבו תנאי השפה נוימן הוא הומוגני, אנו מקבלים:

ל- אנו קוראים וקטור הכוח (force vector) ולאינטגרל אנו קוראים מטריצת המסה (mass matrix).
פירוק זה עשוי להיות יקר חישובית בבעיות סטטיות, אך בבעיות דינמיות הוא מתבצע באופן טבעי מתוך פיתוח המשוואות.

בניית פונקציות בסיס של FEM

שיטת אלמנטים סופיים מגדירה פונקציות בסיס (למטרת הקירוב) באופן מקוטע (piecewise) לאורך תת-תחום, ה”אלמנטים הסופיים”, של כלל התחום. פונקציות הבסיס הן לרוב פולינומים פשוטים מסדר נמוך. אנו מקפידים שמתקיים:

1. פונקציות הבסיס חלקות מספיק כך שהן חלק מ-.
2. פונקציות הבסיס הן פולינומים פשוטים למקוטעין, המוגדרות אלמנט-אלמנט.
3. פונקציות הבסיס יוצרות בסיס בו כאשר , וגם . בנוסף לכך, לכל ו- מחוץ לאלמנטים עם אותה הצומת .

קבוצה של פונקציות אפשריות לבסיס מוגדרות כ:

כאשר וגם

אחרת, . הנגזרת של הפונקציה היא:

וגם

הפונקציות מסודרות כך שה-”פסגה” של הפונקציה ה--ית מתלכדת עם הצומת ה-. אחת מהיתרונות של בחירת פונקציות זו הוא שהיא מאוד נוחה לאינטגרציה.

בסיס לאלמנט סופי חד-ממדי. למעלה, דוגמה לרשת אחידה. למטה, דוגמה לרשת לא אחידה. (Zohdi, 2018).

הערה:

1. התחום ל- הוא מה שאנו קוראים לו אלמנט סופי (finite element).
2. אם אנו מחלקים את התחום ל- אלמנטים, ולכל אלמנט אנו מבצעים קירוב לינארי (linear approximation), אז יש לנו צמתים (nodes) לאורך התחום.

אינטגרציית גאוס

נחזור טיפה על אינטגרציית גאוס. גאוס שם לב שניתן לבצע אינטגרציה על פולינום מסדר בדיוק לחלוטין עם נקודות דגימה ומשקלים מתאימים. לכן, אנו מגדירים מחדש את הפונקציה על תחום מנורמל כך ש:

כאשר הוא היעקוביאן של הטרנספורמציה. לעומת רוב שיטות האינטגרציה, באינטגרציית גאוס אין את הדרישה שהדגימה על הפונקציה תתרחש בנקודות במרווח שווה ביניהן. כעת, אנו יכולים למשל לבצע אינטגרציה על ביטוי מסדר שלישי (וגם סדר נמוך יותר) עם בדיוק נקודות, כי . לכן

• עבור סדר שלישי :
• עבור סדר שני :
• עבור ביטוי לינארי :
• עבור קבוע (): ישנם כעת ארבעה נעלמים שאנו צריכים למצוא. הפתרון שמקיים את כל הדרישות הוא ו- . במקרה הכללי יותר, ה--ות הן שורשים של משהו שנקרא פולינומי לג’נדר, שזה משהו ידוע, ויש לנו טבלה להן.
  

  אינטגרציית גאוס. (Zohdi, 2018).

טבלה לאינטגרציית גאוס. (Zohdi, 2018).

דוגמה:

נביט באינטגרל הבא:

האינטגרל הוא מהצורה:

עם הטרנספורמציה

כאשר נבצע את טרנספורמציה זו נקבל

כאשר . באינטגרציה בשלוש נקודות, נקבל (הפתרון המדויק הוא ):

טרנספורמציה גלובלית/מקומית

אחד מהיתרונות של שיטת אלמנטים סופיים הוא שהחישוב ניתן לביצוע אלמנט-אלמנט. את מטריצת הקשיחות הגדרנו כ:

ופונקציית העומס:

אנו מחלקים את התחום לאלמנטים ,ואז נפרק את החישובים (אינטגרלים על ) לאלמנטים (אינטגרלים על ), , ( עבור element) כאשר:

וגם:

כאשר וגם .

דוגמה:

כדי לבצע את החישובים באופן שיטתי, אנו רוצים להמיר את הקואורדינטות של כל אלמנט לתחום לטובת אינטגרציית גאוס. לפיכך, אנו צריכים לבצע את הטרנספורמציה הבאה, מהמערכת צירים המקומית למערכת הצירים המרחבית האמיתית,

כאשר הוא המיקום האמיתי במרחב של הצומת ה-, ו- .

טרנספורמציה לאלמנט סופי חד-ממדי. (Zohdi, 2018).

שיטות מיפוי אלו (הטרנספורמציות) לרוב נקראים מפות “פרמטריות”. אם סדר הפולינום של פונקציות הצורה גבוהות כמו קירוב גלרקין על האלמנט, הוא נקרא מפה איזו-פרמטרית. אם הוא נמוך יותר, אז סאב-פרמטרי, ואם הוא גבוה יותר אז הוא נקרא סופר-פרמטרי.

מאפיינים דיפרנציאליים של פונקציות צורה

לקירוב לינארי, פונקציות הצורה המקומיות הן מהצורה:

יש להם את התכונות הבאות:

• לאלמנטים לינאריים יש לנו בסיס הכולל שני צמתים, ולכן שתי דרגות חופש.
• פונקציות הצורה () ניתנות לחישוב באופן יחסית די פשוט. הם יחידה בצומת המתאימה ואפס בכל שאר הצמתים.

הערה:

אנו לא באמת מחשבים את הפונקציות ; אנחנו למעשה מתחילים מ- ואז ממפים (transform) אותם לתחום של הבעיה האמיתית. לכן במטריצת הקשיחות או במטריצת המסה שנראה בהמשך, כל הביטויים חייבים להיות מוגדרים במונחים של קואורדינטות מקומיות (local coordinates).

כעת נגדיר מספר ביטויים בסיסיים שיעזרו לנו בחישוב הערכים המקומיים כמו ו- . למשל, כאשר נבצע את הטרנספורמציה אסור לנו לשכוח שיש גם להכפיל את האינטגרנד ביעקוביאן. במקרה החד-ממדי:

ונשים לב גם ש:

כעת נוכל לפתח ביטוי לנגזרת של (לפי כלל השרשרת):

נציב ב-, ונקבל בעזרת אינטגרציית גאוס:

וגם:

דוגמה:

פיתוח הצורה החלשה:
נתונה הבעיה הבאה:

צורתה החלשה היא (לאחר אינטגרציה בחלקים):

נפרק את התחום לשלושה אלמנטים, ולכן נקבל ארבעה צמתים:

• על אלמנט , מתקיים ו- וגם:

  • על אלמנט , מתקיים ו- וגם:
  • על אלמנט , מתקיים ו- וגם:

  אנו מבצעים את החישובים אלמנט-אלמנט. כל החישובים ב- שייכים לאלמנט , בעוד כל החישובים ב- שייכים לאלמנט , וכל החישובים ב- שייכים לאלמנט .
  בניית המערכת הדיסקרטית:
  עבור אלמנט , כדי לחשב את , נביט בביטוי הבא עבור ערכי :

  ספציפית עבור , הביטוי הוא:

  כאשר אנו מקבלים המון ביטויים (כמו השלישי) שמתאפסים כי פונקציות הבסיס שם הן בתחום האלמנט הראשון. הביטויים כמו מכפילים את הביטויים , המגדירים את המיקום שלהם בתוך מטריצת הקשיחות הגלובלית. אם נחזור על הפעולה עבור עם , נקבל את הערכים הבאים בתוך מטריצת הקשיחות הגלובלית (שהיא מטריצה מסדר כי יש צמתים):

  כלומר, מטריצת הקשיחות המקומית היא רק:

  באותו אופן, וקטור העומס המקומי של אלמנט הוא

  כך שבאופן וקטורי:

  נחזור על תהליך זה על כל שלושת האלמנטים ונקבל את מטריצת הקשיחות ווקטור העומס הגלובליים:

  נשים לב שבאלמנט האחרון יש לנו גם תנאי שפה שאנו חייבים להתייחס אליו בוקטור העומס:

  כמובן שאת כל האינטגרלים נחשב בעזרת אינטגרציית גאוס - כלומר, נשתמש ב- ו-.

אלמנטים ממעלה גבוהה

בהמון מקרים, אם אנו יודעים שהפתרון האמיתי הוא חלק, נעדיף לבצע קירוב ממעלה גבוה יותר. בכללי, אם הפתרון המדויק לבעיה הוא חלק, עבור רשת סבוכה מספיק ( גדול מספיק), הפתרון ממעלה גבוהה יותר מדויק יותר. באותו אופן, אם הפתרון האמיתי לא חלק, עבור רשת סבוכה מספיק, הפתרון עם הפונקציות צורה הלינאריות יותר מדויק.

שלושה אלמנטים ממעלה גבוהה עם צמתים.

כדי לבנות קירוב ממעלה שנייה לכל אלמנט, אנו נצטרך כעת צמתים בכל אחד מהם במקום . כמו מקודם, פונקציית הבסיס חייבת להיות יחידה בכל צומת שהיא שייכת אליו, ואפס בכל השאר. לכן, לאלמנט ממעלה שנייה כללי:

• לצומת : מה שמניב:
• לצומת : מה שמניב:
• לצומת : מה שמניב:

הקשר בין ו- יהיה:

בנוסף, למערכת המשוואות שנבנית מפונקציות אלו עדיין יש את הצורה הלינארית:

כאשר הוא מספר הצמתים בחד-ממד. נביט בדוגמה עם אלמנטים, מה שנותן צמתים. נפרק את האינטגרל לאלמנטים:

עבור אלמנט , עם , אנו צריכים לחשב:

מה שנותן:

בצד ימין של , עבור , אנו צריכים לחשב

ולכן

נחזור זאת על ונקבל:

נחזור זאת על אלמנטים ו- ונקבל:

תרגילים

שאלה 1

הביטו במוט חד-ממדי תחת הכוח המפורש (הצירי) עם תנאי שפה דיריכלה הומוגניים על השפה () ותנאי שפה נוימן הומוגניים על השפה (). הניחו ש- ו- .

סכמת הבעיה.

סעיף א’

רשמו את המשוואות השולטות בבעיה.

פתרון:
לפי , או לפי משוואות תנועה של מוט במתיחה:

סעיף ב’

הגדירו את השאריות על התחום ועל תנאי השפה נוימן.

פתרון:
נסמן את הפתרון המשוער ב- . כמו בשאלה מתרגול קודם:

בנוימן:

סעיף ג’

מצאו את הצורה החלשה של הבעיה.

פתרון:
לפי שיטת גלרקין על המשוואה ועל התנאי שפה:

נרשום את הצורה החלשה (אינטגרציה בחלקים):

נפתח את הביטוי השני כדי לראות שיש ביטויים שמתבטלים:

התנאי שפה ב- מבטלים אחד את השני.
תנאי שפה נוימן הומוגניים תמיד מבטלים אחד את השני בצורה החלשה. לכן קוראים להם לפעמים תנאי שפה טבעיים, והם מתקיימים באופן טבעי. לפיכך, אם למשל אנו לא דורשים תנאי שפה דיריכלה מסוימים, אז זה כאילו אנחנו מניחים תנאי שפה נוימן.
תנאי השפה ב- מתאפס בגלל פונקציות הבסיס שנשתמש בהם. אנו נבחר בפונקציות בסיס כך שכל אחד מהם, חוץ מהראשון , נעלם ב- . הערך של ב- הוא , אבל כאשר נציב תנאי שפה דיריכלה למטריצה , כפי שראינו בשאלה מתרגול קודם, השורה המתאימה (הראשונה) של המטריצה לא תילקח בחשבון.
נישאר עם:

נשים לב שבצד שמאל מופיע מטריצת הקשיחות ובצד ימין וקטור העומס, כמו במשוואה .

נרצה כעת לדון בפונקציות הקירוב/בסיס. כאן , אבל בכללי אנו רושמים:

כאשר הם נקודות הצומת. כעת, בביטוי בצד ימין של נקבל:

האינטגרל שקיבלנו נקרא מטריצת המסה. נציב בחזרה ב-:

כאשר במקרה הזה הוא פשוט יחידה. ניתן לרשום את ביטוי זה כ:

כאשר היא מטריצת הקשיחות, הוא וקטור ההזזות, הוא וקטור הכוח, ו- היא מטריצת המסה. נשים לב גם ש- היא מטריצה סימטרית ומוגדרת חיובית.

סעיף ד’

בחרו פונקציות בסיס לשיטת אלמנטים סופיים.

פתרון:
נבחר בפונקציות בסיס שהן רציפות למקוטעין כפי שמוסבר בבניית פונקציות בסיס של FEM:

פונקציות הבסיס שאנו בוחרים בשאלה.

למשל עבור פונקציות הבסיס יהיו:

זהו פיתוח גלובלי, כדי לפתח באופן כללי נבחן עבור אלמנט כללי - כלומר באופן מקומי.

סעיף ה’

פתחו את הטרנספורמציה/מיפוי למערכת קואורדינטות מקומיות.

פתרון:
בעצם מה שאנחנו עושים זה לקחת כל אלמנט ולבצע עליו טרנספורמציה:

טרנספורמציה מקומית של אלמנט.

במקרה שלנו:

נרצה למצוא את , כלומר, את עבור אלמנט מסוים. לפי משוואה :

כאשר .
לחישובים בהמשך, נצטרך לחשב עוד ערכים: היעקוביאן:

כאשר אנו קוראים את כ-”היעקוביאן של הטרנספורמציה של אלמנט ”.

נסמן את הנגזרת לפי עם , ונראה כי:

הנגזרת לפי :

סעיף ו’

מצאו את הקשיחות ה-מקומית ומטריצת המסה ה-מקומית של כל אלמנט.

פתרון:
לפי משוואות

ולכן:

מטריצת המסה:

לכן:

נבחר לחלק ל- אלמנטים. זה אומר ש-. לכן:

על כל אחד מהם אנו יכולים לבצע אינטגרציית גאוס עם נקודה אחת (האינטגרנד לינארי) עם ו-. נקבל:

לכן המטריצות המקומיות (נציב ):

סעיף ז’

פתחו ביטוי לקשיחות ה-גלובלית ומטריצת המסה ה-גלובלית.

פתרון:
חילקנו ל- אלמנטים, ולכן יש לנו פונקציות בסיס כפי שהצגנו ב[[#שאלה 1#סעיף ד’|סעיף ד’]]. לפי הדוגמה, כאשר נשלב את המטריצות המקומיות לגלובליות נקבל (למטריצת המסה אנו מבצעים את אותם הפעולות כמו למטריצת הקשיחות):

סעיף ח’

מצאו את הקירוב לפתרון בשיטת אלמנטים סופיים.

פתרון
מבחינת וקטור הכוח, הוא פשוט . מערכת המשוואות היא:

נתחיל לחשב:

כאן אנחנו כבר צריכים להיזהר ולזכור שיש לנו תנאי שפה מסוג דיריכלה ב- . ניזכר שבחירת ה--ים בקצוות לא שרירותית, כך שהמשוואות שפיתחנו עבורן לא תקיפות, ובכלל, לפי איך שאנו מקרבים את הפתרון , צריך להתקיים:

לפיכך, נישאר עם מערכת המשוואות המצומצמת:

נקבל ש-:

כעת אנו יודעים את הערכים של הפתרון בצמתים, ונוכל לבצע אינטרפולציה כדי לקבל את הפתרון המקורב. נשים לב שלפי [[#שאלה 1#סעיף ד’|סעיף ד’]]:

ולכן הקירוב:

נציב ערכים ונקבל:

פתרון מדויק ופתרונות בשיטת אלמנטים סופיים עם מספר אלמנטים שונים.