Final Report
Microgravity Fluid Mechanics and Heat Transfer Computation
( http://www.cfdlab.ae.utexas.edu/nasa_hpcc/microG/ )
Gary Cox and Oliver Pozo
Project Advisors: Dr. G.F. Carey and Dr. R.O. Stearman
Project Assistants: Bill Barth, Alexandre Ardelea
Facilities: CFDLab, Dept. Aerospace Engineering, UT Austin
ASE 363Q – Spring 1999
Department of Aerospace Engineering
The University of Texas at Austin
May 10, 1999
Memorandum
Date: May 10, 1999
TO: Dr. Ronald O. Stearman
FROM: MicroG Research Group
RE: Microgravity Fluid Mechanics and Heat Transfer Computation Midterm Progress Report
Dear Sir,
The MicroG Research Group is pleased to present our final report on the Microgravity Fluid Mechanics and Heat Transfer Computation Project. Our team was tasked with simulating fluid motion and heat transfer under microgravity conditions using specially-designed software on a parallel-computing cluster. In addition, we began preliminary research into designing boundary conditions and geometries which will create a specified flow and heat transfer. Included in this report is a brief summary of the project background, a detailed discussion of our individual member’s tasks and the theories underlying our methodology. This report is concluded by a description of the work completed.
We would like to announce that our team made satisfactory progress towards our goals. We successfully compiled, configured, and ran the code which computes the fluid behavior and carried out supporting simulation studies for our analysis and parametric studies. In addition, we have been able to visualize the results using various graphical software packages. We were successful in our goal of designing boundary conditions and fluid geometries in order to induce a specified flow and heat transfer for simple 2-D cases.
While much has been accomplished, there is much room for future work to be carried out in continuation with this project. We hope that our work will serve as an aid to other teams performing research in this area in the future.
Sincerely,
Gary Cox
Oliver Pozo
Abstract
In a microgravity environment, such as the International Space Station, where the effects of gravity are greatly reduced, scientists can conduct research that achieves results not possible in ground-based laboratories. Computer simulation can help understand the basic fluid physics, as well as help design the experiments or systems for the microgravity situation.
In our research, we are specifically targeting the use of computation in order to simulate coupled fluid flow and heat transfer in a microgravity environment. Through the use of parallel computing facilities at the Computational Fluid Dynamics lab of the Aerospace Engineering Department, UT Austin, we are able to run multiple scenarios in order to find those flows that are the most interesting for physical experimentation. Computer simulations of idealized cases allow us to better understand the physical processes behind many of the observed physical phenomena in a variety of microgravity environments and non-linear flows. Sensitivity studies to various parameters can be carried out and preliminary design studies of fluid thermal systems were determined to be feasible through designing optimization problems on simple heat-transfer cases.
Acknowledgements
The MicroG group would like to thank the following individuals and organizations for their assistance and sponsorship of this project:
Dr. C.F. Carey – Project Sponsor/Coordinator
Dr. R.O. Stearman – Team Sponsor
William Barth – Graduate Advisor/Technical Assistant
Alexandre Ardelea – Post-Graduate Advisor/Technical Assistant
NASA – Sponsor of ESS Grand Challenge Project NCCSS-154
DARPA – Sponsor for grant DABT63-96-C-0061
Table of Contents
1.0 Introduction
*1.1 Description of the Project and Objectives
*1.2 Experimental Research
*1.3 Member’s Roles
*2.0 Background Theory
*2.1 Introduction
*2.2 Mathematical Formulation
*2.3 Parallel Implementation
*2.4 MatLab as a Design Tool
*2.4.1 Background on the Toolkits
*2.4.1.1 PDE Toolkit
*2.4.1.2 OPT toolkit
*3.0 Technical Discussion & Results
*3.1 Parametric Studies
*3.1.1 Description of Parameters and Domain
*3.1.1.1 Fluids and Initial Conditions
*3.1.1.2 Domain and boundary conditions
*3.1.2 Description of Test Cases and Results
*3.1.2.1 Design Case 1
*3.1.2.2 Design Case 2
*3.2 Preliminary Design and Optimization with MatLab
*3.2.1 Background on the Toolkits
*3.2.1.1 PDE Toolkit
*3.2.1.2 OPT toolkit
*3.2.2 Problem Setup
*4.0 Conclusions
*References
*APPENDICES
*List of Figures and Tables
Figure 1 - Drop Tower Canister being decelerated in the deceleration pit *
Figure 2 - The "Vomit Comet" *
Figure 3 - Fluids Combustion Facility -- International Space Station *
Figure 4: Beowulf cluster network topology *
Figure 5: Description of Domain *
Figure 6: Linear Temperature Distribution *
Figure 7: Circular Temperature Distribution *
Figure 8: Design Case 1 Boundary Conditions *
Figure 9: Design Case 1 Roll Formation *
Figure 10: Design Case 2 Boundary Conditions *
Figure 11: Design Case 2 First Time Step *
Figure 12: Design Case 2 Time Step 25 *
Figure 13: Design Case 2 Time Step 300 *
Figure 14: Design Case 2 Time Step 1000 *
Figure 15: Design Case 2 Fictitious Fluid Time Step 2 *
Figure 16: Domain Layout with Mesh *
Figure 17: Initial guess for optimizer and the resulting solution *
Figure 18: Optimization in progress *
Figure 19: More Progress for the optimizer *
Figure 20: Final solution for optimizer *
Figure 21: Current Process Flowchart *
Table 1 -- Project Schedule *
The Microgravity Fluid Mechanics and Heat Transfer Computation senior analysis and design project is related to the UT research work supported by NASA on "Scalable Parallel Finite Element Computations of Rayleigh-Benard-Marangoni Problems in a Microgravity Environment". Dr. G.F. Carey, Director of the Computational Fluid Dynamics (CFD) Lab in the Aerospace Engineering Department at the University of Texas at Austin is supervising the project. There are several other people involved in the larger project that have been of assistance to us during the start-up process. NASA’s grant was awarded to the CFDLab in 1996 and since then they have been working on research related to microgravity applications. This research has made it possible for the project members to develop a finite element code that can be run in parallel processing mode. This code, MGFLO, developed at the CFDLab uses finite elements and an iterative process to compute solutions to fluid mechanics and heat transfer equations to solve problems in, for example, a microgravity environment. The code allows the user to vary gravity, study different fluids, investigate fundamental fluid mechanics, and change the initial conditions, as well as the boundary conditions. Another feature of the code is the ability to solve the above mentioned equations in either a coupled or uncoupled manner. More information about the formulas and the way the code runs will be provided in the Background Theory section, Section 2.0 in this report.
1.1 Description of the Project and Objectives
This project will use the MGFLO code to analyze and design fluid/thermal systems in a microgravity environment. This design process entails developing an understanding of the fluid/thermal applications involved in the systems. The project will be complemented with a feasibility study of the implementation to the MGFLO code of an optimization loop for designing fluid/thermal systems, which will be done with MATLAB.
The objective of the project is to carry out application tests and parametric studies for heated fluids in a microgravity environment. We also plan to investigate control of the flow using applied thermal input. There will also be some related MATLAB finite element/optimization work.
To obtain background information we began by investigating sources in the World Wide Web (Web) and the UT Engineering Library for literature research. We read several technical papers, published by CFDLab members and other researchers. Many of these papers and information refer to work NASA is doing in this area. NASA Glenn’s (formerly known as NASA Lewis) Web site provided a valuable source of information on current research activities and on earth-based microgravity experiments. (Actually, the Marshall Space Flight Center is the lead center for NASA’s Microgravity Research, but the microgravity research in fluid physics is implemented at Glenn.)
In addition to their earth-based microgravity physical experiments, NASA Glenn also houses the Computational Microgravity Laboratory (CML) where tens of experiments similar to the ones we are conducting are being carried out. Some of the more interesting experiments included Bubble Dynamics on a Heated Surface and Protein Crystal Growth. Following the literature study, we began work on applications tests. This process involved compiling the code and running test case scenarios with idealized fluids. More information about the project and the results will be discussed in Section 3.0 later in this report.
NASA is the main group that conducts experimentation in microgravity in the US. They have several research facilities that can allow physical experiments to be exposed to a microgravity environment (ranging between 2.2 seconds and 14 days of exposure). In fact in the near future (as early as October and at least by late January 2000) they will be able to house experiments, in a microgravity environment, for extended periods of time by using the International Space Station.
The most readily available facility for microgravity testing is the drop towers that are located at NASA Glenn. These drop towers can expose experiments to a microgravity environment for times ranging between 2.2 to 5.2 seconds. In the figures (Figure 1) below one can see the deceleration of a test canister in the 2.2-second drop tower.
Figure 1 - Drop Tower Canister being decelerated in the deceleration pit
Another research facility that NASA uses for longer periods of microgravity exposure is the "Vomit Comet". This KC-135 aircraft is housed at NASA Johnson in Houston and is used for experiments that need about 30 seconds of exposure to microgravity. The pilot, by flying the plane in a parabolic path accomplishes to recreate a near microgravity environment inside the plane. Each time the plane takes off the ground, it flies about 40 parabolic curves of 30 seconds each.
When the experiments being conducted need a semi-extended period of microgravity exposure, NASA puts them in the Space Shuttle inside their USML (United States Microgravity Lab), located in the payload bay. There is a great amount of scientists that want to get their experiments into the Space Shuttle’s USML schedule, but only a few experiments can be selected, due to space constraints, each time the Space Shuttle flies. In the near future though, when the ISS is habitable and the FCF (Fluids and Combustion Facility) is installed in the US Laboratory module in the Space Station, more complicated and lengthy experiments will be able to be conducted in microgravity.
The FCF is composed of three main racks. The first rack is the Combustion Integrated Rack (CIS), and will accommodate combustion science principle investigations. The second rack is the fluids integrated rack (FIR), scheduled to launch in 2002 and will house fluid physics principle investigations. The third rack, still being designed, will be launched in 2003, completing the Fluids Combustion Facility on the ISS.
Other facilities that are available to NASA, as mentioned in the previous section of this report, are the CML which deals with computer based microgravity research. Also many Universities in the US get involved with NASA sponsored programs. For example the University of Texas at Austin is involved in many areas with NASA projects. This microgravity project is part of NASA Grand Challenge project given out to nine Universities in the US. The Pharmacy Department at the University of Texas is also involved with NASA and has flown experiments on both the "Vomit Comet" and the Space Shuttle.
Figure 3 - Fluids Combustion Facility (FCF) -- International Space Station (ISS)
The effective control of fluid/thermal systems in microgravity in the near future will allow the growth of extremely pure crystals, as well as drugs that might help cancer disease.
The senior project team is composed of two members: Gary Cox and Oliver Pozo. We will collaborate with Dr. Carey and graduate and post-graduate researchers in this project.
Gary Cox:
Literature Research
Assistant in running the code
Assistant in visualization of results
Programmer for preliminary design phase (MATLAB)
Oliver Pozo:
Compiling and Running the Code
Visualization of Results
Assistant in Literature Research
Assistant in programming for the preliminary design phase
Background Information Research
WEB design
A quick description of each of the roles is specified in the following paragraphs.
Literature research was conducted at the UT Engineering Library as well as on the WEB. Interesting papers and WEB sites that were found and used in this report can be found listed in the Reference section at the end of the report.
The MGFLO code was compiled with the initial conditions pertinent to our project. With that done, the code was ready to run in the parallel cluster at the CFDLab. Then several preliminary and informative tests were run where the boundary conditions were varied. The different boundary conditions include temperature conditions, slip/no-slip conditions as well as flux/no-flux conditions. Some of the parametric constants that can also be modified include the gravity vector, the surface tension and the viscosity. For example, surface tension was varied from –0.01 to-10.0. For these selections, the code was run in parallel on the Beowulf cluster (described in Section 2.3 in this report). Results were visualized and put into a movie format with TecPlot. These test cases are described in depth in Section 3.0 in this report.
Background information that was gathered prior to the more detailed literature research consisted of several meetings with Bill Barth and Alexandre Ardelea. During these meetings we visited the Non-Linear Dynamics Lab in the RLM building at the University of Texas at Austin. During this visit we saw several physical experiments that were being conducted on Rayleigh-Benard-Marangoni effects in thin layer fluids. Bill Barth gave us a presentation on the CFDLab’s computer cluster and on the MGFLO code that we were going to use during most of our design project. This preliminary stage was important in order to know what kind of technical papers to look for and where to start looking in the library and on the WEB.
With the rapid expansion of the Internet for technical transfer and discussion of experimental results, it is important that all presentations and technical reports be accessible from the Internet. This important feature is referred to as WEB design in the previous list that describes each member’s roles. All the memos and progress reports are available in HTML format as well as the original format they were written in. Links to interesting sites on the internet as well as a comprehensive bibliography can be found in the MigroG site, which can be found at http://www.cfdlab.ae.utexas.edu/nasa_hpcc/microG/ .
An understanding of the behavior of heated fluids in microgravity conditions is important for designing useful experiments for the space shuttle and the international space station. In addition, such understanding is important for the future design of thermo-fluid systems and machinery that might be employed in microgravity environments. There is a wealth of knowledge of terrestrial behavior heated of fluids, with experimental data that spans generations. However, very little is known of the microgravity behavior of fluids due to the relative difficulty of obtaining experimental data under such conditions.
Some examples of coupled heat and fluid flow of which we are interested includes hot-forming manufacturing, crystal growth, chemical flow processes for pharmaceuticals, volcanic lava flows, and natural and forced convection systems. In a terrestrial setting, these flows are characterized by natural convection. In natural convection, buoyancy is the dominant factor driving the flow. Benard performed a series of experiments in the early 1900’s on heated fluid layers. These experiments later became known as the Rayleigh-Benard problem after Rayleigh provided a mathematical analysis of buoyancy-driven flow.
The role of thermocapillary surface tension in influencing fluid flow wasn’t appreciated until later. Usually surface tension is associated with the curvature of the surface as in the minimal surface problem for soap films, droplets, and interfaces between fluids. Because the surface tension coefficient can vary significantly with temperature, variations of temperature on the surface will cause surface tension gradients which will then produce shear stresses that move the fluid. This effect, called the Marangoni effect, is now known to be a much more significant contributor than buoyancy to the effects first observed by Benard in the thin fluid layers, so now the problem is called the Rayleigh-Benard-Marangoni (R-B-M) problem. [1]
In microgravity conditions, buoyancy effects are negligible. Therefore the thermocapillary surface tension effects become the primary driving force behind the fluid flow. Control of the fluid behavior can be exercised by the application of direct or indirect surface heating. This heating causes controlled changes in the surface tension properties of the fluid, and therefore the shear forces inside the fluid which will direct its movement. In addition, thermocapillary shear flow plays a key role in the propagation of flame and other combustion fronts in microgravity conditions.
Terrestrial experiments in thermocapillary effects have been limited to the study of thin fluid layers. In thin layers, buoyancy forces become less important with respect to the surface tension effects. While much valuable research has been done using such setups, it is obvious that restricting experiments to thin-layer geometry limits the applicablity to general geometry cases. Mathematical modeling and computer simulations, on the other hand, are not restricted to such thin layers, so any arbitrary geometry can be simulated in addition to thin layers.
In the following section we follow precisely the formulation in reference [1].
In describing how we go about modeling fluid flows, it is necessary to briefly develop the mathematical framework involved in thermo-fluid flow. Consider 3-dimensional transient fluid flow and heat transfer for a viscous incompressible fluid with a free surface. The mathematical model for such flow is described by the Navier-Stokes equations with the Boussinesq approximation, and the convective heat transfer equation. Due to the variations in surface tension from temperature gradients, there is a shear stress at the free surface. We only consider a fluid geometry with flat surfaces, as a deformable free surface is not currently implemented in the model. The governing viscous incompressible flow equations are as follows:
[1]
[2]
In flow domain W where u is the velocity field, t is the stress tensor and is specified by Stokes hypothesis for a Newtonian fluid, f is an applied body force, g* is the microgravity vector, T* is the reference external temperature, T is the fluid temperature and b is the thermal coefficient. The no-slip condition applies at the walls, so that u equals uwall at the boundaries.
There is a shear stress due to surface tension effects at the free surface. On the horizontal free surface, the tangential shear stress is given by
![]()
[3]
with a similar expression for t zy where g (T) is the surface tension coefficient and T is the temperature.
The relationship describing the heat transfer within the fluid is given by the following:
[4]
where k is the thermal conductivity of the fluid, r is the density, cp is the heat capacity, and Q is a heat source term.
Taking the above equations, integrating by parts, introducing Stokes Theorem, and performing additional mathematical manipulation yields the following result:
[5]
Now, we introduce the discretization of elements and the finite element basis functions. The resulting pressure and temperature expansions are given by the following:
[6]
[7]
[8]
where k is the velocity component index and ujk , pj and Tj are the nodal values. Taking the previous results and converting to matrix form yields the following relations:
[9]
[10]
[11]
With the base relations now written in a form easily dealt with in a numerical manner, a standard q method for numerical integration is used. This method employs forward, midstep, and backward-Euler integrators. Both coupled and decoupled algorithms can be utilized, either solving temperature and velocity fields in parallel or in serial fashion.
In the present work a simple non-overlapping decomposition of the finite mesh is used. This means that the processor interfaces coincide with a subset of element faces and the nodes on those faces are shared by adjacent processors. This in turn implies that the element calculations can be easily parallelized over the processor partitions. Such a method provides an efficient parallel strategy, as it allows for communication requests to be initiated for the subdomain border element calculations while computation begins on the interior subset. With that initiated, the border element computations can be completed. Taken as a whole, the strategy will permit complete communication overlap provided that there are sufficient elements in the interior subset.
When the global mesh is partitioned across processors, the "send list" of nodes is also set up. These are nodes on the subdomain boundaries and therefore are shared by neighboring processors. A local and global numbering system is used for identification of each node. The implicit equivalence of local node values avoids sending global node numbers across the network, which saves valuable bandwidth. However, it also complicates the forming of the send lists. Fortunately, since the send lists only need to be generated once for the static processor partitions, it is not a major concern for the current implementation of the program.
The program is designed for fast scalable parallel computation of steady and transient solutions to coupled incompressible viscous flow with heat and mass transfer. Parallel efficiency is achieved by careful implementation of MPI (Message Passing Interface) and customized communication software. Most of the computation time is spent solving large sparse systems, so much time has been spent in optimizing the solver. Matrix-vector (MATV) products are repeatedly needed throughout the solving process, therefore special consideration was given in this area. The optimized assembly-language solver BLAS MATV product routines are used, which result in roughly doubling the performance of the computations.
The resulting software is called MGFLO, which is designed for unstructured mesh finite-element simulations on irregular partitions. While designed to run on many different parallel-computing platforms, continued work has been done with a Beowulf-type cluster comprised of 16 Intel Pentium-II processors for the investigation of parallel computing on workstation clusters. Each node is equipped with a single processor, 128 MB of RAM, and a dedicated 100 Mb Fast Ethernet network connection. These workstations comprise a relatively-low-cost cluster built from commodity of parts. In various applications, the cluster’s performance has been shown to scale linearly with the number of processors utilized up to the maximum of 16 that are present in the cluster. The investigation of workstation parallel-computing has shown that a low-cost alternative to dedicated supercomputers exists for preliminary studies, and that the investigation of smaller fluid problems via software such as MGFLO on these low-cost clusters is viable.
Figure 4: Beowulf cluster network topology
MatLab was chosen as the programming environment for the optimization design since it has many built-in tools that prove invaluable for rapidly developing an application in controls and optimization. The two toolkits that were principally used were the Partial Differential Equations (PDE) toolkit and the Optimization (OPT) toolkit. Together, these toolkits allow us to create a 2-D or 3-D flux/flow problem (such as heat transfer), and optimize the inputs in order to minimize or maximize a specified output. In the future, we would like to implement such a feedback and optimization system into MGFLO in order to design particular fluid flow patterns as well as heat flow by manipulating the boundary conditions. However, for now we have restricted ourselves to a simpler case in order to validate the concept.
In beginning this phase of the design, it was decided that a problem in heat transfer through a 2-D domain would be easy to rapidly develop, yet would adequately demonstrate the capabilities of the optimization routines as a proof-of-concept for later adaptation and implementation into the larger MGFLO project. Several ideas on how to go about creating a suitable problem were discussed, and in the end a problem dealing with a moving heat source on an edge whose location was to be optimized in order to minimize or maximize the flux leaving another edge was decided upon. In order to arrive at a non-trivial, "interesting" answer, a non-symmetrical geometry for the problem was decided upon.
2.4.1 Background on the Toolkits
Before beginning the actual design of the optimization routines, it was necessary to first research the capabilities of the tools on hand in order to learn what was possible and what wasn’t within the constraints of the programming environment that MatLab provides. The problem relies principally on two toolkits: The PDE and the OPT.
The PDE provides the necessary routines, as well as a graphical user interface (GUI), that allow arbitrary geometries to be set up in various flow/flux problems. Incompressible fluid flow, heat flow, and many other "genres" of problems are available for solving via default equation settings that automatically adjust the governing differential equations and variable names to match those commonly used in each respective discipline. For our design problem, we chose to use the heat transfer setting as that setting provided the correct default environment for specifying the heat transfer-related parameters that we would need for the investigation. The functioning of the PDE toolkit itself is fairly straight-forward. The toolkit takes a specified geometry and creates a finite-element mesh that overlays the geometry, thus dividing the domain into many sub-elements. The fineness of the mesh can be increase by the user in order to get more accurate results at the expense of more processing time. Boundary conditions are then specified by the user at each edge of the domain. In our case, we had the option of specifying a flux condition, a fixed temperature condition, or a combination of the two. After the mesh is created and the boundary conditions set, the steady-state solution of the problem can then be calculated using other routines of the toolbox that solve the governing equations specified by the problem type (heat transfer, fluid flow, etc.). After solving, the user can then visualize the results in a number of ways through color coding and use of vector arrows. In the case of our heat transfer problem, temperature gradients and heat flow could be visualized through the color coding and vector arrows. By manually adjusting the boundary conditions and re-running the solve routine, a parametric study can be conducted on the problem of interest.
The optimization toolkit of MatLab allows a user to take an input and an output function, and let MatLab automatically change the value of the input in search of a maximum or minimum for the output function. This routine works similarly to a root-finding algorithm, in that the input function is generally dependent on a single variable that can be changed, such as a position. The OPT toolkit does not include a specialized GUI, although the results can be visualized through the graphical tools available through the base MatLab application. When used in conjunction with the PDE toolkit, the OPT toolkit allows for a rapid and automated parametric study that results in a specified output being either minimized or maximized by finding the optimal value of a specified input. This can be an enormous timesaver, as calculating the answer analytically for arbitrary geometries can prove to be extremely difficult and time-consuming. In addition, by having an automated loop adjusting the desired boundary condition, much time is saved over manually changing the boundary condition and running the solve routine of the PDE. Of course, the OPT toolkit isn’t restricted to use solely as a companion to the PDE toolkit, but its use as such a companion is definitely one of its more useful and interesting applications.
3.0 Technical Discussion & Results
In this section the progress made by the MigroG team will be discussed. The schedule that the MicroG project has been following can be found in the following Table.
Table 1 -- Project Schedule
|
Design Steps |
Calendar |
|||
|
February |
March |
April |
May |
|
|
Background Research |
||||
|
Literature Research |
||||
|
Compiling/Running of Code |
||||
|
Visualization of Results |
||||
|
Mid Term Presentation |
||||
|
Mid Term Report |
||||
|
Studies of Underlying Methodology |
||||
|
Parametric Studies/Applications 'Matrix' |
||||
|
Supporting Studies on Finite Elements and Optimization ToolKits |
||||
|
MATLAB Design phase |
||||
|
Final Presentation |
||||
|
Final Report |
||||
A day to day description of Table 1 is described below:
Background Research: February 1 – February 10, 1999
Literature Research: February 10 – April 28, 1999
Compiling/Running of Code: February 15 – May 1, 1999
Visualization of Results: March 8 – May 3, 1999
Mid Term Presentation: March 31, 1999
Mid Term Report: April 6, 1999
Studies of Underlying Methodology: April 5 – May 3, 199
Parametric Studies/Applications
'Matrix': April 5 – May 3, 199
Supporting Studies on Finite
Elements and Optimization ToolKits:April 5 – May 3, 199
MATLAB Design phase: April 5 – May 3, 1999
Final Presentation: April 28, 1999
Final Report: May 3, 1999
For the parametric studies the MGFLO code was used. MGFLO is a fluid mechanics and heat transfer iterative code, that uses modified Navier-Stokes equations (discussed in Section 2– Background Theory, of this report). The code allows for variations in the initial conditions as well as boundary conditions. The initial conditions refer to the fluid/thermal parameters that define the fluid that is tested (e.g. water, silicone oil or any other number of fluids). On the other hand, boundary conditions refer to the values that define the boundaries of the domain in which the fluid is (e.g. free surface, temperature conditions, slip/no-slip conditions, flux conditions, etc.). In the test cases that we ran, we chose two different fluids, and several initial and boundary conditions.
3.1.1 Description of Parameters and Domain
3.1.1.1 Fluids and Initial Conditions
The first fluid that we tested was a fictitious fluid in which most of the parameters (initial conditions) were set to unity. With the initial conditions set to unity we started to vary the boundary conditions in order to fulfill a design-space region. The initial conditions for this fictitious fluid were as follows:

All of the parameters, as stated above, are fictitious and without dimensions per se.
The second fluid that we used with the MGFLO code was water. The initial parameters that define water are described below.

With these initial conditions, we started to run cases on the Beowulf cluster (discussed in Section 2.3 in this report) in the CFDLab for different scenarios which included different boundary conditions and different domain sizes.
3.1.1.2 Domain and boundary conditions
For our domain we used a unit cube or deviations from it (e.g. half a cube, quarter size cube –in height, etc.). The domain boundary conditions for the scenarios that we ran are better understood with the visual conception of the domain represented in Figure 5 below.
Figure 5: Description of Domain
There are several different boundary conditions that can be put on the different faces of the domain. Below we will describe all the boundary conditions that we used in our experiments (which are a mere representation of the capabilities of the MGFLO code) and in Section 3.1.2 we will refer to them as we actually use them in the different cases.
Free surface: this term refers usually to the "top" surface (or face 5) that is not bounded by any walls (i.e. free) and in contact with (e.g.) air.
Slip/no-slip condition: this refers to the velocities of the particles in the fluid at the wall. With a slip condition normal velocity is zero and tangential velocity is (generally) non zero, while with a no-slip condition the velocities are zero in all directions at that surface.
No flux: no flux refers to the fact that at that surface there is no perpendicular component of the velocity. (This also refers to a wall were no heat is being lost or gained.)
Mixed BC: this condition refers to a convective condition for the flux of the form h(T – Tref), where h is the convective term, T is the temperature at the surface at that given time, and Tref is a reference temperature (e.g. the surrounding temperature).
Now that these terms have been described, we can look at the different temperature distributions that were used in the test cases. To begin with, the fictitious fluid was used to understand how the MGFLO code worked and also to get familiar with the different boundary conditions and their respective descriptions. To do that we used a linear temperature distribution from face 4 to face 2 of the domain cube. This temperature distribution is better understood from Figure 6 below.
Figure 6: Linear Temperature Distribution
Figure 6 is a cross-section in the X-Z plane, right through the middle of the cube (i.e. half way through the Y axis). The left wall is the "hot" wall, with temperature being 2, and the right wall is the "cold" wall, with temperature being 1.
Another temperature configuration that was used for both the water and the fictitious fluid cases was a circular distribution going from face 4 to face 0. Once again Figure 7 below describes this distribution better, and the figure corresponds to a slice of the cube half way through the Y-axis.
Figure 7: Circular Temperature Distribution
3.1.2 Description of Test Cases and Results
With the initial conditions and the boundary conditions defined, the only parameter that is left to define is the gravity vector. The MGFLO code allows for the use of different gravity vector values. We ran several cases with large gravity fields, low gravity fields (a.k.a. microgravity), and also with zero gravity. The more representative of the cases are described in the following two sections. The cases that are presented in this report were conducted in a zero gravity environment so that the flow would be driven by surface tension instead of buoyancy. What is interesting about surface tension driven flow is that on earth-based experiments (where gravity is present) the flow is mainly driven by buoyancy, and most of its effects have already been studied. On the other hand, surface tension driven flows have not been studied to such a great extent and this was an opportunity to do so. As stated in the introduction fluid/thermal systems developed in microgravity can help grow extremely pure crystals as well as drugs that can help cure cancer or other deadly diseases. Other applications, such as advanced propulsion systems, will also be possible from the understanding of surface tension driven flows.
The first test case that we will present in this report is going to be the most representative and interesting case that we ran while learning the use of boundary conditions and initial conditions to model fluid/thermal systems. During this learning process we used the fictitious fluid that was described in Section 3.1.1.1. To this fluid we changed its surface tension as well as the viscosity. We also varied the boundary conditions between walls, from several "non-dimensional" degrees to just one degree difference. Also the slip/no-slip conditions were varied and the size and dimensions of the cube changed. For all these test cases we used the linear temperature distribution between face 4 and face 2 of the cube (as shown in Figure 6).
As we ran and visualized these cases we interpreted how they affected the flow. We saw that the no-slip condition helped to get the flow moving faster (with less iterations) along the surface were the condition was imposed, while the slip condition kept the flow running slower and further away from the surfaces with that condition. The no-velocity flux condition help the fluid flow stronger parallel to that surface and not at all "through" the surface. The different temperature differences between the walls made the flow move faster/sooner for greater temperature differences, or slower for small temperature differences. And finally, the no flux temperature condition ensured that there was no heat lost or gained through the walls
With these concepts in mind we selected a test problem that produces a clockwise spinning fluid, with a roll-nucleus in the upper half of the cube. We wanted to do this in a short amount of iterations and with reduced gravity. Reduced gravity would ensure that surface tension, and not buoyancy, would drive the fluid motion. With these design concepts we started to set boundary conditions and to slightly change some initial conditions. This fluid, being a fictitious fluid, allowed initial conditions (or fluid/thermal parameters) to be changed.
The boundary conditions for this first design case have a "non-dimensional" unit cube as the domain. Faces 0, 1, 2, 3 and 4 have no slip conditions (which means that the velocities are zero at those walls). Face 5 is a free surface and has a no flux condition (which means that there are no vertical, z, velocities associated with it). Finally, face 2 has a temperature of one and face 4 a temperature of 2 (non-dimensional degree units). The temperature distribution is therefore linear. All the faces without a temperature specification had a no-heat flux condition. These boundary conditions are graphically represented in Figure 8 below. The figure is a X-Z slice of the cube half way along the Y-axis.
Figure 8: Design Case 1 Boundary Conditions
With these initial conditions and boundary conditions for this fictitious fluid in a microgravity environment we expected to see a roll with nucleus in the upper half of the cube. If we take a look at the flow of the fluid in this first design case after 20 time steps (each of .001 seconds) we do see the roll that we expected but closer to the upper surface (face 5) as seen on Figure 9 below.

Figure 9: Design Case 1 Roll Formation
The formation of the roll as seen in Figure 9 above is due solely to surface tension effects. The fluid on the right (face 2) is cold, and therefore has a greater surface tension. The fluid on the left (face 4) is hot, therefore having weaker surface tension. The tendency of the fluid with greater surface tension is going to be to pull the fluid with less surface tension towards it. This motion would be from left (face 4) to right (face 2). Due to the slip boundary condition on the top surface (face 5) and the no-slip condition in the bottom surface (face 0) it is evident that the flow is going to move left (towards the cold fluid) at the top surface. Having a closed domain, the fluid has to go somewhere as it is pulled towards face 2 and therefore is move "downwards" towards face 0. Due to the no-slip condition in face 0, the fluid motion tends to stay away from that surface and therefore completes the clockwise roll in the upper half of the cube as predicted.
As we stated before, this first design case was a representative case that was chosen from the preliminary cases in the initial studies. This case also helped us become familiar with the use of boundary conditions and initial conditions. The use of microgravity and being able to predict (design) the behavior of the fluid made this case particularly interesting and representative.
The second parametric test case, and final case, that we are going to present in this report is also the most representative and interesting from the ones we ran with a real fluid, water. The fluid/thermal parameters that describe water were given in Section 3.1.1.1. To this fluid we could not exactly change any of the specific parameters, but we still could increase and decrease the gravity vector. The surface tension could have been changed if desired by "adding" a few drops of soap to the fluid, which would basically not change any other of the fluid parameters of water. The boundary conditions were significantly different from the ones in Design Case 1. The temperature distribution was circular going from hot on face 4 to cold on face 0, as represented in Figure 6.
We ran several test cases with a full cube (i.e. dimensions 5x5x5 cm3) and with temperature differences ranging from 0.1 degrees to 2 degrees centigrade. We also varied the gravity vector from 0 to 10 earth G’s (i.e. 0 to 98.1 m/s2). All these cases gave us insight on the behavior of such a temperature distribution. We were also aiming for a clockwise roll design in this case. Using all these different boundary conditions we were unable to achieve a roll before we got to steady state. The reason for this was that the calculation fails due to a singularity of the upper right corner. In order to simplify the test problem we added a mixed boundary condition on the top surface (in which there was a no-velocity flux condition, and a thermal flux condition of the form h(T-Tref) that was described in Section 3.1.1.2.).
The final design case that we are presenting here was chosen to be of quarter height (in the Y-direction). The dimension of the cube XxYxZ was therefore 5x0.125x5 cm3. We chose face 4 to have a temperature of 20.1 degrees centigrade, and face 0 to have 20.0 degrees centigrade (making it a 0.1 deg C change). Face 5 had the mixed BC (discussed in Section 3.1.1.2.) where h equaled one and the reference temperature, Tref , (in this case the surrounding temperature) was 20.0 degrees centigrade. These boundary conditions can be viewed graphically in Figure 10 below.
Figure 10: Design Case 2 Boundary Conditions
A slice, half way along the Y-axis, through the domain shown in Figure 10, above, would look like Figure 11 below.

Figure 11: Design Case 2 First Time Step
As we stated before, we expected to have the fluid turn clockwise as time progressed. As predicted (designed) we achieved this result after about 30 time steps (each of about 0.2 seconds). The actual scenario was set to run for 1000 time steps and we found it very interesting to visualize the vorticity in the Y-plane. The visualization of vortices in the Y-plane can give insight to how the fluid is behaving at the local level and therefore explain better how the fluid ends up forming a roll. In the following figures we will see three time steps. The first after about 20 time steps, the second somewhere around the 500th time step, and the last figure somewhere near the 1000th time step. The figures show the temperature contours (in color) as well as the Y-plane vortex contours (black lines). The more spread-out the lines are the stronger the vortex is. Most of the vortex action occurs near the top surface (face 5).
In Figure 12 below we can see how the initial vortex occurs near the upper left corner where the heat starts dissipating.
Figure 12: Design Case 2 Time Step 25
In the following figure, Figure 13, we can see how several vortices have been formed in the upper surface and how the ones to the right have joined the vortices formed at the bottom surface, showing therefore the formation of the roll.
Figure 13
: Design Case 2 Time Step 300
Taking a look at the Y-plane vortices formation towards the 1000th step iteration (i.e. 200 seconds) we see that the top and bottom vortices no longer interact and that the system is therefore arriving at a steady state (Figure 14 next page).
Figure 14: Design Case 2 Time Step 1000
One aspect of these three figures that should not be disregarded are the temperature contour regions. Note how at the beginning, Figure 11, the temperature is evenly circularly distributed (going from face 4 to face 0) and as time progresses, the red plume starts to move across the top surface pushing the cold boundary condition towards the bottom surface. This surface temperature variation affects the surface tension in the fluid and is the main driving factor of this second design case. Due to the different fluid/thermal parameters that govern the way water behaves, it is evident from the last three figures that vortex formation and therefore the characteristic roll takes quite some time to get formed (over 200 time steps). On the other hand, this same design case, tested with the fictitious fluid parameters, yields much quicker vortex and roll formation.
In Figure 15 below, we can see the same case as above but with the fictitious fluid (where all the parameters have been set to unity – defined in Section 3.1.1.1.) The figure shown below shows the second time step (i.e. 0.4 seconds), which compared to the 20th time step in Figure 14 above, shows much more vortex action.
Figure 15
: Design Case 2 Fictitious Fluid Time Step 2
In fact, by the second time step with the fictitious fluid there is almost roll formation and a great deal of vortex formation on the top and left surfaces. The reason for this is obvious when looking at the surface tension and the viscosity of the fictitious fluid. These parameters are set so that the surface tension as well as the viscosity are much higher, therefore making the fluid more prone to roll and vortex formation.
All the numerous cases that we ran in order to be able to design the roll and the different vortex formations could have been greatly enhanced, and the process could have been sped up if an optimization loop existed in the code. Basically, we were the optimization loop. We decided what parameters to change and how, in order to achieve the results that we wanted to get. Ideally and optimization loop should be able to do all of this. For example, if we wanted to increase the vortex formation we could have an optimization loop that calculated the vortex formation after X amount of time steps and if the vortices are not large enough, for example change the surface tension until the optimal value is achieved. Such an optimization toolkit is available in MATLAB, and in the next section a feasibility study was done to see if the same optimization toolkit could be incorporated into MGFLO in the near future so that some optimization automation could speed up the design process of fluid/thermal systems.
3.2 Preliminary Design and Optimization with MatLab
The second phase of the project was to design a system for optimization and feedback as a feasibility study for implementing such a feedback and control system into MGFLO sometime in the future. MatLab was chosen as the programming environment for the project since it has many built-in tools that prove invaluable for rapidly developing an application in controls and optimization. The two toolkits that were principally used were the Partial Differential Equations (PDE) toolkit and the Optimization (OPT) toolkit. Together, these toolkits allow us to create a 2-D or 3-D flux/flow problem (such as heat transfer), and optimize the inputs in order to minimize or maximize a specified output. In the future, we would like to implement such a feedback and optimization system into MGFLO in order to design particular fluid flow patterns as well as heat flow by manipulating the boundary conditions. However, for now we have restricted ourselves to a simpler case in order to validate the concept.
In beginning this phase of the project, it was decided that a problem in heat transfer through a 2-D domain would be easy to rapidly develop, yet would adequately demonstrate the capabilities of the optimization routines as a proof-of-concept for later adaptation and implementation into the larger MGFLO project. Several ideas on how to go about creating a suitable problem were discussed, and in the end a problem dealing with a moving heat source on an edge whose location was to be optimized in order to minimize or maximize the flux leaving another edge was decided upon. In order to arrive at a non-trivial, "interesting" answer, a non-symmetrical geometry for the problem was decided upon.
3.2.1 Background on the Toolkits
Before beginning the actual design of the optimization routines, it was necessary to first research the capabilities of the tools on hand in order to learn what was possible and what wasn’t within the constraints of the programming environment that MatLab provides. The problem relies principally on two toolkits: The PDE and the OPT.
The PDE provides the necessary routines, as well as a graphical user interface (GUI), that allow arbitrary geometries to be set up in various flow/flux problems. Incompressible fluid flow, heat flow, and many other "genres" of problems are available for solving via default equation settings that automatically adjust the governing differential equations and variable names to match those commonly used in each respective discipline. For our design problem, we chose to use the heat transfer setting as that setting provided the correct default environment for specifying the heat transfer-related parameters that we would need for the investigation. The functioning of the PDE toolkit itself is fairly straight-forward. The toolkit takes a specified geometry and creates a finite-element mesh that overlays the geometry, thus dividing the domain into many sub-elements. The fineness of the mesh can be increase by the user in order to get more accurate results at the expense of more processing time. Boundary conditions are then specified by the user at each edge of the domain. In our case, we had the option of specifying a flux condition, a fixed temperature condition, or a combination of the two. After the mesh is created and the boundary conditions set, the steady-state solution of the problem can then be calculated using other routines of the toolbox that solve the governing equations specified by the problem type (heat transfer, fluid flow, etc.). After solving, the user can then visualize the results in a number of ways through color coding and use of vector arrows. In the case of our heat transfer problem, temperature gradients and heat flow could be visualized through the color coding and vector arrows. By manually adjusting the boundary conditions and re-running the solve routine, a parametric study can be conducted on the problem of interest.
The optimization toolkit of MatLab allows a user to take an input and an output function, and let MatLab automatically change the value of the input in search of a maximum or minimum for the output function. This routine works similarly to a root-finding algorithm, in that the input function is generally dependent on a single variable that can be changed, such as a position. The OPT toolkit does not include a specialized GUI, although the results can be visualized through the graphical tools available through the base MatLab application. When used in conjunction with the PDE toolkit, the OPT toolkit allows for a rapid and automated parametric study that results in a specified output being either minimized or maximized by finding the optimal value of a specified input. This can be an enormous timesaver, as calculating the answer analytically for arbitrary geometries can prove to be extremely difficult and time-consuming. In addition, by having an automated loop adjusting the desired boundary condition, much time is saved over manually changing the boundary condition and running the solve routine of the PDE. Of course, the OPT toolkit isn’t restricted to use solely as a companion to the PDE toolkit, but its use as such a companion is definitely one of its more useful and interesting applications.
For the design problem, an irregular, non-symmetric geometry was chosen in order to avoid arriving at a trivial, "boring" solution. After the geometry was decided on and defined in the PDE GUI, the mesh was initialized inside the domain, thereby dividing the domain up into a number of discrete sub-elements. The following figure (Figure 16, next page) illustrates the problem setup after the mesh was initialized.
Figure 16: Domain Layout with Mesh
After the mesh was initialized, boundary conditions were specified on all global domain boundaries. By default, all edges were specified as zero-flux. The boundary condition of the top edge was set via the program to have a temperature distribution that has a small, Gaussian-like shape with a maximum amplitude specified by a position variable. One can think of this as applying a low-grade blowtorch to the top surface, and being able to move it back and forth. Finally, the bottom surface was specified as having a fixed temperature, thereby allowing heat flux both in and out.
With all of the parameters for the problem setup having been specified, the problem setup was saved as a data file to be used by the MatLab design optimization program written for the task. The program was written to read in such a data file and reinitialize the program setup into the GUI using the parameters saved in the file. Most operations are then executed with tools in conjunction with the GUI so that the interim steps of the optimization process can be visualized in addition to the final solution. For a production system, it is possible to write the program in order to utilize the command-line tools exclusively, thereby speeding up the solution process considerably. The final solution could then be visualized for inspection by the user. The program specifies the condition to be optimized for (minimum flux along the bottom boundary) as well as defines the boundary condition of the top surface (the moveable Gaussian-like temperature distribution). The optimization toolkit automatically generates the necessary program loops in order to work through the necessary parameters in order to arrive at the solution. The flux at the bottom surface is calculated by taking the flux solved for by the finite-element portion of the PDE toolkit at each node, and using a simple trapezoidal integration routine to add up all the nodal fluxes into the total boundary flux in order to get the output value used by the optimization toolkit. Some of the interim steps taken by the program are shown in the following graphics: (between Figure 17 and Figure 20, next two pages)

Figure 17: Initial guess for optimizer and the resulting solution
Figure 18: Optimization in progress
Figure 19: More Progress for the optimizer
Figure 20: Final solution for optimizer
As you can see, the optimizer starts from an initial guess for the position of the temperature distribution on the top surface, and moves it back and forth until it finds a solution to the function it is to minimize.
3.2.3 Direction of Future Work
While this example of design work has shown its merit in solving relatively simple cases of heat transfer, it is necessary to put it into perspective on how it fits into the grand scheme of the MGFLO project. At the moment, parametric studies using MGFLO are very tedious as they require manual input of boundary conditions, and then inspection of the graphical interpretation of the numerical results obtained from the computer output. Little in the way of design work can be accomplished with such a manual and time-consuming process. The following figure (Figure 21) demonstrates the current process.
Figure 21: Current Process Flowchart
Since the current manual process is very tedious and time-consuming, it is desired that interpretation of results and adjustment of boundary conditions be automated. In much the same way that the simple design project read the flux and adjusted the position of the temperature distribution heat source, it is hoped that a program can be developed for use with MGFLO. In theory, a user could specify a condition to be minimized or maximized, such as vorticity in a particular region of the domain. Initial guesses for the boundary conditions would be entered, and the rest of the execution and subsequent loops would be automated. The optimizer would have fluid flow algorithms that would be able to calculate such factors as vorticity from the velocity field provided in the solution. It would then automatically update the boundary conditions and rerun the MGFLO code to get the next solution in the series on the way to the "optimal" solution. Of course, this is only a conceptual view, and the implementation of such a feedback system could prove difficult in practice even if the block diagram seems simple enough. However, to be able to automate the manual interpretation of the output data and automatically update the boundary conditions would vastly increase the productivity of a researcher trying to design thermal-fluid systems.
Fluid flow in microgravity conditions, dominated by thermocapillary surface tension effects, is interesting to researchers for a variety of applications. Among these applications are crystal growth, flame propagation, and environmental systems. There are a number of organizations conducting experimental research on the subject, with NASA being one of the more notable. However, experimental research is expensive and time-consuming to conduct. Researchers involved in the field, in both academic and design capacities, need ways to speed up the process of conducting their investigations while leaving physical experimentation as a verification tool for the more interesting cases. That way is through numerical computation. Through the use of powerful computers coupled with sophisticated software, researchers can now simulate migrogravity thermo-fluid flows to a high degree of accuracy, and in far less time than would be required to set up an actual physical experiment.
At the University of Texas, much work in the field of simulating microgravity thermo-fluid flows has been done through the Computational Fluid Dynamics (CFD) Lab, a part of the Department of Aerospace Engineering. The resulting software, called MGFLO, is a finite-element code optimized for parallel execution on a wide variety of massively parallel computer platforms. In the course of our own research, we utilized the Lab’s Linux-based Beowulf system, which is a 16-node parallel cluster made from high-performance commodity PC’s.
Through the use of MGFLO and conducting parametric studies on fluids in the microgravity environment, we studied the behavior of fluid flow as it reacts to heat-driven thermocapillary surface tension effects in lieu of the gravity-dominated buoyancy effects that we are accustomed to seeing. The visualization of the numeric results from the computations provided great insight into the fluid flows and how the surface-tension effects could drive the flow in the entire fluid.
In addition to investigating fluid flow in microgravity conditions under a variety of boundary conditions, a feasibility study into the design of thermal systems through the integration with optimization software was conducted. MatLab was utilized in this phase of the project, as it has many built-in tools that prove invaluable to rapidly developing such an example design project. A simple heat transfer case in which the optimal location for a localized heat source was designed and solved using MatLab. While certainly the heat transfer case is an over-simplification for what would be necessary in order to construct a useful feedback loop for MGFLO, it nevertheless served as a valid proof-of-concept of what is possible.
Appendix A – MGFLO Input Code
Case 1 Example Input Code
title,"ASE 363Q Design Problem"
debug,flush=no
!
! Material Properties
!
material,rho=1.0,mu=1.0,k=1.0,cp=1.0,beta=1.0,gz=-0.01,refu=1.0,refl=1.0,dgammadt=1.0,deltat=0.1,iprint=0
!
! The Mesh
!
mesh,elem=8,8,8,min=0.0,0.0,0.0,max=1.0,1.0,1.0,partition=1,1,2,iprint=1
group=0,ux=0.0
group=0,uy=0.0
group=0,uz=0.0
group=1,t=1.0
group=2,surften=-10.0
group=3,t=2.0
group=4,uz=0.0
face=0,group=0
face=1,group=0
face=2,group=0
face=2,group=1
face=3,group=0
face=4,group=0
face=4,group=3
face=5,group=2
face=5,group=4
end,mesh
write,mesh
!
! The Solve
!
nondim
solve,system="[T,U]",soln_init,iprint=1,user=sioil,2.0,1.0
write,soln
transient,fixed,delt=0.001,iprint=2,nsteps=20,eps=1.0E-6,itermax=5000,writesoln=1,1
write,soln
end,solve
flop,iprint=3
stop
Case 2 Example Input Code
!case2-3
!home/pozo/CFD/mgflo/case2-3.in
title,"ASE 363Q Design Problem"
debug,flush=yes
!
! Material Properties
!
material,rho=998.0,mu=0.001003,k=0.5945256,cp=4179.0,beta=0.0004,refu=1.0,refl=0.02,dgammadt=-0.00017,deltat=0.5,iprint=0
!
! The Mesh
!
mesh,elem=8,8,8,min=0.0,0.0,0.0,max=0.05,0.05,0.0125,partition=1,1,2,iprint=1
group=0,ux=0.0
group=0,uy=0.0
group=0,uz=0.0
group=1,t=20.0
group=2,surften=-0.00017
group=3,t=20.1
group=4,uz=0.0
group=5,qconv=1.0,20.0
face=0,group=0
face=0,group=1
face=1,group=0
face=2,group=0
face=3,group=0
face=4,group=0
face=4,group=3
face=5,group=2
face=5,group=4
face=5,group=5
end,mesh
write,mesh
!
! The Solve
!
nondim
solve,system="[T,U]",soln_init,iprint=1,user=pozo1,20.1,20.0,0.05,0.0125
write,soln
transient,delt=0.1,maxdelt=0.2,iprint=2,nsteps=1000,eps=1.0E-6,itermax=20000,writesoln=1,1
write,soln
end,solve
flop,iprint=3
nondim
stop
Appendix B – MatLab Code
%--------------------------------------------------------------------
% problem.m
%--------------------------------------------------------------------
clear all; global iter; iter=0;
%pde_fig=pdespec;
pde_fig=toto2;
xi=-1:0.05:1;
yi=-1;
%for i=1:5
% x0=-1.+(i-1)*0.3;
% qtot=pdeobjf(x0,xi,yi,pde_fig)
%end
xopt=fmin('pdeobjf',-1,1,[],xi,yi,pde_fig)
function outp=pdeobjf(x0,xi,yi,pde_fig)
global iter ; iter=iter+1;
myboundcond=['feval(''Tdistr'',x,' num2str(x0) ')'];
%%%Set first number in following line to edge where variable
%%%boundary condition applies!!!
pdesetbd(10,'dir',1,'1',myboundcond)
% Solve PDE:
pdetool('solve')
% Export solution and grid
u=get(findobj(get(pde_fig,'Children'),'flat','Tag','PDEPlotMenu'),...
'UserData');
h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu');
p=get(findobj(get(h,'Children'),'flat','Tag','PDEInitMesh'),...
'UserData');
e=get(findobj(get(h,'Children'),'flat','Tag','PDERefine'),...
'UserData');
t=get(findobj(get(h,'Children'),'flat','Tag','PDEMeshParam'),...
'UserData');
[qx,qy]=pdegrad(p,t,u);
qyn=pdeprtni(p,t,qy);
qyi=tri2grid(p,t,qyn,xi,yi);
outp=-trapz(qyi)
%outfile=['solution.bin' num2str(iter)];
%fid = fopen(outfile,'wb');
%fwrite(fid,u,'real*4');
%fwrite(fid,qx,'real*4');
%fwrite(fid,qy,'real*4');
%fclose(fid);
%[imag,cmap]=capture;
[imag,cmap]=capture(pde_fig);
outfile=['plot.' num2str(iter) '.jpg'];
imwrite(imag,cmap,outfile,'Quality','100');
%pause;
function fig=pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',9);
pdetool('snapon','on');
set(ax,'DataAspectRatio',[1 1.5 1]);
set(ax,'PlotBoxAspectRatio',[3 2 1]);
set(ax,'XLim',[-1.5 1.5]);
set(ax,'YLim',[-1.5 1.5]);
set(ax,'XTick',[ -1,...
-0.80000000000000004,...
-0.59999999999999998,...
-0.39999999999999997,...
-0.19999999999999996,...
0,...
0.19999999999999996,...
0.39999999999999997,...
0.59999999999999998,...
0.80000000000000004,...
1,...
]);
set(ax,'YTick',[ -1,...
-0.80000000000000004,...
-0.59999999999999998,...
-0.39999999999999997,...
-0.19999999999999996,...
0,...
0.19999999999999996,...
0.39999999999999997,...
0.59999999999999998,...
0.80000000000000004,...
1,...
]);
% Geometry description:
pderect([-1 1 1 -1],'SQ1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','SQ1')
% Fixed Boundary conditions:
pdetool('changemode',0)
pdesetbd(4,...
'neu',...
1,...
'0',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'neu',...
1,...
'0',...
'0')
% Mesh generation:
setuprop(pde_fig,'Hgrad',1.3);
setuprop(pde_fig,'refinemethod','regular');
pdetool('initmesh')
pdetool('jiggle')
% PDE coefficients:
pdeseteq(1,...
'1.0',...
'1.0',...
'(0)+(1.0).*(0.0)',...
'(1.0).*(1.0)',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setuprop(pde_fig,'currparam',...
['1.0';...
'1.0';...
'1.0';...
'0 ';...
'1.0';...
'0.0'])
% Solve parameters:
setuprop(pde_fig,'solveparam',...
str2mat('0','1000','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))
% Plotflags and user data strings:
setuprop(pde_fig,'plotflags',[1 1 1 1 2 1 1 1 0 0 0 1 1 0 0 1 0 1]);
setuprop(pde_fig,'colstring','');
setuprop(pde_fig,'arrowstring','');
setuprop(pde_fig,'deformstring','');
setuprop(pde_fig,'heightstring','');
%fid = fopen('gridinfo.bin','wb');
%fwrite(fid,p,'real*4');
%fwrite(fid,e,'real*4');
%fwrite(fid,t,'real*4');
%fclose(fid);
fig=pde_fig;
function fig=pdemodel
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',9);
pdetool('snapon','on');
set(ax,'DataAspectRatio',[1 1.5 1]);
set(ax,'PlotBoxAspectRatio',[1.5 1 1]);
set(ax,'XLim',[-1.5 1.5]);
set(ax,'YLim',[-1.5 1.5]);
set(ax,'XTick',[ -1,...
-0.80000000000000004,...
-0.59999999999999998,...
-0.39999999999999997,...
-0.19999999999999996,...
0,...
0.19999999999999996,...
0.39999999999999997,...
0.59999999999999998,...
0.80000000000000004,...
1,...
]);
set(ax,'YTick',[ -1,...
-0.80000000000000004,...
-0.59999999999999998,...
-0.39999999999999997,...
-0.19999999999999996,...
0,...
0.19999999999999996,...
0.39999999999999997,...
0.59999999999999998,...
0.80000000000000004,...
1,...
]);
pdetool('gridon','on');
% Geometry description:
pderect([-1 1 1 -1],'SQ1');
pderect([-1.0000000000000002 -0.60000000000000009 0.39999999999999997 -0.39999999999999997],'R1');
pderect([-0.59999999999999998 0 0.39999999999999997 -0.39999999999999997],'R2');
pdecirc(0.39999999999999997,0.39999999999999997,0.20000000000000001,'C1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','SQ1-R1-R2-C1')
% Boundary conditions:
pdetool('changemode',0)
pdesetbd(14,...
'neu',...
1,...
'0',...
'0')
pdesetbd(13,...
'neu',...
1,...
'0',...
'0')
pdesetbd(12,...
'neu',...
1,...
'0',...
'0')
pdesetbd(11,...
'neu',...
1,...
'0',...
'0')
pdesetbd(10,...
'dir',...
1,...
'1',...
'0')
pdesetbd(9,...
'neu',...
1,...
'0',...
'0')
pdesetbd(8,...
'neu',...
1,...
'0',...
'0')
pdesetbd(7,...
'neu',...
1,...
'0',...
'0')
pdesetbd(6,...
'neu',...
1,...
'0',...
'0')
pdesetbd(5,...
'neu',...
1,...
'0',...
'0')
pdesetbd(4,...
'neu',...
1,...
'0',...
'0')
pdesetbd(3,...
'neu',...
1,...
'0',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')
pdesetbd(1,...
'neu',...
1,...
'0',...
'0')
% Mesh generation:
setuprop(pde_fig,'Hgrad',1.3);
setuprop(pde_fig,'refinemethod','regular');
pdetool('initmesh')
pdetool('refine')
pdetool('jiggle')
% PDE coefficients:
pdeseteq(1,...
'1.0',...
'1.0',...
'(0)+(1.0).*(0.0)',...
'(1.0).*(1.0)',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setuprop(pde_fig,'currparam',...
['1.0';...
'1.0';...
'1.0';...
'0 ';...
'1.0';...
'0.0'])
% Solve parameters:
setuprop(pde_fig,'solveparam',...
str2mat('0','1488','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))
% Plotflags and user data strings:
setuprop(pde_fig,'plotflags',[1 1 1 1 2 1 1 1 0 0 0 1 1 0 0 1 0 1]);
setuprop(pde_fig,'colstring','');
setuprop(pde_fig,'arrowstring','');
setuprop(pde_fig,'deformstring','');
setuprop(pde_fig,'heightstring','');
fig=pde_fig