הדף הזה יכיל קישורים לכל הדרוש על מנת לעשות סימולציה נומרית ב-OpenFOAM של הזרימה מעל אזור קרקעי כלשהוא בכדוה”א.

קישור לדף עם נתונים מלאים של געת Bolund מחוץ ל-Risoe בדנמרק –
http://windenergyresearch.org/2010/09/the-bolund-project/

שלב ראשון – נתוני טופוגרפיה (z(x,y

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

• באתר earthexplorer.usgs.gov אפשר לחפש נקודה רלוונטית.
• כשמוצאים אותה, להכנס לכאן, לנתוני ה-SRTM3 ולהוריד את הקובץ המתאים. הנתונים הם בדיוק של 3 שניות קשת מחוץ לארה”ב ו-1 שניית קשת בארה”ב. זה מתרגם לבערך 90 מטר מחוץ לארה”ב ו-30 בתוך ארה”ב. זה לא מידע מספק לאזורים ממש מורכבים, אבל לאזורים רבים, גם הרריים, זה מספק. לאזורים בהם זה לא מספיק – צריך להשיג נתונים משרות מיפוי מקומי (או של המדינה הרלוונטית, כמו מפי בישראל)

עבור ה-Askevien hill (בו אשתמש להדגמת הטכניקה) המשבצת היא N57E008.hgt.zip

שלב שני – העלאה ל-SAGA

אני משתמש בתוכנת SAGA הפתוחה בשביל עיבוד המשטח והעברתו מקואורדינטות קווי אורך ורוחב להיטל על מישור, ותוכנת google earth החינמית בשביל לדעת היכן המיקום המתאים לחיתוך.

• לאחר unzip של הקובץ, העלאת קובץ hgt מתוך SRTM3 ל-SAGA:
  Modules–>File–>Grid–>Import–>Import USGS SRTM Grid
  בוחרים את הקובץ המתאים ולוחצים ok (לשים לב שיש תפריט דומה שמתחיל מ-files ולא מ-moduls. לא להתבלבל)

• מגדירים את המערכת קואורדינטות של הרשת ע”י כניסה ל-Moduls–>Projection–>Set Coordinate Reference System
• תחת EPSGcode
  • : מכניסים 4326 עבור  ((Longitude / Latitude (WGS 84).
  • שאר הפרמטרים נשארים זהים
  • ב-Data Objects מוסיפים לשורת ה-grid את ה-grid שלנו

• מעבר מ-Lat/Long ל-UTM
  צריך ל”שטח” את הכדור בקטע הרלוונטי, כלומר לעשות projection. דרך אחת מתאימה תהיה לעבור ל-UTM. בישראל משתמשים ב-ITM שהוא סט פרמטרים טיפה שונה, אך אותו סוג היטל.
  • בעזרת גוגל ארץ נמצא את האזור של ה-UTM.
    נכנסים ל-options בתוך גוגל ארץ (במק – דרך preferences, בחלונות tools–>options–>UTM). שם תחת האזור שכותרתו “Show Lat/Long” בוחרים את ה-UTM, הלא הוא ה-Universal Transverse Mercador
  • עוברים עם הסמן על הנקודה במפה אותה רוצים למדל. בתחתית המפה יהיה כתוב למשל 29 V 598060 m E 6339670 N elev 112 m
  • המשמעות היא שזה אזור V 29
  • ב-SAGA נכנסים ל-Modules–>Projection–>coordinate transformation grid
  • EPSG code: הפעם בוחרים את זה שמתאים לUTM zone שלנו – במקרה דוגמא, זה 32629
  • בתוך Data Objects, תחת grids, בחירת ה-grid system וה-source.
    ** יש בעיה של זיהוי ה-epsg משום מה ב-saga 2.08. הפתרון הוא להכניס אותם בתפריט ה-30 parameters באופן מפורש. כלומר לבחור תחת proj4 parameter תחת dialog, לבחור ב-projection type את ה- Universal transverse mercador וב-projection settings לבחור את ה-zone הדרוש
  • בשלב זה יופיע תפריט שיבקש את גבולות הגזרה. אז כבר אפשר לחתוך את הרשת לגודל הדרוש לסימולציה. אבל עדיף להשאיר את זה בגבולות המקוריים, ולחתוך בשלב האינטרפולציה.

• בנקודה זו אפשר לעשות את ההחלקה (כפי שיוסבר עוד מעט) ב-SAGA. לחלופין, אפשר להשתמש בתוכנות אחרות כמו pointwise או gridgen למשל (אשמח לגרסא פתוחה של אלו, אולי אצליח לעשות זאת עם salome בעתיד)
• הערה – כדאי לתת לדברים שמות משמעותיים. הבעיה – הדרך היחידה לעשות את זה היא לשמור את ה-grid. כפתור ימני על מה שרוצים לשנות לו שם, ושמירה. השם לא מתעדכן על המסך – הוא יתעדכן אם תשמרו כ-project ותעלו. אבל הוא משתנה ובכל התפריטים הבאים יופיע השם שנתתם.

שלב שלישי – החלקה (עם SAGA)

• אינטרפולציה וגזירה במהלך אחד, מתוך מפת ה-UTM הגדולה
  • Modules–>Grid–>Gridding–>Spline Interpolation–>Multilevel B spline Interpolation (from grid)
  • Maximum level = 14
  • לסמן את ה-grid
  • לסמן את ה-update view
  • עכשיו יתקבל חלון להכנסת הגבולות הרצויים – אז נאמר, פעם אחת מלבן גדול, 30 ק”מ, ופעם אחת קטן, 5 ק”מ. חישוב הקואורדינטות פה.
  • קודם לשנות את ה-cellsize. למשל – ל-50 מטר.
  • אח”כ להכניס את הגבולות Left    583000
    Right    613000
    Bottom    6325000
    Top    6355000

• Modules–>Grid–>Construction–>create constant grid
  • לבחור 5 מטרים (עבור הדוגמא) ואת הגריד האחרון שנוצר (the bspline)
• ייצור 2 מעגלים ונקודה מרכזית
  • לייצר קובץ טקסט עם הנקודה המרכזית – פשוט שורה אחת של 598060  6339670  5
  • Modules–>File–>Shapes–>Import–>import point cloud from text file
  • להצביע לקובץ הטקסט
  • להגדיר את הקואורדינטות בהתאמה, ו-seperator = space
  • Modules–>Shapes–>Tools–>Shape buffer
  • לבחור 2000 מטר (או כל מרחק אחר) ולבחור את ה-shape שמכיל את הנקודה הבודדת
  • לעשות שוב עבור 4000 מטר
  • לשים לב לייצר buffer חדש ולא לדרוס את הקודם

• בחירת פונקציית אינטרפולציה
  • למשל 1-(R-Rin)/(Rout -Rin) (צריך להציג בצורה קרטזית, כפי שיופיע למטה עבור הדוגמא הזו)
  • Modules–>Grid–>Calculus–>Grid generation–>function
  • הכנסת הגבולות של הרשת
  • (sqrt((x-598060)^2+(y-6339670)^2)-2000)/(4000-2000)
  • grid–>calculus–>grid calculator
  • בחירת הרשתות בסדר הבא: רשת מקורית (bispline), רשת constant grid, ולבסוף interpolation function
  • הכנסת הפונקציה הבאה: g1*(1-g3)+g2*g3 (נכון עבור SAGA 2.0.8)

• תפירת הכל ביחד
  • Shapes–>grid–>spatial extent–>clip grid with polygon
  • בחירת הגריד המקורי
  • בחירת ה-buffer של ה-2 ק”מ (הפנימי) כפוליגון
  • ביצוע
  • לחזור על אותו גבר עם הגריד שעבר חישוב אינטרפולציה (השם שלו הוא calculation… אם לא שיניתם אותו) וקוטר של 4 ק”מ (או הקוטר החיצוני שנבחר)

• Grid–>Construction–>merging
• בוחרים את המעגל הפנימי, המעגל החיצוני, וה-constant grid
• להקליד ב-overlapping cells את first value in order of grid list ולדאוג שהסדר ברשימה יהיה מהמעגל הפנימי (הרשת המקורית) לאינטרפולציה, והרשת הקבועה
• יש רשת!

שלב רביעי – המרה ל-STL

• Modules–>TIN–>conversions–>Grid to TIN
• לוקח כמה שניות ארוכות
• Modules–>file–>shapes–>export–>export TIN to stereo litography file (STL)

שלב חמישי – יצירת רשת חישובית

ניתן לעשות את זה עם הכלי snappyHexMesh שמגיע יחד עם OpenFOAM. זה כלי שיוצר רשת לא סדורה (unstructurred) עם ציפוף באזורים מוגדרים, ואפשרות הכנסת שכבות תאים על פני משטחים (boundary layer cells).

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

שימוש ב-Pointwise:

ייבוא הגאומטריה –
import–>database , ובחירת קובץ ה-STL. הרשת המשולשת שנוצרה בשלבים הקודמים תראה ככה במבט Y+ למשל:

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

• קודם כל, נכנס ל-deafaults בתפריט השמאלי (שורת התפריט מכילה משמאל לימין list layers deafaults) ונבחר מספר כלשהוא ל-dimension. למשל 64. זה המספר של הקטעים שינתן לכל מחבר (connector) חדש שניצור. אח”כ אפשר כמובן לשנות ת המספר ולשנות את פונקציית הריווח של הנקודות לאורך המחבר.
• הקואורדינטה של מרכז הרשת היא 598167.367    6339602.511, בגובה 0 (אח”כ אנחנו נעשה project על ה-database שמתאר את הגבעות שלנו, אז הבחירה של הגובה כרגע היא די שרירותית).
• create–>2 point curves
• הכנסת הקואורדינטה המרכזית. לחיצת enter.
• סימון משבצת ה-advanced. מכיוון שכל הקואורדינטות שלנו הן מספרים לא עגולים (במקרה הזה, אבל לא דוקא תמיד. כשיש חופש היכן לבחור את המרכז – למשל, כשעושים סימולציה עבור מציאת האזורים הטובים ביותר בשטח מסוים, אפשר לעגל ואז זה פחות בעייתי) אנו נעבוד עם xyz offset. בשלב ראשון ניצור את המחבר האנכי, ונכניס 0 0 6500 (הדוגמא של הרשת מגיעה מהפרסום של Martinez 2011)

 

 

שלב ששי – הגדרת תנאי שפה והתחלה ופתרון נומרי של משוואות הזרימה

הפתרון נעשה פה ע”י simpleFoam , עבור תנאים ניטראליים (בהם הטמפרטורה הפוטנציאלית קבועה עם העליה בגובה מפני הקרקע)

שלב שביעי- דגימת התוצאות לאורך קווים ומשטחים שונים

עבור הדוגמא שמובאת בדף הזה, ה-Askervein hill, נעשו מדידות לאורך שני קווים, ובמספר נקודות לאורך מגדלי המדידה. בשביל לדגום את התוצאות נשתמש בפקודת sample. עבורה, יש להגדיר את הקובץ system/sampleDict . בשביל זה אנו צריכים לדעת מה הגובה של הקרקע בנקודות אותם אנו רוצים לדגום. בשביל זה – חוזרים ל-SAGA. (הערה – על מנת לעשות את זה באופן עצמאי, אני אכתוב סביר להניח קוד ייעודי לעשות את כל השלבים הרשומים למעלה עם SAGAS ורשומים פה. אולי עם הספריות של SAGA, או עם הספריות המקוריות ש-SAGA עושה בהם שימוש).

• נגדיר קובץ טקסט של הקווים הרצויים עם הקואורדינטות ב-UTM – לפי איזה דיסקרטיזציה שנרצה
• Modules–>File–>Shapes–>Import–>import point cloud from text file
• להצביע לקובץ הטקסט
• להגדיר את הקואורדינטות בהתאמה, ו- seperator = space או מה שבו השתמשתם
• שלב הבא – הוספת ערכי ה-grid.
  • shapes–>grid–>grid values–>add grid values to points
  • בחירת הנקודות שנוצרו מהעלאת קובץ הטקסט הקודם, ואיזה אינטרפולציה שרוצים
• עכשיו שהנקודות הם shape – אפשר לשמור את התוצאה כ-xyz
  • Modules–>files–>shapes–>export–>export shapes to XYZ
  • תחת attributes לבחור את השם החדש שנוצר (כלומר – לא X,Y או z) – אצלי למשל זה “merged grid”, שזה השם של הרשת האחרונה שנוצרה בשלב 3)
  • הקובץ שנוצר אמור להכיל 3 עמודות – X, Y , והאחרונה היא למעשה z
  • לפעמים מתווספת שורה ראשונה עם הערך ” 0 0 -99999″ – פשוט צריך למחוק אותה

עכשיו אפשר להכניס את הנקודות שנוצרו