a. ‘V: .ww». 3.81.4 v... x:.o~y w. ._ 1 viii}... 1 u.H§.5n..e... L a a flaunt... , i 1.? .7. «.1. .11}; IA :1 ...l.. .1... lab 3...: ‘ f5... 1 7 1...: .\. n ‘I ‘x... 'n Ink: 9.2 I'l‘l‘ll. I. H'CHIGAN STA IVERSITY LIBRARIES \lll\llllllllllll\llllllll\ \llllll l 3 1293 01044 5298 ‘\ This is to certify that the thesis entitled THREE DIMENSIONAL MODELLING AND SIMULATION OF THE CURING OF POLYMER COMPOSITES presented by Ananthapadmanaban Sundaram has been accepted towards fulfillment of the requirements for M. S. degree in Chemical Engineering “Wit fl/flfgg Major professor Martin C. Hawley 0-7639 MS U is an Affirmative Action/Equal Opportunity Institution LIBRARY 2 Michigan State. Unlverslty rue: II RETURN BOXto mwciiiékcléu m'fiir mom. TO AVOID FINES Mum on or Moro dd. duo. DATE DUE DATE DUE DATE DUE 1—7- -I_J- ; MSU IaAn Affirmative MM OM Intuition Wane-m THREE DIMENSIONAL MODELLING AND SIMULATION OF THE CURING OF POLYMER COMPOSITES By Ananthapadmanaban Sundaram A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Chemical Engineering 1994 ABSTRACT THREE DIMENSIONAL MODELLING AND SIMULATION OF THE CURING OF POLYMER COMPOSITES By Ananthapadmanaban S undaram Curing of composite materials is an important process in enhancing the properties of the material for its final application. Modeling of the curing process is required to predict the variation of the different properties of the material and help control the process better. A three dimensional dynamic model of the curing process over materials of complex shapes has been developed using the Boundary Fitted Coordinate System for shape modeling. The developed model has been implemented in numerical simulation which can simulate the shape of the material as well as solve the thermal curing process model over the complex shapes. The simulation software has been developed using the C programming language and an output interface has been developed with the MATLAB(R) external library toolbox. The simulation results have been verified with different analytical models for certain cases and a qualitative treatment of different results in the three dimensional domain has also been presented. ACKNOWLEDGMENTS I express my profound gratitude to my advisor Professor Martin Hawley for his guid- ance, encouragement and support throughout the course of this work. I have at different points of time sought the help of a number of people most notably Dr. James McDowell who has provided a number of suggestions that have gone into my work. Dr. Jianghua Wei provided me with an initial understanding of the problem and I am grateful to him for that. I wish to express my thanks to my colleagues in the research group, Valerie Adegbite and Dhulipala Ramakrishna as well as my friend Sanjay Mishra who have helped me greatly on and off the project. I am always indebted to my parents, who have always been my source of strength and inspiration in life. iii TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES LIST OF SYMBOLS 1 Introduction 1.1 Motivation ................................. 1.2 Background ................................ 1.2.1 Previous Modelling Efforts .................... 1.2.2 Fundamental Concepts ....................... 1.3 Objectives ................................. 2 Shape and Process Modelling 2.1 The Process Model ............................. 2.2 The Boundary Fitted Coordinate System .................. 2.2.1 Algebraic Generation Systems ................... 2.2.2 Elliptic Generation Systems .................... 2.3 Process Model Transformation ....................... 3 Numerical Solution 3.1 Finite Difference Formulations ...................... iv vi vii ix @WUJUJ on 11 14 18 21 24 6 7 3.2 The Multi-Grid Method .......................... 3.2.1 Introduction ............................ 3.3 Operators ................................. 3.4 The Multi-Grid Algorithm ......................... Simulation 4.1 Programming Fundamentals ........................ 4.1.1 Data Structures .......................... 4.2 Graphical Output Interface ......................... 4.2.1 The MATLAB External Library .................. Results and Discussion 5.1 Shape Simulation ............................. 5.2 Profile Simulation ............................. 5.2.1 Analytical Verification ....................... 5.2.2 Thermal Cure Simulation ..................... 5.2.3 Numerical Errors ......................... Summary and Conclusions Future Work BIBLIOGRAPHY A User’s Guide A.1 The Simulation Code ........................... A.2 Input Files ................................. Supporting Files 36 41 41 42 46 48 50 50 53 53 57 67 71 75 82 84 87 C Matlab Script Files 98 D C Code For Simulation 107 vi 2.1 3.1 4.1 4.2 4.3 5.1 A.1 A2 A3 LIST OF TABLES fiansformation Relations ........................ 22 The FMV algorithm ........................... 39 Data structure for mesh generation ................... 45 Conduit structure between mesh generation and profile simulation . . 45 Data structure for profile simulation .................. 46 Data For Simulation Runs ........................ 58 Components of Code : MATLAB Script Files .............. 86 Components of Code : C Language Programs ............. 87 Input files for simulation code ...................... 88 vii 1.1 2.1 2.2 2.3 3.1 3.2 3.3 3.4 3.5 4.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 LIST OF FIGURES Overall framework for cure process simulation ............. Boundary Conforming Transformation ................. Interpolation problem in two dimensions ................ Natural to Physical Coordinate Line Mapping ............. Types of finite difference formulations .................. Performance of relaxation methods ................... Effect of relaxing over different grid sizes ................ Restriction and Prolongation Operators ................. The V and F cycles for multi-grid methods ............... Linked-list structure for Multi-grid method .............. Performance of BFCS method on bad initial guess ........... Sections of initial guess .......................... Sections of boundary fitted object .................... Analytical verification for a spherical body with Dirichlet conditions . Analytical verification for a slab with convective boundaries ..... Geometry for simulation runs ...................... Temperature profiles for varying heat transfer coefficient ....... Cure profiles for varying heat transfer coefficient ........... viii 13 15 16 25 29 30 33 39 51 52 52 55 56 59 61 61 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 7.1 7.2 A.1 Effect of addition of convective boundary on temperature profiles . . 62 Effect of addition of convective boundary on cure profiles ....... 62 Effect of unsymmetrical boundary conditions along 7; coordinate line . 63 Effect of unsymmetrical boundary conditions along 5, C coordinate lines 63 Temperature profiles for varying thermal conductivity ......... 65 Cure profiles for varying thermal conductivity ............. 65 Temperature profiles for varying thickness along n-direction ..... 66 Cure profiles for varying thickness along n-direction .......... 67 Effect of time stepping on simulation : Temperature .......... 69 Effect of time stepping on simulation : Cure .............. 69 Adaptive Mesh Generation ........................ 76 Domain Decomposition Method For Shape Modelling ......... 79 Compiling the code ............................ 85 ix LIST OF SYMBOLS Symbol(s) Meaning (x,y,Z) (E, v, C). (i,j,k), (u,v.w) t T To Physical coordinates Natural coordinates Time. Temperature Ambient or surrounding temperature. Initial temperature. Cure or extent of reaction Microwave power absorbed per unit volume Rate of cure reaction. Radius. Element (ij) of shape factor matrix. Shape transformation matrix. Cofactor of element (ij) of M. Finite difference operator Vector of dependent variables. Restriction operator. Prolongation operator. Heat flux Adaptive meshing functions Symbol(s) Meaning J Jacobian of transformation. Density of material. Specific heat. Thermal conductivity Heat transfer coefficient Grid size Number of grid points along {,1}, C. Boundary in physical coordinate system Boundary in natural coordinate system Heat of reaction. Spatial steps in natural coordinates xi CHAPTER 1 Introduction 1.1 Motivation Thermoset polymer composites are processed under a variety of different methods based on the nature of the final products, as well as the characteristics of the different components which make up the composite itself. The common basis for most of these methods involves driving the thermosetting or cross-linking reaction to completion, to obtain better physical and mechanical properties of the cured material. Different phenomena come into play dur— ing the curing process which critically depend upon the process conditions. The reaction is initiated by supplying heat to the material, in the form of thermal heat sources which prop- agate heat by use of temperature gradients, or by focussing heat into the material by use of other energy sources such as microwaves or radio waves. Irrespective of the particular phenomena of heat transfer, the properties of the material being processed could be considered as varying with respect to the extent of the thermoset reaction, as well as the temperature. In practical industrial applications thick section com- posites need to be processed, which have varying temperature and cure gradients along dif- ferent directions which change dynamically as the reaction proceeds to completion. Hence, the critical system parameters of, temperature and cure are dynamic functions in three di- 1 mensional space, and the domain of definition is the geometry of the material. Any process model which describes these variables, does so, by relating the rate of change of these variables to the different physical phenomena at play in the process. These phenomena are in turn related to the different process parameters, depending upon the ac- tual processing methods used. Hence, a generalized model provides the framework for pre- dicting the critical system parameters of temperature and cure independent of the process- ing conditions as well as the geometry of the material. The specifics of these conditions are treated as various inputs, relevant to the process being considered. Two important as- pects of the material affect the nature of the curing process. These are the geometry and the anisotropy of the system with respect to the properties of the material. This work essen- tially deals with the effects of shape and geometry only under isotropic conditions though the properties are allowed to vary as a function of the dependent variables under this frame- work. Controlling the different processes, involves predicting the nature and magnitude of the temperature and cure profiles. The main motivation for this work is the automation of the control of curing processes which is able to take decisions based on generalized predictive methods at a higher level, and base the lower level reactive control on the actual specifics of the process being employed. The next major motivation for the effort, comes from the fact that a generalized model, and, a simulation based on it would essentially encompass all the different phenomena influencing the process, as well as incorporate into it the flexi- bility to include different sub-models. This would retain the specific information about the kinetics and the dependencies of the material properties on process parameters, as well as, the mode of application of heat such as thermal or microwave. Such a simulation would vastly enhance an automated control system which forms an essential unit of an expert sys- tem for the design and processing of polymer composites. It should be noted, that as dif- ferent sub-models become available, like the model for flow, or the model for microwave power absorption, these should be readily incorporated into the generalized model. 1.2 Background 1.2.1 Previous Modelling Efforts The previous modelling and simulation efforts in this area, have been vastly based on one and two dimensional models (Bogetti & Gillespie [1]). These models though sufficient for the purposes of the application they were put to, are insufficient to characterize the com- plexities involved when the processing of thick section composites, of arbitrary shapes are considered . Essentially, by accounting only for one or two dimensions the predictions that are provided by these models do not include information of the gradients in the other direc- tion(s) and hence also the interactions between the gradients in different directions. This phenomenon would be especially amplified in cases when, there is a preferential curing of the composite using different boundary conditions at the different surfaces. Also, in cases where hybrid sources, that is both thermal and microwave are used for the curing process, the different boundary conditions become especially critical though this case has not been handled by the current effort. 1.2.2 Fundamental Concepts The basis of the process model are the energy and material balance equations in three dimen- sions, which are coupled in different ways depending upon the type of process employed for curing, as well as the nature and characteristics of the material being cured. The basic approach for a modelling effort of this type consists of the following steps. 1. Modelling the domain of definition. 2. Model the process in the original domain. 3. Transformation of the process model equations into the modelled domain. 4. Solve the transformed process model. 5. Representation of the solution in the original domain. Normal modelling efforts tend to keep the domain model the same as the physical one or, at least approximately so. This situation is not possible when the domain is complex with irregular boundaries. Extensive information about the exact shape of the domain is pro— hibitive on the time available, and the process models created in the original domain, cannot be easily solved accurately, in their original domain. Hence a method is required which mar- ries the efficiency of the solution or simulation, with the nature of the problem being solved. The coordinates of the model domain should, in essence, reflect the natural dependencies of the process parameters in space. A finite difference approach, tends to model the domain by assuming that the entire do- main is composed of regular patterns of certain shapes such as rectangular, triangular, spher— ical etc, depending upon the accuracy of the representation using such a pattern, as well as the coordinate system used to represent the model equations. A finite element method on the other hand divides the domain into a number of elements that are polynomially depen- dent on regular coordinates. This gives the flexibility of handling complex shapes not suit- ably solved by the finite difference methods. But there is an increase in the computational complexity and once again, very accurate domain representation is necessary for the imple- mentation of the method. The major drawback of the finite element methods is, there is an uniform use of the same type of elements over the entire domain. The use of varying el- ements is not popular and also requires extensive information from the user for an useful implementation. The next major obstacle in using any of the above methods for the model, is the repre- sentation of boundary conditions especially the conditions involving the flow of heat, mass etc. Typically the flow through any element, is determined by evaluating the normal to the element and projecting the flow parallel to it. For an arbitrary shaped body, this involves ap- proximating the normal derivatives and when the dimensionality of the problem increases, the accuracy of such a representation critically depends on the exactness of the boundary coordinate information provided by the user. For processes whose parameters critically de- pend on the nature of such boundary flow conditions, like the curing process, a good process simulation should provide, a flexible and accurate way of representing these conditions nu- merically, without overwhelrnin g the user with the tedious responsibility of coordinate gen— eration for the complex domain. Finally, from a framework aspect any simulation of a model, of the kind being solved in this effort, should include the necessary tools to extend the same to more general condi- tions. In summary, the programming framework should provide constructs, that would help use different sub-models of varying levels of complexity. For example, under the current implementation, the diffusional effects in the mass transfer are considered not limiting and the concentration profile is solely dependent on the reaction under progress. But an addi- tion of the diffusion phenomena, could be done essentially with little work, without going through the entire process of formulation, representation and solution. Since there exists, a fundamental analogy between the various phenomena being solved, i.e mass, momentum and heat , modelling one with the necessary generality, would imply, that any of the others can be implemented in addition to, or independently, under the same programming frame— work. This is a very important aspect of a good simulation and the necessary framework has been provided in the current effort. 1.3 Objectives The previous section brought out the main objectives of the modelling effort. The modelling goals are 1. Evolve a shape modelling framework to model the arbitrary nature of shapes encoun- tered in a curing process. 2. The shape model of the given domain should aid in a simple representation of the boundary conditions over the surface of the domain, as these critically affect the ac- curacy of the simulation solution. 3. Provide a simulation framework using basic constructs, that help to incorporate var- ious phenomena and sub-models. 4. Provide an accurate and at the same time, convenient representation of the process equations for computationally efficient simulation in the modelled domain. 5. Evolve a solution framework, that is both general, in the sense of using different math- ematical tools, and accurate for the problems being solved. Figure 1.1 indicates the overall framework for the simulation effort. All the different aspects represented in bold indicate the work covered under the current effort. Enhancements that could be added to the current simulation framework have also been indicated in the figure. or Thermal cure cycle Process type (Thermal/MW) | Property to state relationships ’ Kinetic expruons Initial geometry Initial mndltiom — _’ ‘— Numerlal Simuhtion Control Mesh sire Structure Final cure sate Retransformatlon to original geometry Temperature New process; . Cure criteria _, Chngmg Geometry ......... Ema}: ““"“"""" Mamm- lraruformed Coordimtes Mesh Generation 4‘ Profile Solver Temperature: Cure ( Transformation ) Density, Spheat, :‘é'fin—ggfig‘fififi'g} Thermal Conductivity EMT” MOM fl Power Model Kinetics Model Flow Model Figure 1.1. Overall framework for cure process simulation CHAPTER 2 Shape and Process Modelling 2.1 The Process Model Consider, a material of arbitrary shape undergoing the process of curing. Heat is supplied to the body in some manner and this causes the startup of the curing reaction within the body. The reaction releases exothermic heat and this further causes changes in the curing pattern within the body. Hence, the entire process of curing is a dynamic one and since, the heat propagation mechanism also involves conduction the curing pattern is also dependent on the location within the material, or in other words the spatial coordinates. There is no flow within the system and shape of the material is constant throughout the process. The basic equation governing this process is the energy balance equation, which accounts for the heat flux entering and leaving a differential element of the body over time. Under the above conditions and assumptions, this equation is written in standard notation as pa??? = —6 - (7+ Q (2.1) In the above energy balance p is the density of the material, c,, is the specific heat, T is the temperature at any point within the material at a given time, t is the time coordinate, q is the heat flux across any differential surface element inside the material, Q is the net rate of heat generation due to different phenomena. The first term in the above equation involves the conduction of heat. The material properties in the above equation can reflect the changes in the curing process by being dependent on the state of the material i.e., its temperature and extent of cross-linking. A Fourier type relation for the heat flux is assumed which relates the energy balance in heat flux terms to the temperature of the system under isotropic conditions. This is then given as .7 = —,.-v"T (2.2) In the above expression the term k is the thermal conductivity of the material and can be a function of the temperature and other process parameters. In the energy balance, Equa- tion 2.1, the net heat generation term, Q specifically depends on the mode of supply of heat to the body. In the case of thermal cure, the application of the heat from the external source to the body is by convection and hence the term Q contains only the exothermic heat due to the cross-linking reaction. In the case of a microwave(MW) application of heat, this term contains two parts. The reaction heat part would be included because the reaction is still present. However, it has to be noted that the form and functionality of the reaction term may be different in the case of microwave curing for the same material, when compared with thermal curing. It is known that excitation of dipoles causes the evolution of heat dur- in g MW curing. This might influence the reaction mechanism itself by affecting the reaction centers at a molecular level leading to a change in the rate equation for the reaction in the MW case. dditionally, there is another term in O which accounts for the energy absorbed due to the microwave. Unlike thermal means, which supply heat directly only to the bound- ing surfaces of the material by convection, the microwave radiation focusses heat into the body. The interactions of the dipoles of the body with the electromagnetic waves causes successive alignment and relaxation of these dipoles leading to an absorption and release of energy. This energy being absorbed from the MW, depends on the dielectric properties of 10 the material which are once again functions of the temperature, and the degree of cure of the material. Hence, the microwave power absorption by the body might be envisioned to be a function of the above mentioned parameters. Expanding the gradient operator in Equa- tion 2.2, substituting into Equation 2.1 and rewriting the heat generation term in terms of its two constituents i.e reaction and MW parts we get the following equation in Cartesian coordinates. 3T 6 8T 6 (9T a or +p(T,X)7'C(—AH) + P,,,(T,X,.1:,y,z,t) (2.3) In the above equation r. is the rate of the thermosetting reaction, expressed as time‘1 , AH is the heat of reaction and is negative for exothermic reaction, P,,, is the power being absorbed from the microwave radiation per unit area and X is the degree or extent of cure. In the case of conventional thermal curing the Pm is zero and this will be assumed through the rest of the Thesis. The reaction rate needs to be accounted, in terms of the different process variables and this is given by an ordinary differential equation as follows. 1 SE : ”(an (2.4) (It The above equation neglects diffusion, and hence explicit spatial dependency of the extent of reaction. The term rc(X,T) is the rate equation of the thermosettin g reaction available as a function of temperature and cure. The coupled Equations 2.3 and 2.4 form the mathematical process model for curing. The above equations are part of the boundary value problem for T and X. The boundary conditions for this problem come from the convective type boundary. The heat exchange at the bounding surfaces of the body being cured can be modelled as a ll convection process. This leads to the following boundary condition, fi~r7=h(T—T,(t)) onI‘ (2.5) where I‘ is the bounding surface of the arbitrary shaped body, ii is the outward normal on the surface I‘, h is the heat transfer coefficient on the surface of the body, T, is the ambi- ent temperature. In deriving the above boundary condition, the effect of radiation has been assumed to be negligible. Equations 2.3, 2.4 and 2.5 complete the process model for the curing process under the assumptions and conditions explained in this section. 2.2 The Boundary Fitted Coordinate System Consider the mathematical model for the curing process as presented in the previous sec- tion. The system of equations is a second order partial differential equation(PDE) coupled to an ordinary differential equation(ODE). The domain of definition for this system is the geometry of the body. Consider now, the boundary conditions to this system of differen- tial equations. They are defined over the entire bounding surface of the body. For a closed simply connected body, this would be the entire outer surface of the body. Since the body is arbitrary in shape, the surface I‘ is also a complex three dimensional surface. Equation 2.5 which is the boundary condition contains a term ii which is the unit outward normal from the surface of the body. For an arbitrary shaped surface, the orientation of this normal vector is also arbitrary. Solution to this system of equations is essentially numerical, due to the nonlinearities involved and lack of functional representations for some of the parame- ters involved. Thus, the solution strategy would basically involve, the discretisation of the domain including the boundary into a number of elements. In such a case, the boundary nor- mal has to be approximated by using the nearest neighbor elements. Such a procedure has to be carried over all the elements on the boundary. This often involves tedious computation 12 which is not very accurate. Besides, the boundary surface also needs to be approximated, usually by use of polynomial fits for the evaluation of the spatial derivatives involved in the boundary condition. These polynomial approximations do not work very well in a number of cases, due to the fact that these approximations themselves have to vary over the sur— face. In other words, different order polynomials may be required at different regions of the boundary for a good approximation, and the decision as to which to use, is not arbitrary. The Introduction chapter explained the different pitfalls of the finite difference and finite el- ement formulations, as far as the shape modelling was concerned. Both these formulations use some, or all of the techniques explained above, for domain and boundary approxima- tion. There is a lot of user effort involved in the decision making for the shape modelling aspect of the problem solving. The criteria used to make the decisions involved in finite element and difference approximations of arbitrary shapes is a complex numerical exercise and no generic technique exists to evolve them. This is by far the most unattractive aspect in using either of these methods for complex shape modelling. Therefore, it is necessary that the complex geometry be represented in a regular manner by a suitable transformation, for the task of solving the process model. The tradeoff in such a case is the additional work, in- volved in re-formulating the process model in the new domain which introduces additional terms in the model. Consider now, the problem of transforming a complex geometry of three dimensions into a regular shape, say a cubical domain. The aim of this representation is to provide a mapping from each point on the surface of the irregular boundary, to a point on the surface of the regular transformed domain. The advantage of the transformation besides the fact that the domain becomes regular is that there is no necessity for normal evaluation on the boundaries of the regular domain. The entire domain being cubic, the normals are along the positive and negative transformed coordinate directions. Thus, given the surface coordinates of the original domain for the process model, a transformed domain is to be generated by some means, which also imposes, continuity along the coordinate directions 13 Original Geometry Boundary Fitted Regular Geometry TI / <——— Physical Coordinates Natural Coordinates Figure 2.1. Boundary Conforming Transformation within the interior. This problem is basically a Boundary Value problem to be solved for the transformed coordinates, also called the Natural Coordinates, given the surface coordinate values of the natural coordinates as a discrete function of the original or Physical Coordi- nates of the system. The resulting solution of natural coordinates as a function of physical coordinates defined over the entire domain is called the Boundary Fitted Coordinate Sys- tem(BFCS). Figure 2.1 gives the pictorial representation of the method. In the analysis to be developed in the rest of the Thesis, the physical coordinates are denoted as (x,y,z) and the natural coordinates are denoted as (5,7),0. The boundary fitted coordinate system technique attempts to generate a set of natural coordinates €(x,y,z), n(x,y,z), ((x,y,z) under certain constraints. As already mentioned, this 14 problem is a boundary value problem and different methods exist which can be applied to obtain a solution. The most common methods which are employed for this purpose are l. Algebraic grid generation methods. 2. Elliptic grid generation methods. Algebraic generation methods use algebraic interpolation techniques along different direc- tions, using boundary information, to obtain a natural coordinate mapping for the given do- main. Elliptic generation methods make use of elliptic equations of the natural coordinates, with the physical coordinates as the independent variables to solve for the same problem. But, the elliptic methods normally require an initial approximate solution which is then re- fined by use of a generating system of elliptic differential equations. That initial guess, is provided by the algebraic generation methods. The next two sections develop the details of the two methods in light of their implementation in the present work. 2.2.1 Algebraic Generation Systems The shape modelling problem requires the generation of a regular rectangular grid, given the correspondence of the grid points on the bounding surfaces of the regular transformed domain, to the original boundary information on the irregular domain. Consider Figure 2.2 which shows the pictorial representation of the problem on a three dimensional surface. The natural coordinates assume integral values from (0,0,0) to (n,m,l) and let (i,j,k) be the coor- dinates of any point on the natural coordinate space spanned by (517,0. The grid points are shown at their corresponding positions in the physical coordinate space. The figure shows a surface over which C is constant but not necessarily a bounding surface. The intermedi- ate points in the figure are the interior grid points on the regular domain. Corresponding to every value (i,j,k) on the regular grid, there exists a value (x(i,j,k),y(i,j,k),z(i,j,k)). The 15 /&/I gu C = constant boundaries (x,y,z) or r specified at all points on boundaries Figure 2.2. Interpolation problem in two dimensions problem is to find the triplets (x,y,z) for every value of (i,j ,k), given the values at all the ex- treme points. Since each of the physical coordinates x, y or z, are exclusive functions of (i,j,k) the problem is nothing, but one of three dimensional interpolation. The method used to perform this three dimensional interpolation in this simulation is called TYansfinite In- terpolation. The basis of this method is independent interpolation along three directions. Consider the three dimensional grid shown in the Figure 2.3. Let (1),, interpolate the vector F(§,n,() along the 7] direction. For the purpose of this interpolation, any standard technique such as cubic splines, B-spline, Hermite interpolation can be used. Transfinite interpolation performs the following operations E1 = our) (2.6) 16 Figure 2.3. Natural to Physical Coordinate Line Mapping E7, = arr—El) (2.7) E; = oar—Eva.) (2.8) 7"(€,77,C) = EI+E72+E3 (2.9) In the steps above, F refers to a vector in the physical coordinate system with components x, y, 2 respectively along the three principal directions. The Figure 2.3 indicates the dif- ferent bounding surfaces on which 25, 7] and C are constant, respectively. Step 2.7 performs interpolation along the C direction across the 11C surfaces. This step yields the value of the dependent vector f’ along the C coordinate lines on all surfaces including the boundaries in- dicated as 1 and 2 in the figure. The next step 2.8 calculates the error (f — E1) in the one l7 dimensional interpolation on the boundaries 1 and 2 and interpolates this error along the 7] direction. It is to be noted, that the sum of the vectors El and E; exactly matches the bound- ary coordinate information on both C=constant and n=constant surfaces. The next step in- terpolates the total error (f —- El — 5;) due to the two previous interpolations, along the C coordinate lines, using the values of this error vector on the boundaries 3 & 4 where C is constant. As before, the sum of all the three vectors El , E}, E3 exactly matches the supplied coordinate information on all the bounding surfaces. Thus the transfinite interpolation, uses standard one dimensional interpolation techniques to perform three dimensional grid gen- eration. The advantage of this technique is that, one can use any known one dimensional interpolation technique depending upon the nature of the shape and the constraints on the continuity involved. The code developed for algebraic grid generation in this effort makes use of linear interpolation as only two points, one on each opposing boundary is specified for every coordinate line. Sometimes, even in a three dimensional algebraic grid genera- tion only the edge and corner point coordinate information is available. In such cases a two dimensional transfinite interpolation is performed on all the six bounding surfaces first and then a three dimensional transfinite interpolation is performed using the different steps as explained above. For the general case of three dimensional problems algebraic grid gener- ation does not provide a good mapping, in the sense that sometimes the interior point map— pings are not accurate and might appear outside the boundaries. Thus, the algebraic grid generation is not a good stand alone technique for mapping physical coordinates to natural coordinates. In the present formulation an elliptic generating system which makes use of an algebraic grid as an initial guess and refines the solution using continuity and boundary constraints on the natural coordinates has been used to perform the grid generation. 18 2.2.2 Elliptic Generation Systems Elliptic generation systems are natural coordinate systems generated by solving for a set (as many as the dimensions involved) of coupled elliptic partial differential equations, con- structed as a boundary value problem, using the boundary physical coordinate information as the constraining boundary conditions. The idea behind the elliptic generation systems is to evolve natural coordinate surfaces along which one coordinate remains constant, which is similar to a problem of generating stream functions. In this case however the boundary conditions on these natural coordinates or streamlines are the boundary coordinate informa- tion of the object. Elliptic generation systems are guaranteed to provide one to one mapping when the generating equations are Laplacian (Mastin & Thompson [2]). The fundamental basis of the Laplace generating equations is the Eulerian equation for the error (Thompson et. a1. [3]) given by the following expression. I:/// 2 (WWW (2.10) 3’ o.=t.n.c where d),- are the manual coordinates. This error is minimized under the following conditions for the natural coordinates given as the Laplacian system of equations, V26 = 0 V21} = 0 VIC _—._ 0 (2.11) The boundary conditions to the above system of equations are as mentioned earlier, the known surface coordinate information of the domain. These are given as 6 : €1‘(Jiay,3)0n F 19 7] = m~(a:,y,z)onI‘ C = (Mar/,2) on F (2.12) It is to be noted that, the functionalities above are normally available only as discrete val- ues at different points. Thus, the accuracy of surface coordinate information determines the spacings of the grid on the transformed domain. The system of Equations 2.11 with the boundary conditions 2.12 form the elliptic generation system of equations. The characteris- tic of the Laplacian system of equations is that, the generated solution has continuous second derivatives and smooth coordinate lines. The system of equations as defined above is defined on the original domain of the body. Since the body is of arbitrary geometry, once again there is a problem of irregular coordinate lines and surfaces. Thus, the equations formulated to avoid the above problem, themselves have to be solved in a similar domain. This is overcome by reformulating the above equa- tions with a change in the dependent variables. As already stated equations 2.11 and 2.12 provide the solutions C(x,y,z), n(x,y,z), C (x,y,z). By switching the dependent and indepen- dent variables the equations can be transformed to be solved for x(C,n,C), y(C ,1;,C ), z(C,n,C). This is the same as reformulating the problem with a change of the coordinate system from the covariant base vectors to contravariant base vectors. However the two coordinate base vectors need not be parallel as it happens in orthogonal coordinate systems. The advantage of this transformation is the fact that the new set of equations are defined on a regular rect- angular geometry, and hence standard numerical techniques can be used to solve these dif- ferential equations, without the need for approximate interpolation or fitting. The necessary analysis that has to be used to obtain this new set of equations for the three dimensional case can be found in [2] . Basically, this involves determining the matrix of the transformation, in terms of the derivatives of the natural coordinates with respect to each of the physical 20 coordinates. Upon such a transformation the following set of differential equations result. and?“ + 2012179, + 20131‘“ + 20231‘“ + 033$“ 1' 0 0113/55 + 20123751; + 2013966 + 202sync + Oaaycc = 0 011265 + 201225,, + 20133“ + .2023an + C1332“ = 0 (2.13) where the aij are the shape factors of the transformation and are defined as 3 011 = Z fimifimj; 1:12.13: j=1.2.3 (2.14) m=l and fijk is the cofactor of the element in the position (j,k) of the matrix M, which is defined by the expression 1‘5 In .T( M: y: yo yc (2'15) ~ ~ 2: an «C The boundary conditions 2.12 are transformed by the change of dependent variables as, :r = .m(C,77.C) onQ y=yn(~£~.n~C) onfl :: = :g({,1),() on Q (2.16) where Q is the regular boundary(six bounding surfaces) of the regular domain as indicated in Figure 2.1. The coefficients in the differential Equations 2.13 are variable as they are functions of the dependent variables, but the nature of the equations is still elliptic. That is, the transformation does not change the nature of the original formulation. Any numeri- cal method used to solve elliptic equations with nonlinear coefficients is applicable for this 21 system of equations. The solution to the above equations yields the values of x, y, 2 at the different grid points on the regular or natural coordinate system. But the original process model is defined over the irregular domain. The mathematical formulation of the process model needs to be transformed into the regular domain before it is solved. The transfor- mation of the different terms in the process model 2.3 and the process boundary conditions (Equation 2.5) depends on the different derivatives involved. The next section details these transformations and presents the process model in the natural coordinates. 2.3 Process Model Transformation The process model presented in Section 2.1 is defined over the physical coordinate system. The next section presented the transformation of the domain from physical to natural co- ordinates. Correspondingly, the process model should also be represented in its equivalent form in the natural coordinate system using the solution to the elliptic generating system (Equations 2.13 and 2.16). The boundary fitted coordinate generation technique actually ob— tains the natural coordinates as functions of the physical coordinates using the contravariant (normal to coordinate surfaces) base vectors of the system. The curvilinear coordinate lines of the three dimensional system are space curves formed by the intersection of surfaces on which one of the coordinates is constant. Thus, along coordinate lines only one coordinate varies and the others are constant. The various operators such as the gradient, the Laplacian, divergence etc can now be redefined in terms of the base vectors. These base vectors are in turn, related to the physical coordinates through the shape transformation matrix M and the shape factors as which were defined in he previous section. Table 2.1 gives the transforma- tion relations for the different operators involved in the original process model in terms of the covariant and contravariant base vectors as well as the relations of these base vectors to the physical coordinates. In the table the variable 5‘ can be C, 17 or C when i = 1, 2 or 3. All 22 Table 2.1. Transformation Relations Name Symbol Transformation/Definition Jacobian J det—M— Covariant base vector (1",- x€.e'i + yes-é + 25.63 Contravariant base vector (1‘ 5 (fine? + Ewe} + @365) Gradient VA 23;, a‘Aé. Divergence(A) V - A 2i, a‘ - Ag. Laplacian V2A 31,- 223:1 ZfflQU-Agé, + 22:1 (V261) Aé’ Time Derivative 2%? % the other variables have already been defined. The exact details of these transformations and their derivations can be found in [3]. The above transformations can be used in Equa- tions 2.3 and 2.5 to obtain the differential equations and boundary conditions in the natural coordinates given by Equations 2.17 and 2.18 below. These equations are more complex than the original model in the sense that additional terms as well as nonlinearities have been introduced. But the nature of the differential equation is still the same and the boundary con- ditions do not contain anything more than the first derivatives. Besides, these equations are defined over a rectangular grid and hence any method employed to solve parabolic PIT~ X)Cp(T, X)?_7_ = 0"(NT€)5 + 022(K3n)n + OSBIKTCM at J when). + mm} + J2 013 {(KT{)( + (KTC )5} + J2 a’3{(“Tfl)c + (“T<)n} + J2 +p(T,X)r(T,X)(—AH) + Pm(T,X) (2.17) 23 th+mleTZanTg = th,(I); 7, =:l:1whenC=00rn th+ 727—61721.” 201ng5, = th,(t); 72 = :t1 when 17 = Oorm th + 72—— “35320ng0 = th,(t); 73 = i1 when C = Oorl (2.18) PDEs over regular domains can be used to solve the transformed set of equations. Once a solution is obtained, it can be mapped back to the original domain using the solution of the shape model. In the above development, it has been assumed that the shape model is time invariant i.e. the shape of the body does not change with respect to time. When flow prob- lems need to be solved this assumption is not valid anymore and the boundaries themselves change with time depending upon the flow parameters such as pressure and viscosity of the system. In such a case, a shape model has to be solved at every step in time with the bound- ary coordinate information calculated from the changing geometry of the body being cured. In this case, the time derivatives of the process dependent variables need to be transformed in the formulation of the process model in natural coordinates. This relation is given as 8A 8A ()1 (a). = (a) W (a?) “‘9’ In the above equation a? is a vector in the physical coordinate system and E is the correspond- ing vector in natural coordinate system. The subscripts attached to the time derivatives indi- cate the variable being held constant in the partial differential equation. The time derivative on the right hand side is at a fixed position in the transformed space i.e, at a given grid point. Similarly the time derivative on the left hand side is at a fixed position in the physical space i.e., the time derivative that appears in the equations of motion. So, in the case of changing geometry the problem can essentially be solved on a fixed grid though the corresponding physical coordinates change over time. CHAPTER 3 Numerical Solution 3.1 Finite Difference Formulations The last chapter dealt with the different differential formulations involved in the shape mod- elling and process simulation aspects of the given problem. Different methods are available for the effective solution to these problems. However the geometry that is of concern in the transformed problem is a regular (cubical) one and hence the finite difference representation affords a simple and effective means of solving the equations over this domain. Several dif- ferent kinds of formulations exist in determining these finite difference representations to be used. Since a finite difference representation essentially divides up the domain into discrete points along the different directions, two different methods exist, depending on the points at which the derivatives are evaluated. In a Cell Centered approach the variable values are assigned to the center of the cells while the derivative values are evaluated on the edges us- ing the cell centered values. In the Vertex Centered approach, the derivative evaluation is performed on the vertices of the cell which are the grid points for this formulation. Fig- ure 3.1 shows a two dimensional representation of the vertex and cell centered approaches. The next step in the finite difference formulation is, the approximation of the differentials using differences. This is done by using a local Taylor’s series expansion about the grid 24 25 CellCenteredSchent Vertex Centered Scheme Figure 3.1. Types of finite difference formulations point under consideration. When the expansions involve only the neighbors on one side of the grid point, one sided differences result. Besides, these Taylor’s expansions are only of first order and hence they have an error term associated with them. But different one sided differences can be added together to eliminate some of the terms in the error to get second or higher order accurate differences. These are the central difference representations which almost always involve the values of the dependent variables on either side of the grid point at which the derivative is being determined. Thus, depending upon the centering of the grid points and the order of the difference formulations, the derivatives are approximated in dif- ferent manners. The differential equations in the current effort have been approximated us- ing a second order, vertex centered approach. The derivatives at the interior points of the 26 domain have been approximated using the following difference formulations. fi+l.j,k - fi—l.j.k filmy, : 2A5 _ fi.j+l.k - fi,j—l.k full“.- "‘ 2A7] fi,'. "' i,', — fcl.—,,. = ’k‘ifmgf” ‘ (3.1) where fé, f,,, fc are the C, n, C derivatives of a three dimensional function f, (i,j,k) is the point at which the derivative is evaluated and A5, A17, AC are the grid spacings along each of those directions. The double derivatives are approximated as, fi+l,j,k - 2ft,j,k + fi—l,j.k f€€|i,j,t; : A{2 fi.j+l.k - 2fi,j,k + fi.j—l.k fnnlg,j,k A712 fi. '.k. - 2fi.',k + fi. .k-l fCCIgch Z J +1 A52 J (3.2) and the cross derivatives are approximated as follows. fi+l.j+l.lc - fi+l,j—l,k - fi—l,j+l.k + fi—l.j-l,k fénlgd'yc : 4A€An _ ft',J+I,k+I_.I.j-I,k+1— fi.j+l.k-l + fi..,—l.k—l foCIch — 4A77AC [fling/c : .fi+1.j.k+l ‘fi—l.j.k:1-IAEA./2+l,j.k-l + fi-l.j.k-l (3.3) Whenever the second derivatives contain a variable coefficient for the first derivative a local average is used to evaluate the coefficient while the above formulation is used to determine the second derivative of the dependent variable at that grid point. Since the partial differ- ential equation also involves a time derivative all the second space derivatives used in the difference model, are averaged over time according to the Crank-Nicolson scheme. This 27 completes the finite difference development for the differential equations at all the interior points. The boundary conditions of the differential problem needs to be transformed to its finite difference approximation. These conditions may or may not involve derivatives depending upon, whether they are Dirichlet or mixed type conditions. The finite difference formula— tions for these two conditions are different depending upon the purpose of this evaluation. In other words, if the boundary condition itself does not involve derivatives but the deriva- tives on the boundary are required for other evaluations, then a simple one sided second order accurate difference is used. The boundary point formulation in such a case for the present problem is as follows. 3fn.1.k _ 4fn-l,j.k 'I' .f71—2.}.k . ft = 2A5 . i = n boundary (3.4) f, = -3fo..-.k +o4fg'j'k - he“; i = o boundary (3.5) Similar equations can be written for the other two derivatives on the boundary. In the case of mixed boundary conditions the boundary conditions themselves involve first derivatives to be evaluated on the boundary. In such cases a fictitious surface of points surrounding the boundary surface on the outside is assumed. The first derivatives to be evaluated at the boundary points are then approximated by normal central differences involving the single interior point neighbor and the fictitious point on the surface outside the boundary, perpen- dicular to the direction of evaluation as in Equation 3.1. The boundary conditions can then be manipulated to get the value of the dependent variable at the fictitious point as a function of the variable values on the boundary and interior points. Now, the differential equations are written for the boundary point in terms of their finite difference formulations. The ap- proximations used above for the second derivatives would now involve function evaluation at points outside the boundary, which are fictitious. The values of the dependent variable 28 at these points have already been established as functions of variable values at the interior and boundary points, and these functionalities are appropriately used to eliminate the ficti- tious point evaluations. Thus, the final set of boundary difference equations obtained in this manner can be added to the existing set of interior point difference equations to solve for the dependent variable values over the entire grid. For, a three dimensional problem involving various cross derivatives, this process of elimination of fictitious points by substitution is a very tedious and time consuming process with lots of opportunities for mistakes, when done by hand. But the basis of any variable elimination process is a simple concept of co- efficient collection. A brief description of how this could be achieved in the coding without cumbersome reformulation is given in the Simulation chapter of the Thesis. 3.2 The Multi-Grid Method 3.2.1 Introduction Solution to a differential equation, involves a number of steps each of which derives from different areas of mathematics. But essentially all differential equations when solved nu- merically are reduced at some stage to an algebraic approximation either by using finite dif- ference or variational formulations. This reduces the differential equations and the bound- ary conditions, into a set of simultaneous linear or nonlinear algebraic equations. The clas- sic techniques that are used to solve these problems involve matrix reduction or inversion, or iterative solution by repeated correction of an initial guess in the case of nonlinear equa— tions. Techniques such as multidimensional Newton-Raphson methods are also popular, but are normally used in a modified form to overcome computational burden involving calcu- lation of derivatives and huge matrix inversions. All the above methods belong to a class of operations called Relaxations. All the methods start with an initial solution or guess and relax the guess over multiple passes or iterations over the given set of equations. This can 29 also be seen as relaxation of an initial error over different iterations. Figure 3.2 shows the |/\/\ /\ AM v \I V“ V Error .......> /\ /\ A n=20; Independent Variable -------- > Figure 3.2. Performance of relaxation methods performance of a typical relaxation method as a function of the number of iterations em- ployed. The figure indicates two essential features. The first one is that, initially the error is highly oscillatory in nature when seen as a function of the independent variable. But as the number of iterations are increased, it becomes smoother and the average amplitude of the oscillations is reduced. Secondly, after a certain number of iterations further relaxations produce negligible change in reducing the remaining components of the error. The expla- nation for this kind of performance is that, the relaxation methods operate efficiently on the highly oscillatory components of the error to reduce them rapidly, until only the smooth components are left. Their performance degrades when they operate on the smoother com- 30 ponents of the error. Different methods have been adopted in the past that reduce this effect of relaxations by modifications to the actual relaxation algorithm. These include using it- erative parameters which act as functions of the number of relaxations and approximating the solution as weighted sums of solutions over successive relaxations. Multi-grid methods are conceptually based on the reduction of the smooth components of the error by ampli- fying their oscillatory nature. Briggs [4], Hackbusch [5], McCormick([6] and [7]) contain different parts of the details for the development and implementation of multi- grid methods for different problems. These references develop the methodology for one and two dimen- sional problems. The extension to three dimensions is fairly straightforward, except some coding detail which is dealt with, in the next chapter. Figure 3.3 shows the effect of reducing /"'\ Grid Size : II V v '> Error 0.....- /\\/ A V Grid Size : 2h Independent Variable -------> Figure 3.3. Effect of relaxing over different grid sizes 31 the number of grid points, in the solution to a set of difference equations. This reduces the number of equations to be solved, for an approximate solution of the original problem. But more importantly, the representation shows clearly that the error which appeared smooth on the finer grid (containing more grid points) is now more oscillatory or contains oscillatory components in the coarser grid (containing lesser number of grid points). In effect, we have amplified the oscillations of the smooth components of the error by using a two grid rep- resentation for the original problem. Hence, any relaxation method initially used to relax the error on the fine grid could be used to do the same on the coarser grid. Multi- grid meth- ods do this by representing the residue or the difference between the exact and approximate solutions on the coarser grid and use the coarse grid relaxation to correct the smoothened solution on the fine grid. Let the finite difference representation of a differential equation do so on a grid of size ’h’. This is the finest grid over which the problem is being solved. The problem including the boundary conditions can be represented in the operator notation as, Lit (uh) 2 fit (36) where the subscript h refers to the grid size. The operator L can be linear or nonlinear in the dependent variable vector u. The term f,, is the finite difference approximation to the source terms or forcing terms (terms which are not functions of the dependent variable). A relaxation of the above set of equations attempts to provide successive approximations to the variables uh. Since the solution is only an approximation there is always a residual error which is given as, 6h = fl. — Lh (Uh) An approximation to this residual is obtained by relaxing the residual problem on a coarse grid say 2h (larger grid size and hence lesser number of grid points), and using that solution 32 to correct the solution to the problem on grid h. The multi-grid algorithm is nothing, but an extension of the two grid algorithm to the coarsest grid. That is, the multi-grid algorithm successively solves the residual problem on coarser and coarser grids until the problem is relaxed on a grid containing only a single interior point to a good enough accuracy. But there still needs to be established a way to go from the solution on one grid to that of another grid. The operators used for this purpose are discussed in the following section. 3.3 Operators The multi- grid method works with three different kinds of generic operators, other than the operators specific to the problem being solved. Two of these operators are used to obtain coarse grid solutions from the dependent variable values on a fine grid and vice—versa. The operator which calculates the coarse grid solution from the fine grid solution, is called the Restriction operator. The other operator which performs the reverse of the preceding oper- ation is called the Prolongation operator (also called Interpolation in some texts). The third operator is the relaxation operator, which can be any one of a number of iterative methods available for difference equations. The exact construction adopted for any of these oper- ators is in general dependent on the problem being solved i.e the order of the differential equation, the nature of non-linearities etc. The trans— grid operators i.e prolongation and re- striction operators are also additionally dependent on the difference formulation used for the solution. Different operators are used depending upon whether the finite difference for- mulation is based on a vertex centered or cell centered approach. Figure 3.4 shows the two- dimensional coarsening and refining of the grid as well as the solution using the restriction and prolongation operators respectively. Hackbusch [5] provides a very detailed discussion on the different kinds of restriction and prolongation operators in two dimensions. The de- velopment to three dimensions is fairly straightforward. Since the problems being solved 33 I 9 F O Q W Restriction Operator i—i—t o o > H—r—r—i '——r—' Prolongation or Interpolation H_+_.._. Operator fl il Grid Size (h) GridSize(2|r) Figure 3.4. Restriction and Prolongation Operators are all second order partial differential equations, the recommended prolongation operator is linear interpolation and the restriction operator is its inverse. For the problem at hand, the following prolongation operator was used. 213mm 2 1112/th V 0 g i,j, k S n, m,l 11.2inij = 0.50 * (uf’J‘k +'rif:fi1,j,k) V (0 < i < n) ugh-+1,“ 2' 0.50 * (nigh. + 1132+“) V (0 < j < m) 113i,2j,2k+1 = 0.50 (1.1ka + iii-23$“) ; 0 < k < l "disgust“ = 0-25 (“f/ji- + “Egan + “iii-+1 + “thin-+1) V (0 < j, k < m» 1) u3¢+1,21,2k+1 = 0.25 (nigh + “3:1,“- + 2133,“, + It:;1’]‘k+1) V (0 < i, k < 71,1) "3i+1.21+1.2k = 0.25 (”flit- + 183th + 1&3in + uhLJ-Jrhk) V (O < i,j < n, m) 34 It __ s‘ l" _ 2h 2h 211 h 211 2h uzr+1.2j+1.2k+1 - 0-130 (ui.j.k + “41.no- + “tour + "arm + Ham-+1).- + ui+l.j.k+1) +0.125(n?3+1‘k+1 +“fiim+1,k+1) V (O < i,j,k < n,m,l) (3.7) where i,j,k are the grid coordinates of any point in space, u is the dependent variable vector of discretized values, h and 2h are the fine and coarse grids used and n,m,l are the num- ber of grid points used along the three principal coordinate directions. The corresponding restriction operator for the current formulation is given as follows. .le _ h ..'_ ..'_ .._ up“ — 1121-31-2“ r — 0, n 07 J —— 0, m. 07 k — 0,1 __ h f1 — “major.- _ It It , [1 f2 — Uzi—1,2j.2k + “2i+1.2j,2t- + uzi,2j—1,2k +2th + all + uh 2i,2j+1,2k "2i.2j,'2k—1 2i.2j,'2k+1 __ h , It . h f3 — ”Zi—l,2j-l,2k + “2.”-1,2j+1,2r- + ”zany-1,21.- +uh. - + uh- . + it"- - (3 8) 21+1,‘2,)+I,‘2k 21—1,2],2k—1 2:-1,2J,2k+1 ' +, II. + h + h ”2i+1,2j,2k—1 u2i+1,2j,2k+1 uz:’,2j—l,2k—1 +llh + uh + :uh 2i,2j—1,2k+1 2i,2j+1,2k—l 2i,2j+1.2k+1 _ , h It It f4 — “2i—1,2j—1,2k—1 + u2i—l.2j-l,2k+l + uzt—l,2j+1,2k—l , It It It +u2i-l,2j+l,2k+l 'i' u2i+l,2j—l,2k-1 + “2i+1,2j—1,2k+1 + Lh + [,h 12t+1,2j+1,2k—1 "‘ 2i+l,2j+l.2k+1 2’1 6 ‘ Ilia-’1‘. 1' (8f1+4f2 +2f3+f4)/I)4 (3.9) Finally, the choice of a relaxation operator has to be made before providing the formulation for the multi- grid method for a set of nonlinear equations. Any standard matrix relaxation algorithm is a candidate. A number of methods exist including the Jacobi iterations, the al- ternating direction implicit method and the Gauss—Siedel method. All these methods have 35 their distinct features and advantages . But, in conjunction with the multi-grid algorithm it has already been noted that the number of relaxations that are performed on a particular grid is small (about five or six). This implies that the solution algorithm is not significantly dependent on the relaxation procedure for its accuracy or computational efficiency. The sta- bility of the multi- grid method is a characteristic that is dependent on the choice of the relax- ation, but all the ones mentioned above are stable for similar kind of problems and hence this is also not a constraining factor on the choice of the relaxation method. The Gauss— Siedel iterative scheme is one of the most popular schemes used as a relaxation operator and the same has been used in the multi-grid implementation for this simulation. Additionally this method can be combined with a local Newton-Raphson(NR) operator to handle highly non-linear equations ([8]). The NR method has not been specifically implemented for the problems solved in this effort, but a tested routine has been included to do this if the need arises at a later date. The Gauss-Siedel method in its implementation provides a number of choices depending on the ordering of the points relaxed. In other words each of these methods relax over a different ordering of the points at each pass. A Red—Black(odd-even) ordering of points was chosen in the implementation. This essentially means, that the odd numbered grid points along each direction are relaxed first followed by the even numbered grid points. This method provides a very simple reduction to the algorithm when applied for a linearized (all coefficients are locally linearized) coefficient finite difference formulation. Consider the second order accurate representation that was proposed in Section 3.1. When derivatives are approximated using that formulation, it is clear that at any point the differ- ence equation is entirely in terms of the nearest neighbors assuming all the coefficients and non-linear terms have been linearized in terms of the nearest neighbors. Thus at any point the dependent variable value at that grid point can be obtained as a function of its nearest neighbors. This implies that all the dependent variables at odd numbered grid points depend on the variable values at even numbered grid points and vice-versa. Hence a complete re- 36 laxation sweep could be performed only over one set of points at a time. This completes the definition and the choice of the different operators which form the backbone of the multi- grid algorithm. Thus, a formal definition could now be made for this algorithm and the vari- ation of the method that has been used to solve the current problem at this stage. The exact methodology is established in the following section. 3.4 The Multi-Grid Algorithm In the previous section the conceptual basis for the multi- grid method was given and a qual- itative structure was established. This section formalizes and details the algorithms that are used in the implementation. The development is standard and could be found in different references already mentioned in the previous section. Each of these texts provides its own structure and modularity to these algorithms. But the algorithm that is to be given has been derived from these various texts by piecing together different structures as required. Multi-Grid algorithms are defined differently for two different types of problems, namely linear and nonlinear difference equations. The essential difference is in obtaining the forcing vector for the coarse grid from the approximate solution on the fine grid. Due to the non-linearities that are present in the formulated problem the nonlinear multi- grid algo- rithm has been implemented. It was mentioned in the previous section that the coarsest grid correction is obtained for a grid with a single interior point. In the case of Dirichlet boundary condition, this implies that the value of the dependent variable at that interior point could be solved for exactly if the equation is linear. If the equation is nonlinear the interior point solution on the coarsest grid can be obtained only approximately by successive iteration. Hence the algorithm which handles a completely nonlinear problem is called the Fully Ap- proximate Storage (FAS) algorithm (Brandt [9]). Consider the difference equation system of the form given in Equation 3.6 including 37 the discretized boundary conditions. The multi- grid method begins with an initial guess or approximation for the dependent variable on the finest grid. The corresponding error in this approximation on the finest grid is given by Equation 3.2.1. A few sweeps of an appropriate relaxation scheme such as the Gauss-Siedel iteration would reduce the error to the smooth components. This relaxation operator is denoted by Sij')(u,f), where r/ refers to the number of sweeps that are performed on the given grid h. Now consider the immediate coarser grid of size 2h. The above error could be approximated on this grid by using a restriction operator as defined above. The restriction operator which restricts a variable vector on grid size h to a grid of size 2h can be denoted as Rh". This provides a first approximation to the error on the grid of size 2h. This can be written as 67-2/1 = RihEn) (310} Thus the forcing term on the coarser grid is given by the approximation f2h : L2];(u‘2h) _ €er (311) The operator L2,, used above refers to the set of difference equations formulated on the grid of size 2h and not an operator obtained by restriction or other means. This implies that the problem on the coarse grid is now defined as that of finding a solution for the system of equations. L2]:(“2h) : f2}; (312) where f2], is as defined by Equation 3.11. There is only one thing that is to be defined and that is the initial guess for 11-; n and this given by a simple relaxation of the smoothed solution on grid h. ’U-zn mm: = Rihhutl (3.13) 38 Now the coarse grid approximation is smoothened by performing relaxations on the coarse grid to obtain the smoothened solution 11-2/1. The fine grid correction is then done as , .._ ll , ‘2/1 to), — uh + [2,, (rt-2;, — Rn (row) (3.14) In the above equation 1.3,, refers to the interpolation or prolongation operator that carries the solution from a coarse grid to the fine grid. The above development provides the basic for- mal structure for a two grid method or algorithm. The same could be extended by obtaining corrections on the grid size 2h by solving for a similarly constructed problem on a grid of size 4h and so on, till the coarsest grid containing a single grid point is reached. Here, one can iteratively solve for the value of the dependent variable as mentioned earlier. This is the essential multi-grid method. Figure 3.5 shows the pictorial representation of how this algorithm works, that is, one passes from a fine grid to a coarse grid smoothening the initial approximations by relaxation, till the coarsest grid is reached and moves from the coarsest to the finest grid correcting the finer grid solutions along the way. This is called a Fully Multi- Grid V (FMV) algorithm indicative of the V shape pass involved. The FMV algorithm is formally defined as n Table 3.1. The FAS algorithm which was mentioned prior to this de- velopment makes use of the FMV cycle at each grid. This means that, both pre-smoothing and correction using the coarser grids are performed on each of the grids, before moving on to the next grid. The graphical representation of this scheme is given in Figure 3.5. It is also called the F—cycle algorithm. In the FAS algorithm one can also control the number of FMV cycles to be performed and hence change the nature of the FAS cycle from the one shown in Figure 3.5. In a general sense, any multi- grid operation is an offshoot or variation of the methods explained in this section, but they operate on different kinds of variables de- pending upon the nature of the problem being solved and the convenience of representation as dictated by its complexity. Besides, use of the right kind of representation distributes the 39 2h ‘11 The V- Cycle (FMV) 21) ‘1: 8h 32:: The F-Cycle Figure 3.5. The V and F cycles for multi-grid methods Table 3.1. The FMV algorithm FMV(h,lIh,f/,) if h is coarsest then u,,=S‘,:°(u,,,f,,) else begin till:S;:(u/Hfh) V211=Rih(uh) f-2n=L2h(uh)'Rih(Ln(Un)'fn) FMV(2h,U2h,f-2h) u,,=u,,+1.’2‘h(v2,,-ri;,) end 40 memory efficiently during simulation and helps to achieve a high degree of computational speed. The next section deals with the different kinds of variable representations used in the simulation and the manipulations that are performed on them to achieve the above goals. CHAPTER 4 Simulation 4.1 Programming Fundamentals The entire code for the simulation has been written in ANSI C with the supporting programs for output graphical generation using MATLAB(R) script programs, which are essentially programs with MATLAB(R) based commands. The development in the previous two chap- ters indicate that, at the outset, any code for the simulation should contain two relatively independent parts, namely 1. Mesh generation 2. Variable profile simulation (Temperature, Cure etc) Though conceptually, the two parts are different in the sense of accomplishing different as— pects of the simulation, both of them involve solving for differential equations of the same kind. The basic difference is that, the mesh generation part aims to solve for a time invari- ant set of partial differential equation and the profile simulation part involves solving for time dependent set of partial differential equations. The domains of definition for the dif- ferential equations are the same but the nature of boundary conditions may differ. For exam- ple, the mesh generation part involves solving Dirichlet boundary conditions as explained in 41 42 Chapter 2, while the profile simulation can involve solution to either Dirichlet or mixed type boundary conditions. Hence, the basic set of routines that have to be constructed to solve the difference equations in both these apparently different parts perform the same functions using different operators and subject to different boundary conditions. Thus, any code written for this purpose, should have the necessary set of generality in its basic set of routines to allow for usage for the two parts as mentioned above. The mesh generation part consists of the following elements 1. Surface approximation 2. Initial volume generation using algebraic methods. 3. Mesh refinement using elliptic generation The surface approximation involves the generation of a surface mesh using algebraic meth- ods, using the edge and or surface coordinate information provided by the user. The elliptic generation method was explained in detail in Chapter 2 . The finite difference representa- tion of the various different terms involved was also presented there. The multi- grid method is used to solve the finite difference representations of the differential equations in both the above parts. This critically determines the various data structures for the representations to be used for the variables to be used in the code. The next section discusses in detail the basic structures used and their relevance to the understanding and implementation of the various methods discussed so far. The Appendices contain the different routines with the notes on the variables used and the details of the exact implementation of the different algorithms. 4.1.1 Data Structures The basic variables to be solved for, in the simulation are the temperature and cure. Besides, the mesh generation step involves determining the physical coordinates (i.e Cartesian) as a 43 function of the natural coordinates. All these variables are hence functions of the natural coordinates, and hence three dimensional in nature at any point of time. The time depen- dency of the temperature and cure variables is taken into account by replacing the current set of variables, for their values at the previous time step but always retaining values one time step old. This is required, because the finite difference representation for the second derivatives in the process differential equations are done using the Crank-Nicholson tech- nique which involves values of the variables at the previous time step. Hence, the basic data structure unit of this code is a three dimensional matrix represented as U(i,j,k) where U is any three dimensional variable under consideration and i,j,k is the grid position in the natural coordinate system. In addition, as explained in Chapter 3 the multi-grid algorithm also imposes some new representation constraints on the code. Besides being three dimensional in nature, the vari- ables to be solved for are also functions of the grid size. The multi-grid algorithm as de- fined earlier requires that all the variables be calculated not only on the finest grid but also on each of the grids up to the coarsest grid containing just one interior point and all the rest boundary points. The logical solution to this problem is to use a series or a linked list of structures, each of the structures corresponding to a particular grid size and containing in- formation about the grid parameters, the values of the required variables as a function of three dimensional space within that particular grid size. The structures are linked in the di- rection of increasing grid size, to facilitate recursive computation of coarser grid variables, in terms of the values of the same variables on the finer grid, initially before relaxation is performed on the grid. Figure 4.1 shows the pictorial representation of a linked-list struc- ture used for the purpose of multi-grid solution of the finite difference equations. Each of the blocks in the figure could be one of the two array structures that are described below. There are two of these structures that are used, one for the mesh generation part, and one for the profile simulation part. The mesh generation code uses the structure in Table 4.1. In L I I I I I I *- PT I l r . r 2.. - I I-JI I I I I | 4|- i I I I I I I kh (Coarsest Grid) Figure 4.1. Linked-list structure for Multi-grid method the above structure ’***u’ is the C notation, for a three dimensional variable. The source terms are calculated for the coarse grid from the variable values and the source term values for the fine grid as explained in Chapter 3. Once the mesh generation is complete, it provides information on the functional depen- dency of the physical coordinates on the natural coordinates. In addition it also provides the local J acobian values at each grid point as well as the shape transformation matrices i.e the a“ at every grid point. Hence a new structure is defined which serves as the conduit passing information between the mesh generation part and the profile simulation part. The point structure accomplishes this objective. In the C language notation this is defined as in Table 4.2. In the above definition, the variable definition D2_array defines the two dimen- sional matrix a at each grid point which contains the shape factor information. The profile simulation part performs the multi- grid algorithm at each time step from an 45 Table 4.1. Data structure for mesh generation structure varray_xyz { int n,m,l; /* Number of grid points along the 3 directions */ double ***x,***y,***z; /* The physical coordinates as three dimensional matrices */ double ***Sx,***Sy,***Sz; /* Source terms in the difference equations, as three dimensional matrices */ struct varray_xyz *next; /* Pointer to the coarser grid */ } Table 4.2. Conduit structure between mesh generation and profile simulation structure point { double ***x,***y,***z; /* Physical coordinate information at each grid point */ struct D2_array ***a; /* Shape factor information at each grid point */ double ***J; /* Local Jacobian evaluated at each grid point */ } initial to a final point, but it uses a slightly different structure. Essentially, the shape infor- mation available from mesh generation via the point structure has to be made available for each of the different grids from the finest to the coarsest. Besides, the profile simulation part involves solving for difference equations in its relaxation part that use time averaged operators and hence information on the previous time values of the variables should also be available. In addition, the process conditions such as heat transfer coefficients have to be made available at each surface point and for each grid size. All these constraints motivates the use of the structure defined as in Table 4.3 for the profile simulation part. The structures 46 Table 4.3. Data structure for profile simulation structure varray { int n,m,l; /* Number of grid points */ struct point *grid_pt; /* Shape information */ double ***multp; /* Time derivative coefficients */ double ***Sf; /* Source term for difference equation */ double ***T,***X; /*Temperature and Cure variables at current time */ double ***Tg,***Xg; /* Temperature and Cure variables at previous time step */ double ***conv; /* Heat transfer coefficient */ struct varray *next; /* Pointer to next level coarser grid */ } seem to carry a lot of computational burden but actually only the temperature and the cure information are updated during every relaxation sweep. All the other information are either constant throughout an entire step or through all the relaxation passes for a given grid. Be- sides the memory allocation is done dynamically that is the memory is allocated only when required and freed when a particular structure is not required anymore. It was stated earlier that the number of relaxation sweeps on a given grid is very small and this greatly reduces the computational time over the traditional single grid methods. 4.2 Graphical Output Interface Mesh generation provides a model of the original domain on which the equations are to be solved. The accuracy of this representation depends on a number of factors including the nature of the shape and the user supplied surface coordinate information. Thus, the user has to be sure that there exists a one to one correspondence between the modelled or transformed 47 domain and the original shape. The boundary fitted coordinate system generated using the elliptic system of Laplacians ensures this, but the accuracy is not guaranteed. Hence, a vi- sual interface has to be provided wherein the user can view the actual shape of the object being modelled, in correspondence to the transformed coordinates. The user should also be able to view different sections to decide if a more accurate surface representation is required for the shape at hand. In addition, when the simulation moves into actually generating the variable profiles there is a necessity to observe the development of these profiles both with respect to time and space. This provides the user with the problem to shape correspondence of the process in a qualitative sense in the least. The simulation though complete for certain purposes is in its developmental stages, from the point of view of the different processes being modelled, as well as the sophistication of the numerical grid generation procedures which can be enhanced. There needs to be a justifiable tradeoff between the nature of prob- lems being solved, the computational efficiency of the grid generation procedure in terms of time and cost and the accuracy of the results required. Such information can be provided only based on extensive testing of the simulation for different problems and verification with experimental results. A graphical output interface achieves the goal of providing the user in- formation of the magnitude and nature of variables at every stage of the simulation. Such an interface has been established using the MATLAB(R) external interface library which pro- vides tools for simultaneous display of variables as the solution evolves. The representation that MATLAB(R) uses for its variables is different from the representation that is being used within the C code for the simulation, the nature of which are dictated by the methods be- ing used and the demands of generality on the simulation. The next section discusses some details regarding the representational differences and the coding that has been included to translate between the two. 48 4.2.1 The MATLAB External Library MATLAB(R) is a mathematical, control and graphics tool which provides a lot of interac- tive tools that can act as interfaces between the user and a C or FORTRAN code that per- forms a complex numerical task. Besides, MATLAB(R) also provides for its own language which operates on these structural units. In addition, it has a vast compilation of routines which perform predefined functions which include manipulation of two and three dimen- sional plots as well as contour maps. Hence MATLAB(R) can perform tasks based on com- mands from external C code through its external library as well as operate using programs written in its own language. The most efficient use of the tools available from this software is made, when there exists a partition between the nature of the conduit created between the C code and the MATLAB(R) workspace, and the kind of tasks controlled by programs writ- ten in the MATLAB(R) language. This partition is the partition between information pass- ing and information processing. In other words the best arrangement is when MATLAB(R) receives only the data passed to it from the code for display and all the structural manipula- tion and control of the graphical procedures such as plotting and sectioning occurs from the MATLAB(R) language based programs. In any simulation, the graphical interface is just a window through which the user is able to perceive at all times the status of the simulation in an easily understandable form. In a simulation of a complex nature, the interface should pro- vide only the minimum demands on the computer time and memory. This has to be borne in mind when evolving criteria for implementation. To centralize the control the, C code was initially given total control over the MATLAB(R) workspace by explicitly specifying the different operations as well as performing different matrix manipulation operations. It was found that from a speed point of view, the partitioning of the tasks as information pass- ing and processing, between the code and the MATLAB(R) script files is much faster than a centralized control of all operations . This is advised for all future implementations using this workspace. 49 Programs based on the MATLAB(R) language are called script files. These script files are functions which take in arguments which are allowed MATLAB(R) structures, and per— form both numerical and graphical operations on them. The basic information containing unit of the MATLAB(R) language is a matrix (utmost two dimensional). All the variables that MATLAB(R) handles are matrices or vectors including string variables. On the other hand the basic structural unit for all the key variables used in the simulation code are three dimensional matrices. The passage from one representation to another is established from within the C code itself. That is, whenever numerical information needs to be passed to the MATLAB(R) workspace the C code manipulates the structure into the appropriate MAT- LAB(R) format (a three dimensional matrix is converted to a vector in this case). Once the information is passed all further operations to be performed on them is done using the MAT- LAB(R) script files. Once the information is made available on the MATLAB(R) workspace, it needs to be manipulated effectively to obtain the relevant graphical output. Since originally three di- mensional matrices are converted into vectors, simple script programs were written to pro- vide independent access along the different coordinate directions as required . Once this is established, further manipulation of the operators is straightforward and are contained in additional script files written for this purpose. For a detailed discussion of the various commands involved and the exact implementation of the MATLAB(R) external library the appropriate manuals [10] & [11] are referred. The Appendices contain details of the dif- ferent script files used and the tasks being performed by them. The next chapter contains different plots generated by the MATLAB interface . In addition, the interface is also capa- ble of generating color maps indicate of temperature profiles and also allows the user who is conversant with MATLAB to manipulate shapes and figures. CHAPTER 5 Results and Discussion 5.1 Shape Simulation The very first step in the profile prediction is the simulation of the shape of the domain. The boundary fitted coordinate system was implemented using the multi-grid method for the finite difference solution. As mentioned in Chapter 2, the shape simulation consists of two parts, the algebraic grid generation and subsequent refinement using the elliptic system of equations. The algebraic generation is often very approximate and serves only as a first guess for the solution to the elliptic generation system. Thus the capability of the shape simulation is critically dependent on the elliptic generating system. In the current effort this capability was tested by providing a bad initial guess deliber- ately over the entire domain with correct coordinate information being provided only on the boundaries. The elliptic generating system was made to take over from that point on- wards. Figure 5.1 shows the results of such a shape generation procedure starting from a bad initial guess. The elliptic generation system performs very well as shown in the figure. It refines the solution at the interior points so that the interior is boundary conforming and also generates a one to one mapping between natural and physical coordinate systems in the process. Figures 5.2 and 5.3 show the section wise comparison between the surfaces of 50 51 I Initial Guess Boundary Fitted Object Figure 5.1. Performance of BFCS method on bad initial guess the initial guess and the corresponding surfaces on the boundary fitted object. In these fig- ures u,v,w indicate the C , r] and C coordinates respectively. These are plots generated by the user interface created using the MATLAB(R) external library. It is evident that the Lapla- cian system of generating equations provides a smooth solution to the shape simulation in terms of coordinate line continuity in the interior. However, the surface coordinates remain at their specified values. The results of the shape simulation include not only the coordinate values but also the shape factors (m,- s) and the Jacobian (J) values at each of the grid points in the natural coordinate system. These values are then used for solving the process model in the natural coordinate system. u = 4 plane 52 v = 4 plane w = 4 plane u = 4 plane Figure 5.2. Sections of initial guess v = 4 plane w- 4p|ane T Figure 5.3. Sections of boundary fitted object 53 5.2 Profile Simulation 5.2.1 Analytical Verification The boundary fitted coordinate system guarantees a one to one mapping under the elliptic generation system of Equation 2.13. But, the solution to this system is numerical and there is limitation on the accuracy of the problem. The original shape modelling problem has been solved using a system of non-linear elliptic system of equations and though apparently there seems to be a good solution to the shape modelling problem the shape factors need to be accurate to obtain a correct solution for the process model in the natural coordinate system. Besides, the boundary conditions of the transformed problem need to be solved accurately too. Consider for example, a uniform grid containing 10 points along the three principal directions in the natural coordinate system. The entire domain contains 1000 grid points out of which 488 are on the surface. Hence accurate solution to the boundary condition is essential for a good solution to the process model especially in the case of mixed boundary conditions. Hence a verification was performed using two simple cases for which there is no internal heat generation and for which analytical solutions are available. These are l. A spherical domain with Dirichlet boundary conditions 2. A slab with mixed type boundary conditions. The spherical domain problem is essentially an one dimensional problem as represented in spherical coordinates. But the three dimensional Cartesian coordinate counterpart for the problem can be solved to get the same solution because all the points on the surface of the body have the same temperature at all times and hence there is no conduction along the 6 and (I) directions of the spherical coordinate system. The problem being solved here is 5. (LT _ (LT 1'2 (7 or) ‘pcpat (5'1) 54 with the conditions T = T,;0_<_7'SR;t=0 T = T0;7'=R;t>0 (5.2) (7T 1— = 0;7‘=0;t>0 ()1' Of the above boundary conditions the symmetry boundary condition at r=0 cannot be repre- sented easily in the transformed process model. But since this is only a symmetry condition it is inherently present in the shape simulation i.e all the shape factors are symmetric w.r.t the radial coordinate. The analytical solution to the above problem is 2 , 00 _ n , . 2 2t T(r,t)=T,-+(To—T,-) (1 +—7§-E( 1) sin(%7-)crp (3‘ 7’ ’7 D (5.3) 2 ":1 n 10ch Figure 5.4 shows the comparison of the analytical and simulation solutions for different ra- dial positions within the material. It is evident that there is small error involved in the simu- lation. The maximum observed error is about 1.5% and there is absolutely no propagation of the error over time. The simulation was carried out at different radii and the error obtained in the other cases were about the same. The error is basically in the dynamics of the problem because over longer periods of time it dies down to zero as is evident from the figure. The numerical simulation was performed using 3 relaxation sweeps before and after correction and an F cycle algorithm as explained in Chapter 3. The next step in the verification is the case of a three dimensional slab with all dimensions comparable, and with convective boundary conditions. One end of the slab on each of the three directions is insulated while the other end is a convective boundary with an applicable heat transfer coefficient. For sim- plicity, it was assumed that all three dimensions were of the same magnitude and the heat transfer coefficient was the same on all the convective boundaries of the body. The problem 55 Sphere r/FI = 0.3 r/R = 0.90 r/FI = 0.67 ........ Simulation titifiti. Analytical I I I I I l l I l O 1 00 200 300 400 500 600 700 800 900 1000 Time (Seconds) ----> Figure 5.4. Analytical verification for a spherical body with Dirichlet conditions in mathematical terms is stated as follows (PT .627" .827" 8T 54 “353+“a—3n+“527="%a H 22 = 0;r=0;Vy,z;t>0 ()3: 'T (,9— : 0;y=0;V;r,:-:;t>0 (5.5) do (Ii—T = 0;Z=0;V.’E,t;t>0 —n(?—T = hC(T—To);;r=L;Vy,z;t>0 0:1: —K.d—T = hJT—To); ;y=L;V;r,z;t>O (5.6) 56 450 - (7 .7.7) o 400 - S 2 a (If 3 a e g -------- Simulation '- 350 - a ........ Analytical 30. ’ i a r r r r r r r -‘ 0 500 1000 1500 2000 2500 3000 3500 Time (Seconds) ----> Figure 5.5. Analytical verification for a slab with convective boundaries 0T —rc0~ = ltC(T—T0);;z=L;V;r,y;t>0 T = T,; t: 0; Vr,y,z (5.7) The analytical solution for the above problem is not obvious but is constructed as the prod- uct of the solutions to three independent solutions of single dimensional problems under variable transformations detailed in Crank [12]. Under these transformations the solution to the above problem is given, as follows. M8 20003 p"—“’ 32 'L Gx(:c,t) = (L) 6.1: (—-———/"h ) 1 ( 12: + 02 '1' (06041671) 1 port 20003 ([37:51) 2 EL ®y(ya t) = 2 2 exp (__m__ 1 (fin; + (1’ + (Y) 008(16171) pCpt m 11 8 57 00 2000.9 (£53) 2 L ®~ z,t : .. . _ It. _ Ah -( ) 112:1 (Hi: + 02 + 0 Co.s(,[i,)t$b( pcpt ) T(;L‘, y, z, t) = To + (T.- — Tn)(r)r(ar, t)(-)y(y, t)(~)z(z,t) (5.8) th where [3,, is the n root of the equation h. L int ’ 71. :': r anon) a 1.. (5.9) The results of the simulation and the analytical solution obtained as above, are shown in Figure 5.5. Once again it is found that the simulation results are accurate and the maximum error is about 2.0%. This indicates that the nature of the error could be due to truncation as it seems to be independent of the nature of the problem. Smaller grid sizes are the solution to this problem but increasing the grid points in all the three directions is computationally inefficient. 5.2.2 Thermal Cure Simulation The last section verified the results of the process model using certain simple cases. Both these cases do not involve any nonlinearities and they were basically heating experiments from the point of view of a chemical process. The advantage of a simulation package is that effects of different process parameter variations can be studied. Since this simulation is only in the developmental stages with regards to inclusion of different sub-models, it is necessary to at least assess the qualitative correctness of the results for a generalized case of thermal cure. A case study to this effect was performed on a general thermal cure pro- cess for an epoxy/ graphite composite using the material property and cure cycle data from literature (Bogetti& Gillespie [1]). The data are given in Table 5.1. In the table the reaction 58 Table 5.1. Data For Simulation Runs Parameter Value p 1.52 g/cm3 c, 0.942 kJ/(gm K) h: 4.457 x 10'3 kW/(Cm K) h 0.0022 kW/(Cm2 K) A, 35.03 x 106/sec A2 -33.5667 x 106/sec A3 3266.67/sec AEI 8.07 x 104 J/mol AE2 7.78 x 10“ J/mol AE3 5.66 x 104 J/mol AH -l98.90 kJ/g rate parameters are for the following rate equation. 1 ‘d—‘ti- = (k1+ o,a)(1— 0)(0.47 — 0); (0 g 0.30) (10 . E _ 1...,(1— 0), (0 > 0.30) k.- = A,e.z:p(—AE.-/RT); (i = 1,2,3) (5.10) The shape on which the simulations were performed is shown in the Figure 5.6(an inverted V). The base case values of the distance between the top and bottom V—shaped surfaces and the length along the C-coordinate direction in Figure 5.6 were, 2.54cm and 7.62cm respec- tively. In the base case the top and bottom V surfaces are convective boundaries while all the other surfaces are assumed to be insulated. This reduces the problem to single dimensional one in natural coordinates. However it is to be noted that the dimensionality is reduced to one, only from the point of view of natural coordinates but not so in a Cartesian sense. Var- ious parameters were then varied from their base case values shown in the Table 5.1, to observe the changes in the temperature and cure behavior at select positions within the ma- Top V Back Front ’w \i Bottome V Left Right Figure 5.6. Geometry for simulation runs terial. In the ensuing discussion temperature and cure profiles imply the time behavior of the two variables. Effect of Boundary Conditions The boundary conditions play a very important role in the nature of temperature and cure gradients observed within the material. The effective heat transfer coefficient given by (th)eff is the critical parameter which determines the nature of the boundary conditions. Three different runs were performed each using a different value of (he/k)...“ . The results are shown as temperature and cure profiles w.r.t time in Figure 5.7 and Figure 5.8. The val- ues z/Z(= 0, 0.76) shown in the figure are the points on the central C C plane of the object, at 60 71:0 and 71:4 natural coordinate values. Several points are to be noted here. An increase in (hC/k).ff leads to decrease in the exotherm within the material. This is expected because, a larger heat transfer at the boundary convects heat away more efficiently compared to a smaller value of the same. However, the rise to the maximum temperature in this case is faster because a high value of the heat transfer coefficient also allows faster convection of heat to the material. For very large values of (hf/km, the boundary conditions become of the Dirichlet type and there is no difference between the ambient and surface temperature of the material. The simulation results are in concurrence with these observations. In a three dimensional simulation especially one involving complex shapes boundary effects are not restricted to a single direction. Different boundary conditions create gradients along differ- ent directions which interact in a complex manner to affect the temperature and cure profiles of the material. To observe this effect, an additional run was performed, with the front and back surfaces of the V shape shown in Figure 5.6, being subjected to convective boundary conditions with the same heat transfer coefficient as in the base case. The results shown in Figure 5.9 indicate that the exotherms are now smaller but are attained quicker than in the base case. This is to be expected, as the addition of a convective boundary plays the same role as an increase in the heat transfer coefficient on a given boundary. However, the exact nature and magnitude of these exotherms is critically dependent on the relative dimensions of the material along these principal natural coordinate directions. To better understand the complexity involved, the simulation was performed with unsymmetrical boundary con- ditions along the different natural coordinate directions. The zero boundary along each co- ordinate line was assumed to be insulated while the other boundary was convective with a (hc/k)eff equal to that of the base case value. The results of this simulation are shown in Fig- ures 5.11 and 5.12. The (u,v,w) directions listed in the figures are the same as (C, 7), C) of the natural coordinate system. The point (4,2,4) is the central point inside the body. Consider Figure 5.11 which shows the temperature profiles along 7) coordinate lines on the central 71C 61 480 \ . —— \ . \ Cure Cycle ‘.> \ 460 ~ 'x. \ Est ‘. \ Fe \ ‘ ‘ I": \ 440 "" '\.\ . ‘ ‘4 \\.\ 1 \ .\\ \ ' . \ 420 ' r“ A ' I L400 . ‘ 5 e 9380 '."- h/k=50:Z/L=O :.i 0.33 h/k=50;z/L=O.76 . T. U g‘sso ...... h/k-roo.z/L-o '1‘ .93 -------- hik- 100;z/L-=O.76 t: 340 00000000 NR = 500 ; z/L = O '; O....... h/k = 500: z/L = 0'76 I 320 300 " 280 I I 1 I l 0 so 1 oo 1 so 20o 250 300 Time (Minutes) ----> Figure 5.7. Temperature profiles for varying heat transfer coefficrent 1 r I ______ _r-:,_ I r 0.9 - J 0.8 — .. 0'7 ” - h/k = so; z/L = o ‘ 0.6i- h/k=50;z/L=O.76 ~ A 1 $05.. ...... h/k=100;z/L=O _ ‘5 0 0,4)- -------- h/k = 100: z/L = 0.76 _ 0.3 — _ 0.2 l— -4 0.1 — a o I I I .l o 50 100 150 zoo 25o Time (Minutes) ----> 300 Figure 5.8. Cure profiles for varying heat transfer coefficient 62 300 480 fi . \ v r r I . _ .-—- " —' v w _/ ‘. " Cure Cycle 460 ~ I ‘. 1’ ‘. ' \ 44o - .’ ‘\ I, .\\ 420 ~ 1 " A ...... i400 £ Base Case :5: 380 '.-.'.'. Z/Z=O.76; S §360 Additional Convective Boundary +9 ...... z/Z = 0.76; 340 320 3oo . ' - 280 1 4 4 ‘ ‘ O 50 1 00 150 200 250 Time (Minutes) ---—> Figure 5.9. Effect of addition of convective boundary on temperature profiles 300 1 1 el.’ 0.9 ~ .’ .4 'I l 0.8 I. .. “I 0.7 ~ .' r .' I A 0'6 h 1‘ Base Case r | _i le=0.76; e 0.5 f p _, 3 i o 0 4 __ j _ Additional Convective Boundary .. i .’ ...... z/Z = 0.76; I . 0.3 "' I} .- .I' 0.2 ~ ,’ _ u' . I.‘ 0.1 - ,:,~‘ .. vi. 0 I I I I L 0 50 1 00 1 50 200 250 Time (Minutes) ----> Figure 5.10. Effect of addition of convective boundary on cure profiles 63 550 ‘l j I U I .-;'\':-. 5’ \\'- i— all \\' 500 I \‘. ~ \ .- 5 ....,, —— \- .._ J \‘V: .4" ”(P A 450 e- ' ‘ i: do” i i g / g 400 [ICUFQ Cyc'e . . . . (u,v,w):(4,2,4) ‘5 .4- -.-.-.- (u.v,w):(4.4,4) 8 .”- - - - - (u.v.w):(4.6.4) a5, ’— 350 300 -‘ 250 I J I I I O 50 100 1 50 200 250 300 Time (Minutes) ----> Figure 5.11. Effect of unsymmetrical boundary conditions along 7) coordinate line 550 j I fir I U . . . . (u,v,w):(2.4.4) -.-.-.- (u,v,w):(6,4,4) - - - - (u.v.w):(4,4.2) ------- (u.v.w):(4.4.6) Temperature (K) --—> 250 I I I I l 0 50 1 00 1 50 200 250 300 Time (Minutes) -———> Figure 5.12. Effect of unsymmetrical boundary conditions along 5, C coordinate lines 64 surface of the body. The exotherm is reduced in magnitude but is reached faster as the grid points move closer to the convective boundary. This is because, the effects of convection are more pronounced for points closer to the convective boundary, than ones farther away for the case in which the other boundary is insulated. Consider the Figure 5.12 which shows the temperature profiles along the E coordinate line and C coordinate lines from the central point of the body. Comparing the temperature-time behavior of the points (2,4,4) and (6,4,4) with the central point (4,2,4), we find, that in the former case the point lies closer to the con- vective boundary in the v-direction while it is farther away from the convective boundary in the u-direction. There seems to be some offsetting of the deviations and the resulting plot for (2,4,4) is almost similar to that of point (4,2,4). In the latter case i.e for point (6,4,4) however, the grid point is closer to the convective boundary in both directions and hence a much lower exotherm. A similar analysis performed on the results for the points (4,4,2) and (4,4,6) shown in Figure 5.12 would indicate that the results would be expected to be simi- lar to those of points (4,4,2) and (4,4,6) respectively, in a qualitative sense. However, it is evident from the results that the offsetting that occurred in the previous case does not hap- pen in this case and there is a pronounced difference in the temperature -time behavior for the point (4,4,2) compared to that of (4,4,4). Similarly, the exotherm observed for the point (4,4,6) is much lower than that observed for the point (4,4,2). The underlying basis for this difference is due to the complex nature of the shape and the relative dimensions along the natural coordinate lines. In thick section composites like the one for which the simulation is performed such effects can only be captured by a three dimensional simulation especially when boundary effects across certain boundaries cannot be neglected. Effect of Conductivity The thermal conductivity of the material determines the magnitude of spatial gradients ob- served within the material. Three values of thermal conductivity including the base case 65 550 r 1 fl 1 1 (9 (D O 30 ,5 450 — i S? E; .3 400 -.-.-.-.- k = k__base/4; 2/2 = 0 g ......... k a k_base/4; z/Z - 0.76 E ----- k = k_base/2: le = 0 i- k - k_base/2; z/Z - 0.76 +++++++++ k - k_b859: Z/Z - 0 ooooooooo k a k_base; z/Z - 0.76 250 1 Cure --> u I 1 O 50 1 OO 1 50 200 250 300 Time (Minutes) ----> Figure 5.13. Temperature profiles for varying thermal conductivity -.-.-.-.- k = k_base/4; 2/2 = 0 —( ......... k - k_base/4; 2/2 =- 0.76 ----- k a k_base/2; z/Z - 0 .1 k =- k_base/2; z/Z :- 0.76 d +++++++++ k - k_base; Z/Z = 0 000000000 k = kflbase; z/Z = 0.76 d I I I I O 50 1 00 1 50 200 250 300 i I I I Time (Minutes) ----> Figure 5.14. Cure profiles for varying thermal conductivity 66 were used with all the other parameters at their base values to study the effects of thermal conductivity. The results are as shown in Figure 5.13. For large values of thermal conduc- tivity, heat is conducted very quickly between the different layers in the material and thus the spatial gradients are small. For lower values of thermal conductivity on the other hand, the heat conduction is slower and hence large differences in the temperature between the different layers occurs, leading to non-uniform curing patterns. These observations are also evident from the simulation results shown in Figure 5.13 and Figure 5.14. Effect of Thickness 550 . r t i 500 ~ Cure Cycle ,5 450 l g g 400 . L - 2.54; z/L — 0.76 ‘9‘; § ...... L—1.27;z/L—0.76 ,‘l’ 350 -------- L — 5.08; z/L — 0.76 300 a l I I I I 2500 so 1 oo 1 50 200 250 300 Time (Minutes) ----> Figure 5.15. Temperature profiles for varying thickness along n-direction Figure 5.15 shows the effect of thickness as defined earlier on the nature of the tempera- ture profiles. As is evident the thickness of the material along different directions also criti- cally affects the heating patterns as well as the nature of temperature and cure profiles within a material. Three different runs were performed in which the thickness of the V shaped com- 67 posite between the top and the bottom V surfaces was varied. The results are shown in Fig- ure 5.15 and Figure 5.16. Increasing the thickness of the material causes larger exotherms. This is because the spatial gradients are large as the conduction phenomena is overshad- owed by the kinetics due to the higher thickness. On the other hand the time to maximum temperature for the material decreases with increase in thickness. This is because the heat conduction takes longer to provide sufficient heat to the different layers which are now far- ther away from the surface. -.-.-.-. L — 2.54; Z/L — 0.76 I I l ' : o 6 -’ ' ‘ 'I l ...... L— 1.27 ;Z/L—O.76 I I l ------- L — 5.08; ZIL - 0.76 I O _—':- .o l I I I 50 1 00 1 50 200 250 300 Time (Minutes) ----> Figure 5.16. Cure profiles for varying thickness along n-direction 5.2.3 Numerical Errors The last section outlined the errors involved in the simulation some simple cases. However, the original problem is highly nonlinear and the only verification that can be performed is through experiment. Nevertheless, the qualitative trends for the general cases give ample 68 evidence for the correctness of the simulation though the exact accuracy cannot be ascer- tained. Even given this limitation, it is evident that the numerical simulation is critically dependent on the time step being used for the simulation. The time step affects the nature of the results by influencing both the stability and accuracy of the simulation. For an arbi- trary nonlinear problem as the one here no theoretical guarantee is available for the stability of any given method. A good method would give reasonably consistent results for differ— ent time steps within the range of its stability. Another factor which affects the accuracy of the results is the grid size. However, for an arbitrary complex shaped body the grid re- finement needs to be performed equally well in all directions for solving the process model accurately. Such a refinement causes a large computational burden which is proportional to the cube of the relative grid sizes in any particular direction. This affects the computational speed of both the shape and process simulations. Hence for a reasonably sized grid, addi- tional accuracy is usually sought by decreasing the time step in the process simulation. To identify the magnitude of influence of changing the time step different runs were performed for a thermal cure simulation with large gradients for three values of the time step within the range of stability. This range was determined by trial and error and it is dependent on the nature of the problem. However, it was found that for the various cases run, a time stepping of a maximum of 10 units in time serves well without problems of instability. The caution that is to be adopted here is, the time stepping cannot be larger then the smallest time of change in the thermal cure cycle. This would lead to certain regions of the cure cycle not being captured in the simulation and hence resulting in spurious profiles. The results shown in Figures 5.17 and 5.18 are for the base case except that the cure cycle and the heat of re- action value have been changed to shorten the region of simulation. The results show that accuracy is not greatly affected by using a large time step in spite of the great amount of nonlinearities involved. The maximum error over the entire range of time and space values in going from a time step of 2 seconds to 10 seconds is about 1.0%. This demonstrates the 69 800 T I T I I m 550 — /, \‘s, z/Z = 0.76 " -> s \ "-.~ xx \x. 500 2/2 = o "-. "‘s _ \\.:’ N /'\ ..~ \__:‘ ~‘\ 1 ‘ ‘ @450 (JUI’W (D ‘5 ’93 §400 ----- dt = 5 secs ‘ ’93 ..... dt = 2 secs 350 -.-.-.-. dt = 10 secs .. 300 d 250 I I l I I O 10 20 30 4O 50 60 Time (Minutes) —---> Figure 5.17. Effect of time stepping on Simulation : Temperature 1 T I I .f’f --’——' fl I ,/ /'/ .’ / 0 9 "" I I .4 i” _ I a O 8 li 2/2 = II 0 7 e I/ a A I1 l i ----- dt = 5 secs m 0.5 — a '5 ..... dt 2 2 secs 0 0.4 P It! -.-.-.-. dt = 10 secs - /i o 3 ~ // .. // O 2 * /' I ~ / -’ / ’ z/Z = o 76 o 1 ~ / / - I ,/ / 0 I___.—-/ .— ’ i 1 i i O 1 O 20 30 4O 50 60 Time (Minutes) ----> Figure 5.18. Effect of time stepping on simulation : Cure 70 ability of the simulation to perform very well within the range of its stability. CHAPTER 6 Summary and Conclusions Curing of polymer materials is an important process in the composites industry because it enhances several properties of the composite material . The end properties of the material are critically dependent on how well the process is controlled as well as the process con- ditions adopted. Thus a general modelling scheme for the process has to take into consid- eration all these factors. The process parameters not only include the ambient conditions but also the inherent properties of the object such as the shape of the material. The shape, as mentioned earlier in the Thesis influences the nature of interactions between the various phenomena at play during the curing process. A differential equation based approach to modelling can include all the important as- pects of the process as well as capture the various physical phenomena at play. The process model developed in this work is based on a system of differential equations. The curing process itself can be either thermal or microwave depending upon the manner in which heat is supplied to the material and the model is general enough to incorporate the curing mech- anism (thermal/microwave). Most of the previous modelling efforts are either one or two dimensional in terms of accounting for the spatial variation of the process variables. Thus they suffer from problems of incomplete process representation in terms of accurately pre- dicting the trends in the temperature and cure profiles. This is especially true for thick sec- 71 72 tion composites and for materials of complex shape. The differential equation based model has been made to include the influence of arbitrary shapes by the use of the boundary fitted coordinate transformation technique. This method has been used to perform the shape modelling whereby the complex geometry of the ob- ject being processed is transformed by a suitable coordinate transform to obtain a natural coordinate system. Such a coordinate system has been generated by solving a system of nonlinear elliptic partial differential equations. Though several techniques could be used to arrive at the boundary fitted coordinates the elliptic generation system is more efficient as it captures the continuity of the coordinate lines across surfaces and it needs minimum user input in terms of making specialized decisions regarding the mesh generation aspect. The elliptic generation systems also readily provide the shape factors and the local jacobians for the coordinate transformation which are then directly used in solving the process model. The shape modelling is fairly general and requires only the surface or boundary information of the coordinates of the object under consideration. The process model that is originally developed in the physical or complex domain is then transformed analytically to another system of differential equations in the natural co— ordinate domain. This transformation casts the original generalized boundary conditions on the regular boundaries defined by the natural coordinate system. Thus the tedious and error prone normal evaluations are avoided on the boundary of the object when calculating for the surface values of the dependent variables. The resulting system of differential equations are basically of the same nature as that of the generating system and hence the same numerical methods that are applicable for the solution of one can be used for the solution of the other. Since all the differential equations are highly nonlinear and the boundary information is almost never available in a closed functional form, numerical methods have to employed to obtain the solution to the system of differential equations in the case of mesh generation as well as temperature profile simulation. Vertex centered finite difference formulations were 73 used to cast the differential equations into a system of algebraic difference equations. The difference equations are generated by the discretization of the regular domain and the approach to solution is by relaxation of an initial guess. Since the degree of discretization governs the accuracy as well as the computational efficiency of relaxation, a method which makes use of multiple gn'd sizes is more effective in reducing the computational burden on a single grid relaxation. Thus multi-grid methods have been used to obtain the solution of the difference equations resulting from the differential formulations for both the shape and the process model. Initially a single grid relaxation method(a1ternating direction im- plicit method) was used to solve the system of difference equation for the shape generation part. When a multi— grid method was implemented for the same problem with a Gauss-Siedel relaxation strategy the computational time was reduced tremendously. The advantage of the multi- grid methods is the flexibility on the choice of the relaxation scheme. Multi- grid methods improve convergence of single grid relaxations by performing fewer relaxations over multiple grids. The same strategy combined with a Crank—Nicolson time based dis- cretization was used for the solution of the energy balance and the material balance differ- ence equations. The simulation code was developed in the C language. Extensive work has been done in the area of boundary fitted coordinate system as well in multi- grid solutions to differential equations. Several versions of the implementation of these techniques exist in the public domain as Fortran codes. The difficulty in all the readily available codes is that the prob- lems they handle are not very generic from a dynamic point of view. The problems they have been developed for, are mostly steady state problems and any dynamic problem has to be recast to use the existing code. This leads to problems in integration. But, the major problem in using the public domain tools is that there is no choice of methods available at every level for the programmer to incorporate into the code. For instance all the boundary fitted coordinate system based methods have been implemented in FORTRAN code using 74 single grid methods as the solution strategies for the difference equations. Besides, the pro- cess methodologies also impose restrictions on the kind of data structures that have to be used. In this regard one seldom finds a thermal cure cycle as part of the boundary condi- tions for instance. Integrating different available codes leads to mismatch of data structures between the component pieces. The most efficient structure from the point of view of a 3D multi- grid method has to be translated into a manipulatable form for the boundary fitted co- ordinate technique and vice-versa. Thus the option of using public domain tools for devel- oping the code for the simulation was found to be less useful and less attractive compared to the development of an indigenous code for the specific purpose at hand. The simulation may involve the use of complex shaped objects. Thus visualization of the results of the shape modelling phase of the simulation is important. The simulation presents the temperature and cure profiles as a function of time. A graphical display of the evolving profiles helps the user to view the trends and decide whether to change the grid pa- rameters for better accuracy. A graphical output interface which displays 3D object images as well 2D plots of temperature and cure profiles and color maps of different sections of the heated object has been developed using the MATLAB external graphics library. Finally, the various segments of the model and the code have been tested using different analytical cases for their accuracy. For a general scenario involving a complex shaped body the qualitative trends of the results of the simulation have been analyzed for cases of different parametric variations and have been ascertained to be accurate for the situations considered. The model as well as the code are in their developmental stages and several issues need to be addressed before they could become powerful predictive tools. An outline and a brief discussion of some of the avenues for future work are discussed in the next chapter. CHAPTER 7 Future Work The previous chapter outlined the capabilities of the simulation as well as demonstrated the accuracy of the method using certain problems for which analytical solutions were avail- able. The main objective of this effort was to provide a generalized framework for predict- in g different properties of a material being cured under different process conditions. The full potential of the different techniques involved in shape modelling have not been exploited in this effort. There are several tradeoffs that need to be considered when doing an extensive numerical shape modelling. Currently a basic modelling approach with allowances for ad- ditional features has been established. There are two main developments that can be made in the immediate future. They are 1. Adaptive mesh generation. 2. Composite and domain decomposed mesh generation Both the above developments are basically offshoots from the current method and essen- tially need the current framework for their implementation. The first of the two develop- ments involves rearranging the mesh without altering the entire grid. Figure 7.1 shows the effect of adaptive meshing. The grids are moved towards a particular coordinate line 104 cally, without rearranging the entire grid pattern. Adaptive meshing is done when accurate 75 76 information is required, for large gradients of the dependent variable in certain regions of the domain. By attracting neighboring coordinate lines to that region, a local refinement of the mesh is obtained in that region, as shown in Figure 7.1. There are several points that Original Mesh D'mtribution Adaptive Mesh Refining Around Grid Point 4 Figure 7.1. Adaptive Mesh Generation need to be considered before actually implementing this method. The first one is that there might not be a necessity to do this because the range of variation in the gradients may not be very large. It is to be noted that adaptive meshing has to be performed at every time step and certain numerical criteria need to be established, like the magnitude of the gradients that has to trigger the re-meshing process. This process is itself iterative and hence consumes some computational power and cuts down on the speed of the method. Usually adaptive 77 mesh generation is performed for steady state problems and hence the re-meshing is done only once. Nevertheless one can foresee the need to do this in the case of microwave curing. Basically, curing of a material in a thermal autoclave is controlled using a thermal cure cy- cle. This control is over time because, the variation in the ambient temperature is w.r.t time. At best, a spatial control is obtained by establishing different heat transfer coefficients at different boundaries. In the case of microwave curing several process parameters could be changed to obtain very accurate spatial monitoring and control of the cure process. Thus any system which performs the control for microwave curing needs to monitor the spatial variation of temperature very accurately. Hence there is certainly a need for a numerical simulation system, which can predict large variations in temperature over space accurately in the local neighborhood of that region. Finally, the implementation of the adaptive mesh- ing procedure needs to be discussed. The shape simulation was performed by solving for the system of Equations 2.13. There are no parameters in this system of equations that involve the distribution of the grid independently without changing the surface coordinate distribu- tion. The adaptive meshing procedure involves solving for the set of equations given by, allies + 20121357, + 2013133 + (122%,, + 2023£qc + 0:331?“ l 1‘? {$EP(63 7’, C) + xllQ(€t 713C) + $(R(€,7I, C)} = 0 0113/66 1' 20123157; + Zarayac + 0223/1221 + 2023ync + Gail/(C 1 , . +.]_2. 3/€[)(€17]1€)+y7)Q(€77’1C)+ch(€77’1Q)} = O “112% + 201225,, '1' 201335( 1' 0222117) + 202.3an + 03:3ch 1 . 1’? Z€P(f, 7” C) + qu<£a 713C)+ 3CR(£i 7]: ()I = 0 (7'1) The boundary conditions are still the same as in Equation 2.16. The above equations include in them three functions P, Q, R which contain the various parameters to distribute the grid in any desired manner throughout the domain. But, the one to one mapping guarantee provided 78 by the Laplacian generating system of Equation 2.13 is not valid any more. But for certain classes of functions P,Q,R the mapping still holds (Thompson [13]). Adaptive meshing is discussed in great detail in [3]. Since the generating system for the shape model is now p(T,X)(Tp(T,X)(T- : ”11(h c)¢+022(“2 l)l+n33(h C)C ()l, J 012{(HT€),,+(“T11)5} + J2 013{(KT€)( + (”Tclg} + J2 a... (on). + <~T<>,,} + J2 +N{T5P(E,71,C) + T.,Q(E, 7], C) + 72315,?!» C)} +p(T, X)7-(T, X)(—AH) + P,,,(T, X) (7.2) Poisson instead of Laplace the process model is also appropriately transformed to include the new functionalities. The process model equations are now written as in Equation 7.2. The rate equation and the boundary conditions of the transformed process model Equa- tion 2.17 still apply in the same form. The above equation has also to be iteratively solved along with the mesh generation equations at every step to obtain the refined solution at the new grid points. The current simulation code includes the terms P, Q, R, to be defined as functions of natural coordinates in its formulation, but, these terms have been set to zero in accordance with the Laplacian system of generating equations. The adaptive meshing pro- cedure could be readily performed once the parameters and the functionalities of P,Q,R are established. The second development in the shape modelling area is the method of domain decom- position. For very complex shapes, which are in reality a composite of different three di- mensional bodies attached together at different boundaries, a single transformation into a composite cube as it is done currently is not very accurate. This problem could be handled 79 by considering a series of contiguous cubes each bearing a unique transformation relation to a part of the original domain. Thus the original domain is decomposed into a number of sub— domains each of which is transformed by the boundary fitted coordinate system. Figure 7.2 shows the pictorial representation of the concept of domain decomposition as it applies to the task of generating a shape model for a three dimensional complex object. Miki & Tak- Domain N 5 Eva/\x: Decomposition Transmitting Boundary Conditions J Original Geometry Figure 7.2. Domain Decomposition Method For Shape Modelling agi [14] have performed a domain decomposition algorithm on three dimensional bodies and the same method could be used to solve the shape modelling problem. But, the major computational burden occurs at the next step, which is the profile generation using the pro- 80 cess model. Since the original domain now consists of different parts each having its own set of shape factors and Jacobian the shape model needs to maintain the continuity of these parameters and this is not automatic. Further, since some of the boundaries of any object or subdomain are contiguous with other boundaries the process parameter conditions along these surfaces needs to be solved iteratively to convergence between the different bodies. Since this step has to be performed at every point in time, this method becomes computation- ally very expensive for an unsteady state problem. In a real time environment where these predictive results are required, this method might involve a lot of time. Besides, the exact way in which the domain is to be decomposed and the number of component domains are decisions to be made by the user. There do not exist, at least at the present time, any gener- alized criteria to make these decisions. Thus, the accuracy of the required solution, the time available for simulation as well as the limitations on the amount of information that could be provided to the code have to be considered before going on to the actual implementation of this method. However, the basis of each transformation explained above is still the same as for a single independent body and hence the basic routines and code developed for this simulation can be used with few or no modifications. From a process point of view there are some enhancements that can be made. Several simplifications have been made in performing the simulation as far the process is concerned and these were outlined in Chapter 2. The most important piece is the microwave power ab— sorption model. A model which incorporates, the dimensionality of the problem, the com- plex nature of shape involved as well as dependence of dielectric properties on temperature and cure needs to be developed and tested. When available, this could be readily incorpo- rated into the model. Several control criteria, such as, mode-switching have been evolved for the microwave curing process. Once a working model for power absorption is available these criteria could be included into the simulation to examine the effectiveness of these criteria for various shapes. The next major effect that has not been included in the present 81 work is the anisotropy of the material properties especially the thermal conductivity. This would greatly affect the formulation of the process model since the thermal conductivity in the vectorial flux balance is now a tensorial quantity instead of a scalar quantity. This would also change the formulation of the transformed process model with the inclusion of cross-derivatives for the temperature as well as derivatives for the Jacobians. The necces- sary formulations in this case are provided in [1]. Since all the necessary primitives such as derivatives of three dimensional quantities are already present in the current framework the simulation for this case can be implemented using the same code with the additional terms in the process model. Finally, the simulation framework provides an interactive link between the different phenomena at play during the curing process. Each of the sub-models that is used is normally tested experimentally, independent from the effects of other phenom- ena. When used together in a simulation module these models may need to be corrected to include the effects of other interacting phenomena. Thus, experimental verification of the results of the simulation, once various sub-models have been incorporated, is necessary for the understanding and control of the curing process by different means. BIBLIOGRAPHY BIBLIOGRAPHY [1] T. A. Bogetti and J. Gillespie, Jr., “Two-dimensional cure simulation of thick ther- mosetting composites,” Journal of Computational Physics, vol. 25. pp. 1—273, 1991. [2] C. Mastin and J. F. Thompson, “Transformations of three-dimensional regions onto rectangular regions by elliptic systems,” Numerische Mathematik, vol. 29, pp. 397— 407, 1978. [3] J. F. Thompson, Z. Warsi, and C. W. Mastin, Numerical Grid Generation : Founda- tions and Applications. Amsterdam: Elsevier Science, 1985. [4] W. L. Briggs, Multigrid Tutorial. Society for Induatrial and Applied Mathematics, 1987. [5] W. Hackbusch, Multi-Grid Methods and Applications. Berlin: Springer-Verlag, 1985. [6] S. F. McCormick, Multi-level Adaptive Methods for Partial Differential Equations. Philadelphia: Society for Industrial and Applied Mathematics, 1989. [7] S. F. McCormick, Ed, MultiGrid Methods. Philadelphia: Society for Industrial and Applied Mathematics, 1987. [8] W. H. Press, Numerical Recipes in F ortran the Art of Scientific Computing. Cambridge University Press, 2 ed., 1992. [9] A. Brandt, “Multi-level adaptive solutions to boundary-value problems,” Math. Comp, vol. 31, pp. 333—390, 1977. [10] The MathWorks Inc., The MAT LAB External Interface Guide for UNIX Workstations. [11] The MathWorks Inc., The MAT LAB User’s Guide for UNIX Workstations. [12] J. Crank, The Mathematics of Difi‘usion. Oxford University Press, 1956. [13] J. F. Thompson, Z. Warsi, and C. W. Mastin, “Boundary fitted coordinate systems for numerical solution of partial differential equations — a review,” Journal of Computa- tional Physics, vol. 47, pp. 1—108, 1982. 82 83 [14] K. Miki and T. Takagi, “A domain decomposition and overlapping method for the gen- eration of three-dimensional boundary-fitted coordinate systems,” Journal of Compu- tational Physics, vol. 53, pp. 319—330, 1984. [15] P. R. Eiseman, “Coordinate generation with precise controls over mesh properties,” Journal of Computational Physics, vol. 47, pp. 331—351, 1982. [16] 1. Woo and G. Springer, “Microwave curing of composites,” Journal of Composite Ma- terials, vol. 18, pp. 387—409, 1984. [17] C. W. M. Joe F. Thompson, Frank C. Thames, “Automatic numerical generation of body-fitted curvilinear coordinate system for field containing any number of arbitrary two-dimensional bodies,” Journal of Computational Physics, vol. 15, pp. 299—319, 1974. [18] J. Douglas, Jr., “Alternating direction implicit methods for three space variables,” Nu- merische Mathematik, vol. 4, pp. 41—63, 1962. [19] W.-H. Chu, “Development of a general finite difference approximation for a gen- eral domain-part 1 : Machine transformation,” Journal of Computational Physics, pp. 392—408, 1971. [20] A. Amsden and C. Hiit, “A simple scheme for generating general curvilinear grids,” Journal of Computational Physics, vol. 11, pp. 348—359, 1973. [21] R. B. Bird, W. E. Stewart, and E. N. Li ghtfoot, Transport phenomena. John Wiley and Sons, 1960. APPENDICES APPENDIX A User’s Guide A.l The Simulation Code The simulation code has been written in C to be executed on a UNIX operating system. Currently the simulation has been run on the SUN and Hewlett Packard workstations. The process of running the simulation on a given workstation consists of two parts. 1. Compilation 2. Execution The compilation step involves the generation of the executable file for the code using the different source files. The execution part is the actual running of the simulation under dif- ferent conditions. Normally compilation is an inbuilt operation for any code. But in this case, UNIX environments are available on SUN and HP workstations which differ in their architecture and computational speeds. Programs compiled under one architecture are not executable in the other. Even the compilation commands are different and the system files to be included in each of these cases reside in different system directories. Besides, the sim- ulation code at this stage is not stand alone in the sense of being complete only with the data specified through the various input files. Functionalities such as the rate equation need to 84 85 sundar:new_mesh> make -f HP.mk cc -g -I/usr/local/matlab/extern/include -c mesh_generation.c cc -g -c memory.c cc -g -c shape_in_out.c cc -g —c matrix_cperations.c cc -g -I/usr/local/matlab/extern/include -c shape_guess.c cc -g -I/usr/lccai/matlab/extern/include —lm -c plot_surfaces.c cc -g -lm -c derivatives.c cc -g -lm -c shape_solve.c cc —g -lm -c multigrid.c cc -g -c profile_simulate.c -I/usr/local/matlab/extern/include cc —g -lm -c properties.c cc -g —c process_conditions.c cc -g -c process_in_out.c cc -g -c solve_cure.c cc —g —c pcwer_model.c cc -g -lm —c temp_grid.c cc ~g -lm -c fictitious.c cc -g -I/usr/local/matlab/extern/include -lm -c plot_profiles.c cc ~g mesh_generation.o memory.o shape_in_out.o matrix_operations.o \ shape_guess.o plot_surfaces.o derivatives.o shape_solve.o multigrid.o \ profile_simulate.o properties.o process_conditions.o process_in_out.o \ solve_cure.o power_model.o temp_grid.o fictitious.o plot_profiles.o \ /usr/locallmatlab/extern/lib/hp700/libmat.a \ —I/usr/local/matlab/extern/include -1m —0 profile_gen Figure A. 1. Compiling the code be specified and these are done at this stage by changing the appropriate routine within the program. Once the code is changed it has to be recompiled and thus this makes compilation an important aspect in the running of the simulation. The UNIX environment provides for certain files and commands which can specify and control the sequence of the compilation operation itself. These are called the make files and have a ’.mk’ extension to their name. Two such files are provided with the simulation, namely SUN .mk and HP.mk for use under the two different architectures. Use of the appro- priate file in the make command compiles the different C programs and links them to form an executable file. Figure A.l shows an example of a compilation session on the SUN work- station. Every time the code needs to be compiled the make process can be used. Another 86 Table A. 1. Components of Code: MATLAB Script Files Script File Comments arrangeshapem Arranges matrix for 3D plot. distributem Distributes matrix for plots between coordinate lines find.pic.m Finds named plot init_display.m Initializes display variables shape_plot.m Generates 3D shapes and maps surf_plot.m Generates surfaces of 3D volume time_tplots.m Generates time plots for time.xplots.m temperature and cure. advantage of make files is the fact that they compile only the updated files and use the ob- ject files from the other previously compiled programs. Thus in a code consisting of many programs only the modified ones are compiled instead of the whole set. The README file provided in the Appendices along with the make files gives the various directories in which the system files i ncluded within the program reside in different architectures (HP/S UN etc). Once the compilation is complete, the executable file is generated and resides in the file pro- file_generate. This filename can be used as a command in the UNIX environment and when used so it executes the simulation code. The user is allowed the option of proceeding to the profile simulation or just stop with the mesh generation aspect. This is essential, because the accuracy of the shape model needs to be ascertained before moving to the predictive part in some cases. The entire code consists of a list of 18 C programs and 8 MATLAB(R) script files (with ’.m’ extensions). Table A.l and Table A2 give the names of the various programs and their area of operation. The details of these programs are presented in Ap- pendix B. As mentioned earlier some of these programs need to be modified by the user, to take care of model changes, and the subroutines that need to be changed in such cases contain comments, that indicate this as well as other routines that need to be changed. 87 Table A.2. Components of Code : C Language Programs Program Comments mesh_generation.c Control program for entire code. shape.in_out.c Reads in surface or edge coordinates and stores in appropriate format. shape_guess.c Performs algebraic guess generation using Transfinite interpolation derivativesc Evaluates different derivatives memory.c Performs dynamic memory allocation and freeing. matrix.c Performs various matrix operations multi-grid.c Contains all routines and operators for the Multi-grid method shape_solve.c Solves elliptic generation system profilesimulatec Control program for profile generation processjn_out.c Read in process and property data process_conditions.c and evaluates process conditions propertiesc and material properties solve_cure.c Calculates cure using rate equation power_model.c Evaluates MW power absorbed (Has to be updated) (Not in use currently) plot-profiles.c Manipulate data structures and plot_surfaces.c and generate plots and maps. temp-grid.c Contains all operations to solve the transformed process model fictitious.c Evaluates boundary conditions using fictitious surface. A.2 Input Files The structure of the code is such that all the user specified data such as surface coordinates, process conditions, rate equation parameters are modified according to the requirements be- fore the code is executed. These is done by the use of input files which are files with a ’.DAT’ extension. The above construction is specifically based on the fact that, in the future the nu- merical simulation module needs to be integrated into an expert control system for the curing process. Data files can easily be updated with property and other information by the control- 88 Table A3. Input files for simulation code File Content PROPERTYDAT Properties of material and kinetic parameters PROCESSDAT Process variables such as heat transfer coefficient, cure cycle information NUMERICALDAT Simulation parameters such as time step, mesh size, number of relaxation sweeps etc. DISPLAY.DAT Plotting options like display time, surface and time plot options. COORD.DAT Surface coordinate specification for complex geometry ling system and this would allow information flow from the different components in the ex- pert system to the numerical module. These files contain different types of data which have to be changed according to user preferences before the simulation is executed. Table A3 gives the list of different input files that are used and the type of information contained in each of them. The Appendices contain a sample of all these data files with the details of the exact sequence of information to be provided to the program. Once the modifications are completed the simulation code can be executed. The input files also contain information regarding the display of the different kinds of plots. Accordingly, the simulation code dis- plays different shape and surface plots of the types shown in Chapter 5. There is minimum user input information to be provided during the execution of the code. There are however, facilities provided for the expert MATLAB user to execute standard MATLAB commands from within the simulation to manipulate the various plots generated. The program provides these options to the user at execution time and whenever shape plots are generated by the program. The construction of all except one of the files is very straightforward. The commented listings in the following pages describe the various input files used and the data presented in 89 each. The surface coordinate information has to be provided in a certain format to be mean- ingfully used by the code for shape modelling. The shape to be generated is envisioned as system of packed surfaces and each of these surfaces is a system of orthogonal or nonorthog- onal coordinate lines. The first step in providing the surface coordinates is to identify the direction of packing of the surfaces and the direction of movement of coordinate lines within each of these surfaces. Thus three coordinate directions are obtained for the surface gener- ation. Along each of the bounding surfaces two of the natural coordinate lines vary and the third is assumed to be constant. The surface information as read by the program is to be provided on six bounding surfaces for the transformed cube. Thus the shape is actually as- sumed to be cut up into six bounding surfaces in any particular manner for providing the surface coordinate information. If C, 7], C are the three natural coordinates, then the input file is constructed in three steps as follows. The variables n, m, l are the number of natural coordinate grid points along C, 7] and C directions respectively. 1. C=0,n-1 (7] =0 to m-l [C =0 to MD 2. 7] = 0,m—l (C = 0 to H [C =0 to n-1]) 3. C = 0,1-1 (C =0 to n-l [7} =0 to m-1]) In the above steps the brackets indicate nesting of the coordinate information to be provided to the program. For example, the first step implies that at either end of the C coordinate line, for every point along the 7] coordinate line, the surface coordinate information at all the C coordinate lines has to be provided. The ordering of this information is important and has to follow the sequence shown in the above steps. It is suggested that while generating the coordinate system a right hand coordinate orientation be maintained i.e., move along coordinate directions in such a manner, as to keep the surface being generated to the right. Once the input files are updated the simulation can be executed as explained in the previous section. APPENDIX B Supporting Files B.1 Make File Setup W MakaEilefihanees The main system dependent change that needs to be nrnade to the make file is the resident directories and programs for the MATLAB external libraries. 8111108.: IN CLUDE.PATH = -I/usr/local/matlab/extern/include libmat.a resides in lusr/local/matlab/extern/lib/sun4/ engineh resides in /usr/local/matlab/extern/include/ HELIX; INCLUDE_PATH = -I/usr/local/matlab/extern/include libmat.a resides in /usr/local/matlab/extern/lib/hp700 engineh resides in lusr/local/matlab/extem/include Salads; INCLUDE_PATH = -I/opt/matlab4.2/extern/include libmat.a resides in lopt/matlab4.2/extern/lib/solZ/libmat.a engineh resides in/opt/matlab4.2/extern/include NOTE: In all the above search paths the directory 'matlab' can be any version of matlab for e.g 'matlab4.0a' or 'matlab4.2' etc. 90 91 W CFLAGS: -g LIBS: -lm IN CLUDE_PATH = -I/usr/local/matlab/extern/include all: prof_compile prof_compile: mesh _generation.o memory.o shape_in_out.o matrix_operations.o \ shape _g'uess.o plot_surfaces.o derivatives.o shape_solve.o multigrid.o \ profile_simulate.o propertieso process_conditions.c process_in_out.o \ solve_cure.c power_model.o temp _grid.o fictitiouso plot__profiles.o ; $(CC) $(CFLAGS) mesh _generation.o memory.c shape_in_out.o matrix_operations.o \ shape _g'uess.o plot_surfaces.o derivatives.o shape_solve.o multigrid.o \ profile_simulate.o propertieso process_conditions.o process_in_out.o \ solve_cure.c power_model.c temp _grid.o fictitious.o plot_profiles.o \ lusr/local/matlab/extern/lib/sun4/libmat.a $(INCLUDE_PATH) $(LIBS) -0 profile _gen profile_simulate.o: profile_simulate.c lusr/local/matlab/extern/lib/sun4/libmat.a grid.h $(CC) $(CFLAGS) -c profile_simulate.c $(INCLUDE_PATH) properties.o: propertiesc $(CC) $(CFLAGS) $(LIBS) -c propertiesc process_conditions.c: process_conditions.c $(CC) $(CFLAGS) -c process_conditions.c process_in_out.o: process_in_out.c $(CC) $(CFLAGS) -c process_in_out.c solve_cure.o: solve_cure.c $(CC) $(CFLAGS) -c solve_cure.c power_model.o: power_model.c $(CC) $(CFLAGS) -c power_model.c temp _grid.o: temp _grid.c grid.h $(CC) $(CFLAGS) $(LIBS) -c temp _g"rid.c 92 fictitiouso : fictitiousc grid.h $(CC) $(CFLAGS) $(LIBS) -c fictitiousc plot _profiles.o: plot_surfaces.o plot_profiles.c grid.h \ /usr/local/matlab/extern/include/engine.h $(CC) $(CFLAGS) $(INCLUDE_PATH) $(LIBS) plot_surfaces.c -c plot_profiles.c multigrid.o: multigrid.c grid.h $(CC) $(CFLAGS) -c multigrid.c shape_solve.o : shape_solve.c grid.h $(CC) $(CFLAGS) $(LIBS) -c shape_solve.c derivativeso : derivativesc $(CC) $(CFLAGS) $(LIBS) -c derivativesc plot_surfaces.o : plot_surfaces.c /usr/local/matlab/extern/include/engine.h $(CC) $(CFLAGS) $(INCLUDE_PATH) -c plot_surfaces.c shape _guess.o: shape _guess.c $(CC) $(CFLAGS)$(INCLUDE_PATH) -c shape _guess.c matrix_operations.o : matrix_operations.c $(CC) $(CFLAGS) ~c matrix_operations.c shape_in_out.o : shape_in_out.c $(CC) $(CFLAGS) -c shape_in_out.c memory.o : memory.c grid.h $(CC) $(CFLAGS) -c memory.c mesh _generation.o : mesh _generation.c gridh /usr/local/mat1ab/extern/include/engine.h $(CC) $(CFLAGS) $(INCLUDE_PATH) -c mesh _generation.c 93 3.2 Input Files WAT yes /* Volumetric shape plot option */ 2 /* 1 - Surface info available from input, 2 - Edge info. available */ /* The next three inputs indicate the number of grid points along Psi, Eta and /* Zeta directions */ 9 9 9 /* The next three inputs are the scale factors for the x,y and z coordinates /* along Psi, Eta and Zeta directions respectively */ 0.0625 0.0625 0.0625 /* The next three inputs are */ /* 1. No of V sweeps */ /* 2. The number of pre smoothing operations for mesh generation */ /* 3. The number of post smoothing operations for mesh generation */ /* respectively */ 2 5 5 /* File into which the final results of mesh generation are stored */ Temp_new.m /* The rest of this file is the coordinate information in surface or edge */ /* See the Appendix A of the Thesis for a description of the format of */ /* this input */ HHOOOOOOOOO Homxiozoric-wioi-ao OOOOOOOOOOO 94 W 5 /* Display time every 'n' time units */ 4 /* The next three inputs indicate the */ 4 /* surfaces for which the color maps are */ 4 /* to be plotted */ yes /* Option for time Vs variable plots */ no /* Option for surface color maps */ Temp /* Indicates the variable to be plotted for */ /* Temp - Temperature */ /* Cure - Cure */ /* Both - Temperature and Cure */ /* Choosing this option also saves the generated data */ /* into matlab binary '.mat' files. */ 95 W 300 /* Final time upto which simulation needs to be performed */ 5 /* Time step :- delt */ 5 /* The last three parameters are the number of pre and post smoothing 5 /* relaxation sweeps and the number of V iterations respectively. */ 5 5 /* Time step (delt) */ 96 W y /* Type of conv coeff y => constant on each of the 6 surfaces 11 => different at each point on each of the 6 surfaces */ 0 0 0.0018 0.0018 0 0 2 / * Number of steps in thermal cure cycle*/ 0 /* Startup time */ 375.15 /* lst step end temperature */ 900 /* Second step startup time */ 450.15 /* Second step and temperature */ 4800 /* Third step startup time */ 450.15/* Third step ending temperature */ 303.15 /* Initial temperature */ 0 /* Initial cure */ 10000 /* Final time for thermal cure cycle*/ 97 WW /* Property data : This is subject to change depending on the functionalities /* used within the program for the different properties */ 4.457e-3(Thermal Conductivity) 0.940(Specific Heat) 0(Fiber volume fraction) 1.520(Density of material) /* RxN data */ 35.033e6 (klo) 8.07 e4 (E1) -33.5667e6 (k20) 7.78e4 (E2) 3266.6667 (k3o) 5.66e4 (E3) 0.60 (X _ge1) 0.47 (B) -175.20 (DELTA H) APPENDIX C Matlab Script Files mgam; init displax,m % This is an initial setup file for Matlab plots. Edit the indicated portion % of this file before choosing plot option in the COORD.DAT file. global tfinal tdisp imagp Tt Xt time Tin Xin x y z; tfinal=0; tdisp=0; % Edit the following matrix 'imagp' . Each row of the matrix indicates the % natural coordinate point of the shape. Simultaneous plots of all these % points is done during the simulation. imagp=[ 4 0 4 4 2 4 4 4 4 4 6 4 4 8 4 1; Tt=[]; Xt=D; time==[]; 'I‘in=[l; Xin=[]; x=[]; y=[]; z=[]; 98 99 Bmamafinm % This function finds the picture/figure with a given title/name. function handle = find_pic(string) f=1; g=0; while(g~=1), str1=get(f,'name'); ifistrcmp(string,str1)==1) g=1; else f=f+1; end; end; iflg==1) handle=f; else handle=gcf; end; 100 W % This function arranges the dimension of the 3D matrix to reflect different % primary directions. function [xx,yy,zz] = arrange_shape(x,y,z,n) nr=max(size(x))/n; xx=reshape(x,nr,n); yy=reshape(y,nr,n); zz=reshape(z,nr,n); 101 E I . I. l O] I % The 3D volume is visualized as a stack of surfaces. These surfaces can be % stacked along different natural coordinate directions. % The following function redistributes these stacks along a given direction % specified by the index argument. function f = distribute(x,n,m,l,index) % x - 3D volumetric matrix % n,m,l - 3D dimensions % index - 1,2 or 3 specifying one of the 3 natural coordinate directions, % psi, eta or zeta. if (index == 1) f=x; else s=size(x); ifindex == 2, count_new=1; jumpi=m*l; jumpj=m; for j=1:m, for k=1:l, for i=1:n, count_old=(j-1)*jumpj+(i-1)*jumpi+k; if min(s) > 1, if s(1) < 8(2), f(1:min(s),count_new)=x(1:min(s),count_old); else f(count_new,1:min(s))=x(count_old,1:min(s)); end; else f(count_new)=x(count_old); end; count_new=count_new+1; end; end; end; else count_new=1; jumpi=m*l; iumpj=m; for k=1:1, for i=1:n, for j=1:m, 102 count_old=(i-1)*jumpj+(i-1)*jumpi+k; if min(s) > 1, if s(1) < 3(2). f(1:min(s),count_new)=x(1:min(s),count_old); else f(count_new,1:min(s))=x(count_old,1:min(s)); end; else f(count_new)=x(count_old); end; count_new=count_new+1; end; end; end; end; end; 103 W % Plots the different surfaces of the 3D volumetric matrix. function h = surf_plot(x,y,z,i,n,c) if nargin == 5, c=abs(sqrt(x."2+y."2+z.“2)); end; s=size(x); nr=max(s); m=round(nr/n); if 8(1) > 8(2)) xg=reshape(x([ltnrl.i).n,m); yg=reshape(y([l:nr],i),n,m); zg=reshape(z([ltnrl.i),n,m); cg=reshape(<:([linrl,i).n.m); else xg=reshape(x(i,[1:nr]),n,m); yg=reshape(y(i,[1:nr]),n,m); zg=reshape(z(i,[linrl).n.m); cg=reshape(C(i,[1:nr]),n,m); end; surflxg,yg,zg,cg); 104 W % Recursively calls the surf_plot function to recreate the volume of the % complex shaped body. function h = shape_plot(x,y,z,n,m,l,index,c) if nargin == 7, c=abs(sqrt(x."2+y."2+z."2)); end; x=distribute(x,n,m,l,1); y=distribute(y,n,m,l,1); z=distribute(z,n,m,l,1); nr=round(max(size(x))/n); [xx,yy,zz]=arrange_shape(x,y,z,n); c=reshape(c,nr,n); for i=1:n, surf_plot(xx,yy,zz,i,m,-c); if i==1, hold on; end; end; x=distribute(x,n,m,l,2); y=distribute(Y»n,m,l,2); z=distribute(z,n,m,l,2); nr=round(max(size(x))/m); [xx,yy,zz]=arrange_shape(x,y,z,m); c=reshape(c,nr,m); for i=1:m, surf_plot(xx,yy,zz,i,l,-c); ifi==l, hold on; end; end; x=distribute(x,n,m,l,3); y=distribute(y,n,m,l,3); z=distribute(z,n,m,l,3); nr=round(max(size(x))ll); [xx,yy,zz]=arrange_shape(x,y,z,l); c=reshape(c,nr,l); for i=lil, surf __plot(xx,yy,zz,i,n,-c); if i== , hold on; end; end; 105 W % This is a function which saves time Vs temperature details % and plots the information for different grid points as specified in the input % files. % The saved data is in the form of a matrices containing gridpoint % information (mats) and also the time Vs variable information (Mat). function h = time_plots(t,c,n,m,l) ' global tfinal tdisp imagp Tt time Tin Xin x y 2; if t 1, for s=1:size(imagp,1). d(s)=locate(c,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); %plot(t,d(s),'o'); end; else d=locate(c,imagp(1),imagp(2),imagp(3),n,m,l); %plot(t,d,'o'); end; Tt=['I‘t;d]; time=[time;t]; plot(time,Tt); if abs(t-tfinal)<=1.0e-6, if size(imagp,1) > 1, for s=1:size(imagp,1), xd(s)=locate(x,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); yd(s)=locate(y,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); zd(s)=locate(z,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); end; - else xd=locate(x,imagp(1),imagp(2),imagp(3),n,m,l); yd=locate(y,imagp(1),imagp(2),imagp(3),n,m,l); zd=locate(z,imagp(1),imagp(2),imagp(3),n,m,l); end; mats=[xd' yd' zd']; Mat=[time Tt]; save pltT.mat mats Mat; end; 106 BMW % This is a function which saves time Vs cure details % and plots the information for different grid points as specified in the input % files. % The saved data is in the form of a matrices containing gridpoint % information (mats) and also the time Vs variable information (Matx). function h = time_plots(t,c,n,m,l) global tfinal tdisp imagp Tt Xt time Tin Xin x y 2; if t<1.0e-8, set(gca,'NextPlot','add'); end; if size(imagp,1) > 1, for s=1:size(imagp,1), d(s)=locate(c,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); %plot(t,d(s),'o'); end; else d=locate(c,imagp(1),imagp(2),imagp(3),n,m,l); %plot(t,d,'o'); end; Xt=[Xt;d]; time=[time;t]; plot(time,Xt); if abs(t-tfinal)<=1.0e-6, if size(imagp,1) > 1, for s=1:size(imagp,1), xd(s)=locate(x,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); yd(s)=locate(y,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); zd(s)=locate(z,imagp(s,1),imagp(s,2),imagp(s,3),n,m,l); end; else xd=locate(x,imagp(1),imagp(2),imagp(3),n,m,l); yd=locate(y,imagp(1),imagp(2),imagp(3),n,m,l); zd=locate(z,imagp(1),imagp(2),imagp(3),n,m,l); end; mats=[xd' yd' zd']; Mat=[time Xt]; save pltX.mat matx Matx; end; APPENDIX D C CODE FOR SIMULATION W This is the main program which calculates the mesh for a given complex body by using various input files as indicated and also passes the results of the mesh generation procedure on to the profile simulator which calculates the temperature and cure profiles. #include finclude #include #include"grid.h" #include'engineh" char gra[5]; /" Choice for graphical output ‘/ char plot_choicel 10]; /* Temperature/Cure/Both profiles ‘/ struct varray_xyz head; /"‘ XYZ grid structure for multigrid method for shape generation '/ Engine ‘eptr; /" Matlab Workspace pointer ‘/ char map151,ttx15]; /" Options for map plots and time- temperature-cure plots ‘/ int N ,M,K; /‘ No of mesh points along Psi,Eta and Zei directions respectively ‘/ int n0,nl,n2; /" Multigrid iteration parameters for shape generation 1’/ l‘ The following are different files of input and output for data transfer inp_f -- Shape and coordinates input "/ FILE ‘inp_f,"new_f,*sc_f,"bc_f,"'r1_f,‘r2_f,"‘r3_f; FILE "numal_input,“proc_inp,’prop_f,*disp_f,‘mat_lab; struct varray top; I" Head list for varray ‘/ double scale.x,scaleqy,scale_z; /* Scaling parameters for x,y and z coordinates for the shape to be modelled */ double ”’transfeLcoeff; l’ Heat transfer coefficient matrix "/ double ”thermal__cycle_data; /“ Thermal cycle data matrix ’/ double T_initial,X_initial; /" Initial temperature and cure *I double Cp_f,Cp_m,rho_f,rho_m,kf,km; /" Specific heat, density and thermal conductivity parameters for property data ‘/ double vol_f,delh; /" Volume fraction and heat of rxn respectively ‘/ double k lo,k20,k3o,E1,E2,E3,X_gel,B; /" Rxn rate equation parameters "'/ int no,nt0,ntl,nt2; /" Iteration parameters for multing method ‘/ int psi,eta,zei; /’ Display coordinates for graphical output (natural) ‘/ int bou nd_condition; /" Nature of boundary condition 1. Dirichlet 2. Mixed ‘/ 107 double final_time,t_display,delt; P Final time for simulation, time of display and step time */ main() 1 extern char gralfil; extern int N,M,K; extern double scale_x,scale_y,scale_z; struct point ’grid_pt,‘set_grid(); double *"mem_alloc_3D(); double ssox,sssy’sssz; char mat_file[20]; int i,ch; /* Open input file to read in coordinates of surfaces or edges ’/ inp_f=fopen("COORD.DAT”,”r”); fscanflinp_f,"%s",gra); fscanfiinp_f,"%d",&ch); fscanflinp_f,"%d",&N); fscanfiinp_f,"%d",&M); fscanflinp_f,"%d",&K); fscanflinp_f,"%ll”,&scale_x); fscanfiinp_f,"%lf',&scale.v); fscanflinp_f,"%lf',&scale_z); fscanflinp_f,"%d",&n0); fscanflinpj,”%d",&n1); fscanflinp_f,"%d",&n2); fscanflinp_f,"%s",mat_file); mat_lab=fopen(mat_file,"w"); /* Allocate memory for x,y,z values ‘/ x=mem_alloc_3 [)(N,M,K); y=mem__alloc_3 D(N,M,K); z=mem_alloc_3D(N,M,K); read_input(Ch.X.y,z); fclose(inp_0; P bc_f=fopen("check.m","w"); print_input(bc_f,x,y,z,N,M,K); fclosemc.0; ‘/ /" If display is required open matlab workspace ‘/ iflstrcmp(gra,”yes")==0) ' eptr=eng0pen(); surfaee_approx(ch,x,l‘1,M,K); surface_approx(ch,y,N,M,K); surface_,approx(ch,z,N,M,K); I‘ Set appropriate structures for multigrd method for mesh generation */ Struct_set(x,y,z); l‘ new_f=fopen("0ut.m”,"w"); print_out(new_f,head.x,head.y,head.z,N,M,K); fclose(new_f); */ iflstrcmp(gra,"yes")==0) [ plot_surfl l,x,y,z,(doub1e "*)NULL,N,M,K); vie w_options(x,y,z,(double “"‘)NULL,N,M,K, l ); print_plot_option8(); l FAS_xyz(&head); sc_f=fopen("O.m","w"); print_out(sc_f,head.x,head.y,head.z,N,M,K); felose(sc_l); iflstrcmp(gra,"yes")==0) { plot_surfi2,x,y,z,(double ‘**)NULL,N,M,K); view_options(x,y,z,(double ***)NULL,N,M,K,2); print_plot_optio ns( ); I /"‘ Set structures to pass shape information to profile generation program */ grid_pt=set_grid(N,M,K); update _grid(grid_pt,x,y,z,N,M,K); free__list(); printf("Mesh generation complete. Beginning thermal simulation. \n"); /"‘ Call profile simulation */ profile _gen(grid_pt,N,M,K); fclose(mat_lab); I" If Matlab workspace is open close it and shut down 1’/ ifistrcmp(gra,"yes")==0) engClose(eptr); ) E . ll . . I This program contains routines for all the different operators for the multigrid method as well as the relaxation procedure for three dimensional elliptic generation difference equations with dirichlet boundary conditions for mesh generation. #include #include #include #include "grid.h" f" Restriction operator for multigrid algorithm */ Restrict(fp,f,n,m,l) double *"fp,**"f; /" fp : coarse grid matrix of dependent variabes f : fine grid matrix of dependent variables "/ int n,m,l; /" Mesh size */ 1 double term1,term2,term3,term4; int pas; int iik; p=(n+l)/2; q=(m+1)/2; r=(l+1)/2; for(i=0;in>3)&&(ptr->m>3)&&(ptr->l>3)) l temp=(struct varray_xyz *lmallodsizeoflstruct varray_xyz»; temp->ndintXptr->n+1)/2; temp->m4intXptr->m+1)/2; temp->l=(int)(ptr->l+l)/2; temp->x=mem_alloc_3D(temp->n,temp->m,temp->l); temp->y=mem_alloc_3D(tempo>n,temp->m,temp->l); temp->z=mem_alloc_3D(temp->n,temp->m,temp->l); temp->Sx=mem_alloc_3D(temp->n,temp->m,temp— >1); temp->Sy=mem_alloc_3D(temp->n,temp->m,temp- >1); temp->Sz=mem_alloc_3 D(temp->n,temp->m,temp- >1); temp->next=NULL; ptr->next=temp; ptr=ptr~>next; 1 return; 1 I“ Fully multigrid V algorithm for mesh generation */ l“ Details of the routine are in FMV algorithm of the program temp _grid.c */ FMV_xyz(v) struct varray_xyz *v; { extern int nl,n2; int iJJt; struct varray_xyz *temp; double sssfx’sssfy,ssafz,sssflx,sssf1y,sssflz’sssmx’sssmy,s surf-22; double "”rx,""ry,*"rz,trunc,err,norm_xyz(); double ssstx’sssty’ssstz,sssvx’sasvy’sssvz; int vn,vm,vl,tn,tm,tl; double "*mem_alloc_3D(); P printf("FMV_xyz begins here\n"); */ vn=v~>n; vm=v->m; vl=v->l; =V->X; VY=V->Y; vz=v->z; ifiv->next==NULL) l /" printf("FMV_xyz : (‘.oarsest Grid \n"); */ Relax_xyz(vx,vy,vz,v->Sx,v—>Sy,v~ >Sz,vn,vm,vl,n1+n2); 1 else 1 Relax_xyz(vx,vy,vz,v->Sx,v->Sy,v->Sz,vn,vm,vl,n1); temp=v->next; tn=temp~>n; tm=temp->m; tl=temp->1; tx=temp~>x; ty=temp->y; tz=temp->z; f1x=me m_alloc_31)(tn,tm,tl); f1y=mem_alloc_3 D(tn,tm,tl); flz=mem_alloc_3 D(tn,tm,tl); fx=mem_alloc_3D(tn,tm,tl); =mem_alloc_3 D(tn,tm,tl ); f z=mem_alloc_3 D(tn,tm,tl ); rx=mem_alloc_3D(tn,tm,tl); ry=mem_alloc_3 D(tn,tm,tl); rz=menLalloc_3 D(tn,tm,tl); =mem_alloc_3 D(vn,vm,vl); f2y=mem_alloc_3 D(vn,vm,vl); i‘Zz=me m_alloc_3 D(vn,vm,v1); Restrict(rx,vx,tn ,tm,tl ); Restrict(ry,vy,tn,tm,tl ); Restrict(rz,vz,tn,tm,tl); A_xyz(fit,fy,fz,rx,ry,rz,tn,tm,tl ); A_xyz(l‘2x,f2y,f2 z,vx,vy,vz,vn,vm,vl ); Restrict(f1x,f2x,tn,tm,tl); Restrict(f1y,l‘2y,tn,tm,tl); Restrict(fl z,l‘22,tn,tm,tl); mat__sub(fx,fx,flx,tn,tm,tl); mat_sub(fy,fy,f1y,tn,tm,tl ); mat_sub(fz,fz,f1 z,tn,tm,tl ); v->tru nc= norm_xyz(fit,fy,fz,tn,tm,tl ); Restrict(tx,vx,tn ,tm,tl); Restrict(ty,vy,tn ,tm,tl); Restrict(tz,vz,tn,tm,tl); mat_copy(rx,tnt,tn,tm,tl); mat_copy(ry,ty,tn,tm,tl ); mat_coPy(rz,tz,tn,tm,tl ); A _xyz(flx,f1y,f1 z,tx,ty,tz,tn,tl,tn); A _xyz(f2x,f2y,f2 z,vx,vy,vz,vn,vm,vl); mat__sub(f2x,f2x,v->Sx,vn,vm,vl); mat_sub(f2y,f‘2y,v->Sy,vn,vm,vl); mat_sub(f22,f2z,v->Sy,vn,vm,vl); Restrict( fx,f2x,tn,tm,tl ); Restrict( fy,1‘2y,tn,tm,tl); Restrict(fz,122,tn,tm,tl); mat_sub(temp->Sx,flx,fx,tn,tm,tl); mat_sub(temp->Sy,f1y,fy,tn,tm,tl ); mat_sub(temp->Sz,f1 z,fz,tn,trn,tl ); for(i=0;i<=1;i++) FMV_xyz(te mp); mat_sub(f1x,tx,rx,tn,tm,tl ); mat_sub(f1y,ty,ry,tn,tm,tl); mat_sub(f1 z,tz,rz,tn,tm,tl); lnterpolate(i‘2x,f1x,tn,tm,tl); Interpolate(f2y,fly,tn,tm,tl); Interpolate(l22,flz,tn,tm,tl); mat_add(vx,l‘2x,vx,vn,vm, v1); mat_add(vy,f2y,vy,vn,vm, v1); mat_add(vz,l‘2z,vz,vn,vm,vl); Relax_xyz(vx,vy,vz,v->Sx,v->Sy,v->Sz,vn,vm,vl,n2 ); free_3 D(rx,tn,tm,tl); free_3 I)(ry,tn,tm,tl ); free_3 D(rz,tn,tm,tl); free_3 D(fx,tn,tm,tl); free_3D(fy,tn,tm,tl); free_3 D(fz,tn,tm,tl); free_3 D(flx,tn,tm,tl); free_3D(f1 y,tn,tm,tl); free_30(f1 z,tn,tm,tl); free_3 D(i‘Zx,vn,vm,vl); free_3 D(l2y,vn,vm,vl); free_3 D(f‘Zz,vn,vm,vl); 110 ) l’ printf("FMV_xyz (N = %d) ends here\n",vn); "/ return; i /" FAS algorithm for mesh generation */ /’ Details of the different steps are given in FAS routine in the program 'temp_grid.c' */ FAS_xyz(v) struct varray_xyz *v; I extern int n0; int i,j,k,cnt; struct varray_xyz "temp; double g; double essk’sssfy,stsfz’ssaflx’sssny’sssflz,sssmx,sssf2y,s "[22; double "*rx,"’ry,"‘rz,err,trunc,norm_xyz(l; double “*mem_alloc_3D(); double asstx,ssaty’ssstz,ssavx’assvy,tssvz; int vn,vm,vl,tn,tm,tl; vn=v->n; vm=v->m; vl=v->l; vx=v->x; vy=v->y; vz=v->z; if(v->next==NULL) l /* printf("FAS_xyz : Coarsest Grid \n"); */ for(i=1;i<=n0;i++) FMV_xyz(v); ) else i l‘ printfl"FAS_xyz (N = %d) begins here \n",vn); */ temp=v->next; tn=temp->n; tm=temp->m; tl=temp->l; tx=temp->x; ty=temp->y; tz=temp->z; Restrict(tx,vx,tn,tm,tl); Restrict(ty,vy,tn,tm,tl); Restrict(tz,vz,tn,tm,tl); f1x=mem_alloc_30(tn,tm,tl); fly=mem_alloc_3D(tn,tm,tl); flz=mem_alloc_3D(tn,tm,tll; fx=mem_alloc_3D(tn,tm,tl); fy=mem_alloc_3D(tn,tm,tl); fz=mem_alloc_3D(tn,tm,tl); rx=mem_alloc_3D(tn,tm,tl); ry=mem_alloc_3D(tn,tm,tl); rz=mem_alloc_30(tn,tm,tl); l2x=mem_alloc_3 D(vn,vm,vl); f2y=mem_alloc_3D(vn,vm,vl); f22=mem_alloc_3D(vn,vm,vl); mat_copy(rx,tx,tn,tm,tl); mat,copy(ry,ty,tn,tm,t1); mat_copy(rz,tz,tn,tm,tl); A _xyz(f1x,fly,f1z,tx,ty,tz,tn,tm,tl); A_xyz(l2x,l2y,f2z,vx,vy,vz,vn,vl,vn); mat_sub(f‘Zx,I‘2x,v->Sx,vn,vm,vl); mat_sub(t‘Zy,f‘2y,v->Sy,vn,vm,vl); mat_sub(l‘22J2z,v->Sz,vn,vm,vl); Restrict(fx,12x,tn,tm,tl); Restrict(fy,f2y,tn,tm,tl); Restrict(fz,l22,tn,tm,tl); mat_sub(temp->Sx,flx,fit,tn,tm,tll; mat_sub(temp—>Sy,f1y,fy,tn,tm,tl); mat_sub(temp->Sz,f1z,fz,tn,tm,tl ); FAS_xyz(temp); mat_sub(flx,tx,rx,tn,tm,tl); mat_sub(fly,ty,ry,tn,tm,tl); mat_sub(flz,tz,rz,tn,tm,tl); Interpolate(f2x,f1x,tn,tm,tl); lnterpolate(l2y,f1y,tn,tm,tl); Interpolate(l22,flz,tn,tm,tl); mat_add(vx,i‘2x,vx,vn,vm,vl); mat_add(vy,i2y,vy,vn,vm,vl); mat_add(vz,f2 z,vz,vn,vm, v1); /*if(vn==9) print_out(vx,vy,vz,vn,vm,vl); */ err: l; cnt= l; s=0; while(cnt<=n0)f" &&(err>=0.333*v- >tru nc)&&(fabs((g-err/v.>trunc )/g)> 1.09-4))”/ 1 111011“: 1 ) g=err/v->trunc; else s=1; FMV_xyz(v); A_xyz(f2x,f2y,l2z,vx,vy,vz,vn,vm,v1); mat_sub(f‘2x,v->Sx,i‘2x,vn,vm,vl); mat_sub(l‘Zy,v->Sy,f‘2y,vn,vm,vl); mat_sub(f‘Zz,v->Sz,12z,vn,vm,vl); err=norm_xyz(12x,l2y,f‘2z,vn,vm,vl); /* printf("FAS_xyz (N = %d): Error (Old) : %6.4lf Error (New) : %6.4lf \n",vn,g,err/v->trunc); */ cnt++; ) free_3 D(rx,tn,tm,tl); free_3 D(ry,tn,tm,t1); free_3D(rz,tn,tm,tl); free_30(fx,tn,tm,tl); free_3IXfy,tn,tm,tl); free_3 D( f z,tn,tm,tl); free_3D(flx,tn,tm,tl); free_3D(f1y,tn,tm,tl); free_30(fl z,tn,tm,tl); free_30(f2x,vn,vm,vl); free_3D(f2y,vn,vm,vl); free_3D(f22,vn,vm,vl); ) /" printfl"FAS_xyz (N = %d) ends here \n",vn); ’/ return; I Was This program calculates the first approximation to the shape of the body using transfinite interpolation in two and three dimensions. #inclu de lll #include #include #include "grid.h" /* Transfinite Interpolation Routine */ transfinite_3D(f,n,m,l) double "*f; int n,m,l; 1 double “"Fl,*“mem_alloc_3 I)(),int_plt(); double "*F2; double ""F3; int i,j,k; Fl=mem_alloc_3D(l,n,m); /* Correction along j(eta) direction */ for(k=0;k #include #include #include "grid.h" int cyclic(i) P Routine which returns the cyclic value ofi in (1,2,3) ’/ int i; 1 int cyl; switch(i) 1 case 0: l cyl=l; break; return(cyl); ) set_matrix(X,Y,Z,n,m,l,iJ,k,mat) /* Sets the jacobian matrix at any grid point in the transformed coordinate */ double sssxsesy’sssz; int n.m,1,iJ.k; double ”mat; ( double f_psi(),f_eta(),f_zei(); double *"temp_ptr,sum; int 9.9; for(p=0;p<3;p++) l switch(p) 1 case 0: l temp_ptr=X; break; 1 case 1: l temp _ptr=Y; break; 1 case 2: { temp_ptr=Z; break; l l for(q=0;q<3;q++) 1 switch(q) 1 case 0: 1 matl p ][q ]=f_psi(temp_ptr,iJ,k,n.m.l, I ); break; ) case 1: 1 matl p][q l=f_eta(temp_ptr,i,j,k,n,m,l, l 1; break; 1 case 2: 1 mat[p][q]=f_zei(temp_ptr,ij,k,n,m,l,1); break; 1 l l l /* for(p=0;p<3;p++) 1 sum=0; for(q=0;q<3;q++) sum=sum+pow(fabs(mat[p][q]),2.0); for(q=0;q<3;q++) matlpllql=matlpllqllsqrflsumk l "/ return; i calc_beta(Mat,b) /"‘ Calculate the beta matrix using 'M' matrix values at a given gridpoint */ double ”Mat,"b; /* Mat - Matrix 'M' (for the transformation) , b - Returned matrix of betas ’/ 1 int i,j,k,l,m,n; for(i=0;i<3;i++) 1 k=cyclic(i); m=cyclic(k); for(i=0;i<3;i++) 1 l=cyclic(j); n=cyclic(l); blileI=Matlk][ll*Mat[m][n]-Mat[k][n]*Mat[m][l]; l I return; 1 double calc_alpha_jac(a,X,Y,Z,n,m,l,iJ,k) /"‘ Calculation of alpha values using beta matrix */ double "a,***X,"**Y,"*Z;/* X,Y,Z - x,y,z arrays; i,j,k - location; a - returned matrix of alphas "/ int n,m,l,i,j,k; 1 double sum,”mat,"b; double jac,Jacohian(); 114 int oipiq; mat=(double **)malloc(3*sizeof(double *)); b=(double ")malloc(3"‘sizeof(double ‘1); for(o=0;o<3;o++) 1 mat[o ]=(double *)malloc(3*sizeofldouble)); blo]=(double ‘)malloc(3"'sizeof(double)); ) set_matrix(X,Y,Z,n,m,l,i,j,k,mat); calc_beta(mat,b); jac=lacobian(mat); for(o=0;o<3;o++) 1 for(p=0;p<3;p++) 1 sum=0.0; for(q=0;q<3;q++) surn=sum+b1q1101"b1qllpl; 81011p1=811mi ) ) free_2 D( mat,3,3 ); free_21)(b,3,3); return(jac); 1 double Jacobian(mat) /“ Calculates J acobian(J ) given the transformation matrix (M) J: l M l */ double ”mat; /‘ mat - Transformation matrix (M) */ 1 double ret; ret=mat10][ 0]*(mat[ 1][ 11*mat12][2]- matl 11121‘m3t12ll 1))-mathll 11*(mat111101‘mat12ll2l- matl 1l[2]*mat[2][0])+mat[0]l2]"‘(matl l][0]"‘mat[2][ 1]- matl 111 l]"‘mat[2][0]); return(ret); } /‘ Setting Boundary Conditions Over The Volume ‘/ set_boundary(fx,fy,fz,fxl,fy1,le,n,m,l) double team’sssfy,sssfz; double "*fxl,""'fyl,""le; int n,m,1; 1 int step_i,stepJ,step__k,i,j,k,p; for(p=1;p<=3;p++) 1 switch(p) 1 P Along Psi Direction */ case 1: { step_i=1; stePJ=1; step_k=l- 1; break; 1 /* Along Eta Direction */ case 2: l step_i=n-1; stePJ=1; step_k=l; break; 1 /* Along Zei Direction */ case 3: { step_i=1; stspasm-1; step_k=1; break; 1 ) for(i=0;ix=mem_alloc_3 D(n,m,l); grid_pt->y=mem_alloc_3 D(n,m,l); grid_pt->z=mem_alloc_3D(n,m,l); grid_pt->J=mem_alloc_3D(n,m,l); grid_pt->a=(struct D2_array ’**)malloc(n*sizeof(struct DZ_array “1); for(i=0;ia[i]=(struct D2_array ")malloc(m"‘sizeoflstruct [)2_array *)); for(i=0;ia[i]Lj]=(struct D2_array *)malloc(l"‘sizeof(struct D2_arrayl); for(k=0;ka[i]Li][k].c=(double **)malloc(3*sizeofldouble ")); for(p=0;p<3;p++) grid_pt->a[i1Li 11k ].c[ p1=(double ‘)malloc(3"‘sizeof(double)); ) ) l return(grid_pt); ) /“ Updation of grid values given the coordinate values at each grid p0 int and the mesh size */ update _grid(grid_pt,fx,fy,fz,n,m,l) struct point *grid_pt; double sssfi’sssfy,sssfz; int n,m,1; 1 int iJ,k,p,q; double *‘an,norm(),jn,**A_norm(); mat_copy(grid_pt->x,fir,n,m,1); mat_copy(grid_pt->y,fy.n.m,1); mat_copy(grid_pt->z,fz,n,m,l); for(i=0;iJ[ilLillkl=calc_alphajac(grid_pt- >a[i lLi llk].c,fit,fy,fz,n,m,l,iJ,k); l ) return; 1 double "A_norm(grid_pt,n,m,l) struct point *grid_pt; int n,m,1; 1 int iJ,k.p.q; double *"temp,"a,”*mem_alloc_3 D(),**mem_alloc_2 DO; a=mem_alloc_2D(3,3); temp=meuLalloc_3D(n,m,l); for(p=0;p<3;p++) 1 for(q=p;q<3;q++) 1 for(i=0;ia[ i ][j 11k].c1 p 11 q 1; l l a] p ][q ]=norm(temp,n,m,l); 1 ) free_3D(temp,n,m,l); return(a); 1 /" The mesh generation operator evaluation for each of the three ellip tic equations involved in the generation I"/ A _xyz(outx,outy,outz,ux,uy,uz,n,m,l) double *"outx,"‘**outy,*"outz;/* Output matrices corresponding to the three equations */ double "‘ux,"“uy,"*uz; /"‘ Input values for the three coordinates x, y, 2 as a function of grid points in the natural coordinates ‘/ int n,m,1; /* Mesh size (required for multigrid operation) ‘/ 1 int iik; double brl,fx2,fx3,fy1,fy2,fy3,fz1,fz2,f23,jac,P1),Q(),R(),**al; double f_psi(),f_eta(),f_zei(),f_psi_psi(),f_zei_zei(),f_eta_cta(),f _psi_eta(); double f_eta_zei().f_zei_psi(),calc_alpha_jac(),**mcm_alloc_2 DO; set_boundary(ux,uy,uz,outx,outy,outz,n,m,1); al=mem_alloc_2D(3,3); for(i=1;itol)) 1 if((n==3)l I (m==3) l I (1==3)) 1 mat_copy(tempx_err,ux,n,m,l); mat_copy(tempy_err,uy,n,m,l); mat_copy(tempz_err,uz,n,m,l); error=norm_xyz(ux,uy,uz,n,m,l); } cnti=1; begi=2; /* Red-Black Sweeps ‘/ while(cnti<=2) 1 /" Psi direction */ for(i=begi;i #include #include "grid.h" #include "engineh" /‘ Profile generation program to evaluate temperature and cure profiles as a function of time */ profile _gen(grid_pt,n,m,l) struct point ‘grid_pt; f" Shape information for every grid point on finest grid */ int n,m,1; /" Mesh size over finest grid */ 1 extern int psi,eta,zei,bound_condition; int i,bc; extern char gra15]; double ‘*"‘T,”*X,td,"“mem_alloc_3 11(1; extern double final_time,de1t,t_display; extern struct varray top; extern Engine ‘eptr; extern int N,M,K; extern FILE ‘numal_input,"proc_inp,“prop_f,"disp_f; double t,find_min(),find_max(),MinT,MaxT; /" Read in property, process and simulation parameters for multigrid method */ proc_inp=fopen("PR()CESS.[)AT',"r"); prop_f=fopen("PROPERTY.[)AT","r"); numal_input=fopen("NUM . DAT","r"); read_process_input(grid_pt,n,m,11; f" Read in display options for profile simulation program ‘/ disp__ =fopen("DISPLAY.DAT","r"); read_profile_plot_options(n,m,l); iflstrcmp(gra,"yes"1==0) initialize_workspace(grid_pt,N,M,K); T=mem_alloc_31)(N,M.K); X=mem_alloc_3D(N,M,K); t=0; td=t_display; /"‘ Set initial conditions */ 1nitial_Condition(T,X,N,M,K); printf("Time : %6.31f; ",t); MinT=find_min(T,N,M,K); MaxT=find_max(T,N,M,K); printf("MinT : %6.4lf; Max’l‘ : %6.41f \n",MinT,MaxT); if(strcmp(gra,"yes")==0) Display_variables(t,grid_pt,T,X,N,M,K); bc=bound_condition; /" Begin simulating over time until final time "I while( f abs(t-final_time)> 1.0e-6) 1 I‘ Set and allocate appropriate structures for multigrid method *‘/ arrange_struct(t,grid_pt,T,X,n,m,l,bc); /" Perform FAS algorithm at current time ‘/ FAS(t,&top,bc); /"' Increment time "/ =t+delt; /" If it is time for display, display the results and show graphically if required ’/ if(fabs(t-td)<= l .0e-6) 1 printf("Time : %6.31f; ”,t); MinT=find_min(top.T,N,M,K); MaxT=find_max(top.T,N,M,K); printf("Min'I‘ : %6.4lf; MaxT : %6.4lf \n",MinT,MaxT); iflstrcmp(gra,"yes")==01 Display_variables(t,grid_pt,top.T,top.X,N.M,K); td=td+t_display; 1 118 /"' Update next guess as the current temperature and cure profiles "/ mat_cop)’(T.top.T,n,m,l); mat_copy(X,top.X,n,m,l); 1 /* Print final values at all grid points into file for fututre use */ print_values(grid_pt,top.T,top.X,N,M,K); print_temp(grid_pt,T,X,n,m,l); mat_lab_options(); print_plot_options( 1; free_3I)(T,N,M,K); free_3D(X,N,M,K); return; 1 print_temp(grid_pt,T,X,n,m,l) struct point ‘grid_pt; double "*T,"*X; int n,m,1; 1 int i,j,k; double r; for(i=0;ixl n- 1 1U 1U 1.2.0)+pow(grid_pt- >y[ n- 1 NJ lljl,2.0)+pow(grid_pt->z1n- l 11) 11,11,201); printf("R = %6.5lf; T[%d %d %d] = %6.4lf; X[%d %d %d] : %7.51f \n",r,n-1,j,k,T[n-l][j][k],n-1,j,k,X[n— 11111110); 1 1 return; 1 print_shape_temp(grid_pt,T,X,n,m,l) struct point "grid_pt; double *"T,*"X; int n,m,1; 1 double "‘x,”"y,"‘"z,”a,jac; extern FILE I"mat_lab; int i,j,k,lirn,s; double r; x=grid_pt->x; y=srid_Pt->y; z=grid_pt->z; mat_lab=fopen("TempV.m","w"); for(i=0;i<3;i++) 1 for(j=0;j<3;j++) 1 for(k=0;k<3;k++) 1 r=sqrt(pow(x[i][j][k],2.0)+pow(ylillj11kl,2.0)+pow(z[i1L1] [k],2.0)); fprintflmat_lab,"(%d,%d,%d); x:%6.41f y:%6.41f z:%6.4lf r:%6.41f T:%6.31f X:%6.3ll\n",i,j,k,x1111.i11k1,y1111111kl,zlilljllkl,r,Tlilljllk 1. XML) 111(1); 1 1 1 for(i=0;i<3;i++) 1 for(i=0;i<3;j++) 1 for(k=0;k<3;k++) 1 a=grid_pt->a[ i ][j ][k ].c; jac=grid_pt»>J[il[j11kl; fprintflmat_1ab,"(%d,%d,%d);a1 l:%8.7lf; a 12:%8.7 If ; a13:%8.71f; a22:%8.71f; a23:%8.7lf; a33:%8.71f; J:%8.7lf\n",i.i,k,a[01[0],a[O][ 11310112131 111 11,al 11121.al 2112lJac); 1 1 1 for(i=8;i>5;i--) 1 for(i=8;j>5;j--) 1 for(k=8;k>5;k--) 1 r=sqrt(pow(x[i]1j][k],2.0)+pow(y[i]le[k],2.0)+pow(z[i]le [k],2.0)); fprintfimat_lab,"(%d,%d,%d); x:%6.41f y:%6.4lf z:%6.4lf r:%6.4lf T:%6.31f X:%6.3lf\n",i,j,k,x[i ][j 11k1,y111Li 11k l,z[illj 11k1,T.'1‘1i1Li 11k 1, Xlillj 111(1); 1 1 1 for(i=8;i>5;i--) 1 for(i=8;j>5;j--) 1 for(k=8;k>5;k--) 1 a=grid_pt->a[ i ][j ][k ].c; jac=grid_pt—>J111U11k1; fprintflmat_lab,"(%d,%d,%d);al 1:%8.71f; a 12:%8.7 1f ; a13:%8.71f; a22:%8.71f; a23:%8.71f; a33:%8.71f; J!%8-711\D".iilik,8101101,31011 11.8101121,81111 11,81 11121,a1 21121Jac); 1 1 1 return; 1 print_va1ues(grid_pt,T,X,n,m,l) struct point *grid_pt; double *"T,*"X; int n,m,1; 1 double tttx,*t*y,t**z; double "a; extern FILE *mat_lab; int iik; double r,fmd_min(),find_max(); x=grid_pt->x; y=srid_pt->y; 119 z=grid_pt->z; fprintflmat_lab,"tmp =[ "1; for(i=0;ia[ i ][j ][k ].c; fprintflmatJab,"%d %d %d a1 l:%8.71f a22:%8.7lf a33:%8.71fa12:%8.71fa13:%8.71fa23:%8.71f J1%8-711\n".i.1,k,8101101,81 111 11.3.1 21121.81011 1],a[0][21,a[ 11121.81'id_ptp>JIilLi11kl); 1 1 1"/ return; 1 This part of the code calculates the interior and boundary difference operators of the process differential equations and also relaxes the set of equations at the different grid levels using the FAS algorithm. #inclu de #include #include #include "grid.h" #define gas_const 8.314 /" Calculate L_h(u_h) given u_h */ Construct_L(t,grid_pt,conv,Opr,mul,T,X,n,m,l,bc1 double t; /* Time */ struct point *grid_pt; /* Shape factors for current grid 'h' ‘I double I"""‘conv; /“ Transfer coefficient matrix for current 11' "/ double "*Opr,"*mul,***T,"*X; /" Opr : Return operator mu] : Time derivative coefficient T,X : Dependnet variables (u_h) */ int n,m,l,bc; f" Mesh sizes and natuer of BC */ 1 double *“Th,*"mem_alloc_3D(1,L_interior(1; double boundary__operator(1; int i,j,k,bound(),m1,m2,m3; /" Allocate memory and calculate space derivative coefficients */ Th=mem_alloc_3D(n,m,l); fill_conductivity('Ih,T,X,n,m,l1; for(i=0;i01&&(k<1-11)) 1 ml=bound(i,n); m2=bound(i,m1; Opr[ilU11k1=T1i11j11k]-1.0/2.0*(T[i-m1 1U11kl+T1ilU- m211kl); 120 1 ifl((i==0)l l(i==n-l))&&((1>01&&(j01&&(i0)&&(i01&&(k<1- 1 )1) Opr[i ][j ]1 k ]=boundary_operator(t, l,grid_pt,T,X,Th,con v,mul.iJ,k.n.m.l); if(((i>01&&(i01&&(k01&&(i01&&(ja[i1[j][k1.c; jac=temp->J111U11k1; den=density(T[ilLi11k1,X111U11k1); /* Time derivative contribution "'/ terml=mulli1L111kl"'1‘11]Li 11k]; /"‘ Space derivative contributions */ term2=delt*a[ 0][01*f_psi_psi(Th,T,i,j,k,n1/pow(iac,2.0) rm3=delt*a[ 1 11 1l‘f__eta_eta(Th,T,i,j,k,m1/pow(jac,2.0 1; term4=delt*a[ 2 ][ 2 ]*f_zei_zei(Th,T,iJ,k,l)/pow(iac,2.0); P Cross derivatives */ term5 =delt‘a1 0 ][ 11*(f_psi__eta(Th,T,i,j,k,n,m1+f_eta__p si(Th,T,i,j,k,n,m11/pow(jac,2.0); term6=de1t*a[ l][2]*(f_eta_zei(Th,T,i,j,k,m,l)+f_zei_eta (Th,T,iJ,k,m,l))/pow(jac,2.01; term7=delt*al 01121*(f_zei_psi(Th,T,iJ,k,1,n1+f_psi_zei( Th,T,iJ,k,1,n)Vpow(jac,2.0); l" Shape-Source term contributuions "‘l term8=delt*kav*(P(i,j,k1*f_psi(T,i,j,k,n,m,1,21+Q(i,j,k1 *f_eta(T,iJ,k,n,m,1,2)+R(i,j,k1’f_zei(T,i,j,k,n,m,l,2)1; l“ Include following terms if using an implicit method term9=delt*(- delh)‘kinetics_rate('l‘[i][j 11k],X[i11j11k]1"'(den- rho_f*vol_f1; term0=delt*MW_power('1‘1i]U11k1,XIilLi11k1); "/ /* Total is sum of all contributions */ ret=term1- 0.50‘(term2+term3+term4+term5+term6+term7+ter m8); return(ret); 1 P Calculate off-diagonal neighbour contribution to relaxation source terms ‘/ double L_off_diag(grid_pt,mul,’1h,T,X,i,j,k,n,m,l) struct point *grid_pt; /"‘ Shape structure for current grid size ‘I double "*mu1,*"Th,""T,*”X; /* Time and space derivative coefficients and dependent variable matrices */ int i,j,k,n,m,l; /* Grid point coordibates(i,j,k1 and mesh sizes (n,m,l) */ 1 extern int N,M,K; struct point 'temp; extern double delh,vol_f,rho_f,delt; 121 double ret,term2,term3,term4,term5,term6,term7,term8,ter m9,term0,den,h,g,r; double "a,jac,kav,kinetics_rate(); double density(),f_psi(),f_eta(),f_zei(),f_zei_eta(1,f_eta_zei(),f_ psi_psi(),f_eta_eta(); double f_zei_zei(),f_psi_eta(1,f__eta_psi(),f_zei_psi(),f_psi_zei(1, MW_power(),kinetics_rate(),P(),Q(),R(1; /* Calculate current grid spacings ‘/ h=(double)(N- 11/(n-1); g=(doub1e)(M-11/(m-l); r=(double)(K- 11/(1-11; temp=grid_pt; /" Calculate isotropic average of space derivative coefficients */ kav=l.0/6.0*(Thli+l 11111k1+Th1i- lllj11kl+'1‘hlillj+111k1+'lh[illj- 111kl+ThIilLi11k+l1+ThlilL111k-1l1; a=temp->a[i]1j ][k].c; jac=temp->J[ilL111k1; /"‘ Neighbours along Psi direction */ term2=delt*0.5"‘('l'hli+lle11k1+Th1i- 11U11k11‘8101101*('1‘1i+11U11k]+'1‘1i- lllj 11k1V(p0W(h,2.01"‘pow(iac,2.01); /"' Neighbours along Eta direction ‘7 term3=delt*0.5‘(Thli]Lj+l 11k1+Th[i][j- 111k])*a[11111‘('1‘1i11.1+111kl+'1‘111U-1lllePOWUaC*S.2.0); /* Neighbours along Zei direction ‘/ term4=delt*0.5"(Th[ilL111k+l1+Th1ilUlIk- l11*812112l*('1‘111L111k+l1+Tlillj11k-111/pow(jac*r,2.0); /" Cross Derivative terms do not involve current grid point so they are entirely off diagonal */ term5=delt*a[01[ 1 ]*(f_psi_eta(Th,T,iJ,k,n,m1+f_eta_p si('Ih,T,i,j,k,n,m1)/pow(jac,2.01; term6=delt"a[ 1][2]*(f_eta_zei(Th,T,iJ,k,m,l)+f_zei_eta (Th.T.iJ.k.m.l))/p0W(iaC.2-0); term7=delt"a[ 0112 1*(f_zei_psi('1h,T,iJ,k,l,n)+f_psi__zei( Th,T,iJ,k,l,n))/pow(jac,2.0); /"‘ Shape factor contributions (currently zero) ’/ term8=delt*kav‘(P(iJ,k1*f_psi(T,i,j,k,n,m,l,21+Q(i,i,k) *f_eta(T,iJ,k,n,m,1,2)+R(i,j,k)*f_zei(T,i,j,k,n,m,l,2)1; /* Include following terms if using implicit method for nonlinear terms term9=delt“(- de1h1*kinetics_rate(T[ i ][j ][k],X[ i ][j ]1k]1"‘(de n- rho_f“vol_f); term0=de1t*MW_power('I‘li1U11k1,X1i11J11k]); */ l‘ Total is sum weighted by Crank-Nicolson average weight */ rat:- 0.50*(term2+term3+term4+term5+term6+term7+ter m8); return(ret); 1 /* Calculate the source terms independent of current dependent variable values */ set_rhs(t,grid_pt,conv,rhs,mul,Tg,Xg,n,m,1,bc) double t; /" Current Time ‘/ struct point ‘grid_pt; /“ Shape structure ’/ double "‘"conv,"*rhs,"*mul,***Tg,”*Xg; /* Coefficient and Prev. time step dependent variable matrices */ int n,m,l,bc; /‘ (n,m,l) : Mesh sizes; bc : Nature of BC */ 1 extern int N ,M ,K; extern double delt,delh,vol_f,rho_f; int c,i,j,k; double ht,cp,thermal_cyc1e(1,density(1,kinetics_rate(1,MW_po wer(1,h,’1‘s,rho; double coeff,source,off,"a,jac,t1,t2,t3,t4,t5,t6,t7; double density(),f_psi(),f_eta(1,f_zei(),f_zei_eta(),f_eta_zei(1,f_ psi_psi(1,f_eta_eta(); double f_zei_zei( ),f_psi_eta(),f_eta_psi( ),f_zei_psi( 1, f_psi_zei( 1; double P(),Q(),R(),heat_transfer_coeff(),kav,rxn__mw; double "*Th,"*mem_alloc_3D(1,boundary_operator(1,bou rid ary_const_temp(),thermal_conductivity(1; double compute_boundary_off_diagonal(1,combined_boundar y_source(),compute__boundary_coefficient(1; /* Calculate space derivative coefficient values at previous time ’/ Th=mem_alloc_3D(n,m,l); fill_conductivity('Ih,Tg,Xg,n,m,l); l‘ Calculate ambient temperature from cure cycle data ‘/ Ts=thermal_cycle(t+delt1; /" Over all the grid points */ for(i=0;iO1&&(j0)&&(i01l l(k<1-l))) rhs[i11j11k]=0.00; 122 /‘ For pure surface points (non-corner and non-edge) */ /’ Psi constamnt surface "/ ifl((i==01 l l(i==n-1))&&((j>0)&&(j 01&&(k01&&(i0)&&(k<1-111) 1 rho=density(Tg1U1j11kLXglilUl1kl); rxn_mw=delt*((rho- vol_f"rho_f1*kinetics_rate(Tg1i]Lj][k],Xg[i ][j 11k]1*(- delh)+MW_power(Tglillj11k1,Xgli11j11k111; c=2; source=2.0‘combined_boundary_source(t,c,grid_pt,Tg ,Xg,’1h,conv,mul,i,j,k,n,m,l1; off=compute_boundary_off_diagonal(t,c,grid_pt,Tg,Xg ,Th,conv,mul,i,j,k,n,m,l); coeff=mu11ilLi11kl- compute_boundary_coeffrcient(t,c,grid_pt,Tg,Xg,'I‘h,co nv,mul,i,j,k,n,m,l); rhs[i ][j ][k ]=coeff‘Tg1 i 11.) 11k l-source-off+rxn_mw; 1 /* Zei = constant surface ‘/ if(((i>0)&&(i01&&(jal i ][j 11k 1.c; jac=grid_pt->Jlilljllk1; kav=1.0/6.0*('I‘h[ilLi+1]1kl+'I'h[illj- 111k1+Th1i+1 1L1 11kl+Thli- 11[j11k1+'1‘h[ilLi11k+11+Th1i11jllk-1 1); rho=density(TglilLi11k],Xg[ilLi]1k1); /"‘ Space derivative contributions "/ t1=a[01[0]*f_psi_psi(’1h,’1‘g,i,j,k,n1/pow(iac,2.01; t2=a[ ll[ 1]*f_eta_eta('lh,Tg,i,j,k,m1/pow(jac,2.0); t3=a12ll2]‘f_zei_zei(Th,Tg,i,j,k,11/pow(jac,2.0); t4=a1 01[ l ]‘(f_psi_eta(’l‘h,Tg,i,j,k,n,m)+f__eta_psi(Th,Tg ,iJ,k,n,m11/pow(jac,2.0); t5=a1 1 ll 2]*(f_eta_zei('1h,Tg,i,j,k,m,l1+f_zei_eta(’Ih,Tg, i,j,k,m,l1)/pow(jac,2.01; t6=a101121*(f_zei_psi('1h,Tg,i,j,k,l,n1+f_psi_zei('lh,Tg,i, j,k,l,n)1/pow(jac,2.0); l’ Shape contributions */ t7=kav"'(P(ij,k)*f_psi(Tg,i,j,k,n,m,l,21+Q(i,j,k)*f_eta(T g,i,j,k,n,m,l,2)+R(i,j,k)‘f_zei(Tg,i,j,k,n,m,l,2)1; /* Total is sume of different contributions including time derivative and contributions from nonlinear terms */ rhsl ilLi 11k ]=mulli 111 11k l‘Tgi i ][j ][k 1+delt*((rho- vol_f’rho_f)*kinetics_rate(Tg{i ]Li ]1k1,Xg[ i ][j ][k 11*(- delh1+MW_power(Tg[illj11k1,Xglil1j11k111+0.50*delt*(t1 +t2 +t3 +t4 +t5 +t6+t7 1; 1 1 1 1 free_3I)(Th,n,m,l); return; 1 f" Relaxation Routine For Temperature and Cure . This could be used as a relaxation routine for a single variable in a 3-D parabolic different ial equation. T (Temperature) is the main variable being solved for. In the case of only one variable pass a 'NULL' pointer in the place of X & Xg "/ find_next(t,grid_pt,conv,mul,T,X,Tg,Xg,Sf,n,m,l,nu,bc 1 double t; /"‘ Time ‘/ struct point *grid_pt; /"‘ Structural Parameters ast current grid */ double sssconv’sssmu],sssT,*ssx’sssTg,assxg’ssssf; /" conv : Convective heat transfer coefficents mul : Coefficient of time derivative T : Temperature; X : Cure; Tg,Xg : Temperature and cure evaluated at previous time step ; Sf : Sum of all source terms (sum of non dependent variable based terms and terms evaluated at previous time step ; */ int n,m,l;/"‘ Mesh sizes in three primary directions */ int nu,bc; /" nu : No of iterations to be performed bc : Indicates kind of BC ( 1 : Constant/Dirichlet 2 : Mixed/General/Robin ) ‘/ 1 extern double final_time,delt; extern int N,M,K; double "‘"I‘h,h,g,r,denm,t1,t2,norm(),boundary_surface_tem p(),corner_point_temp(),edge_point_temp(1,L_interior (),solve__kinetics( 1; int iJ,k,count,cnti,cntj,cntk,begi,begj,begk,1imit,ml,m2,m 3,bound(); double error,errord,Lop,dnl_dT,dLop,*"mem_alloc_3D(1; double "‘*’old_T,T_max,T_min,"a,jac,nonlinear_derivative() ,temp_check( 1; double thermal_conductivity(1,L_off_diag(),bou ndary_const_t emp(1,find_bou ndary_temp(1,ret,thermal_cyc1e(1; double fmd_min(),find_max(1,tol,rho,rxn_mw; /* Allocate memory for storing thermal conductivity/Coefficient of space derivative data ‘/ 'I‘h=mem_alloc_3D(n,rn,l); fill__conductivity('1h,T,X,n,m,l1; /* Calculate space steps in three principal directions for current grid*/ h=((double1(N- 1))/(n-11; g=((double)(M-1)1/(m-1); r=((double)(K- 1))/(1-1 1; errord=1.0; count=l; limit=nu; tol=l.0e-4; if((n==3)l | (m==3) I l(l==3)) 1 old_T=mem_alloc_3D(n,m,l1; errord=1; 1 else errord=1.0e-8; /* Gauss-Seidel Red-Black Relaxation Sweeps */ /* Outer most Loop Controlling no of sweeps ‘/ while((count<=limit1 I l(errord>tol1) 1 cnti=1; begi=0; if((n==3)l |(m==3)l l(l==3)) mat_copy(old_T,T,n,m,l); /* Red or Black Sweeps (l1"‘/ while(cnti<=21 1 l" Psi Direction */ for(i=begi;i01&&(j 01&&(i01| l(k01&&(j01&&(k01&&(i01&&(k<1-1))1 ret=find_boundary_temp(t,2,grid_pt,conv,Sf,Th,mul,T ,X,Ts,Xs,n.m,l.iJ,k,bc); ifl((i>01&&(iO1&&(j J1i11.111k1; a=grid_pt->alil1j11kl-C; denm:pow(jac,2.01; l’ Contribution from time derivative */ tl=mullilljl1k1; /" Contributions from space derivatives */ t2=0.50*delt*(al01101“(Th[i+11U11k1+'1‘h[i- 111,1 11k11/pow(h,2.01+a[ 11111*('1h[il[j+111k1+'1‘hli]1j- 1][k1)/pow(g,2.01+a1211 21*(1'111110 llk+11+Thlillj ][k- 111/p0w(r,2.01)/denm; /" Contribution from neighbouring points ‘/ Lop=L_off_diag(grid_pt,mul,Th,T,X,i,j,k,n,m,l1; ret=(StIilijllk1-Lop1/(t1+t2); 1 /" Check for consistency of calculated variable */ Tlillj][k]=temp_check(t,ret,Tg,Xg,mul,i,j,k1; 124 /" Calculate coupled cure variable using 4th order RK method */ Xlilljllk1=solve_kinetics('1‘1i][j11kl,X[ilUllkl.TslilUllk l, Xglilljlfkl); /* Update Thermal Conductivity 1"/ Thlillj11kl=thermal_conductivity(TIi11,1 l11(1,X1i11j11k1); 1 begk++; cntk++; 1 1 begj++; cntj++; 1 1 begi++; cnti++; 1 if((n==31l l(m==3)| l(l==3)) 1 mat_sub(old_T,T,old_T,n,m,l1; T_max=find_max(T,n,m,l); T_min=find_min(T,n,m,l1; errord=norm(old_T,n,m,l1/(T_max-T_min); 1 count++; 1 /* Free allocated memory for future use */ free_BD(Th,n,m,11; if((n==311 l(m==3)| l(l==3)) free_3D(old_T,n,m,l1; return; 1 P The Fully Multigrid V algorithm for a single variable ‘/ FMV(t,v,bc) double t; /* Time */ struct varray ‘v; /“ The V Structure containing the variables and shape 1 parameters for the current grid */ int bc; /" Nature of BC 1. Dirichlet, 2. Mixed */ 1 extern int nt1,nt2; /"‘ FMV iteration parameters */ int inkk; struct varray I“temp; /" V structure for the next level corase grid "7 double star’sssn’sssm; double "*rf,tru nc,err,norm(); double asstf’sssvf,sssvx,ssstx,sssrx; int vn,vm,vl,tn,tm,tl; double ***mem_alloc_3D(1; vn=v~>n; vm=v->m; v1=v->l; vf=v->T; =v->X; I" Perform at] relaxations if the current grid is coarsest i.e contains only one interior point "/ ifiv->next==NULL) 1 l“ printf("Time : %1f; FMV : Coarsest Grid Begins \n",t); ‘/ find_next(t,v->grid_pt,v->conv,v->multp,vf,vx,v- >Tg,v->Xg,v->Sf,vn,vm,v1,nt1,bc); 1 else 1 find_next(t,v->grid_pt.v->conv,v->multp,vf,vx,v- >Tg,v->Xg,v->Sf,vn,vm,vl,nt 1,bc); temp=v->next; tn=temp->n; tm=temp->m; tl=temp->1; tf=temp->T; tx=temp->X; I" Use restriction operators to restrict current grid variables to coarser grid */ Restrict(tf,vf,tn,tm,t11; Restrict(tx,vx,tn,tm,tl); f1=mem_alloc_3D(tn,tm,tl); f=mem_alloc_3D(tn,tm,tl); rf=mem_alloc_311(tn,tm,tl1; rx=mem_a110c_3[)(tn,tm,tl1; f2=mem_alloc_3 D(vn,vm,vl1; mat_00py(rf,tf,tn,tm,tl1; l’ Calculate the operated values L_2h(u_2h1 (Coarse grid operator) "/ Construct_L(t,temp->grid_pt,temp->conv,f,temp- >multp,tf,tx,tn,tm,tl,bc); /" Calculate the operated values L_h(u_h1 (Current grid operator) */ Construct_L(t,v->grid_pt,v->conv,f2,v- >multp,vf,vx,vn,vm,v1,bc1; /“ Restrict current grid operated solution fl = l_h- 2h(L_h(u_h11 "/ Restrictf f1 ,f2,tn,tm,tl1; /"' f = L_2h(u_2h1-l_h-2h(L_h(u_h)1 */ mat_sub(f,f,fl,tn,tm,tl1; v->trunc=norm(f,tn,trn,tl); /‘ f1 =L_2h(u_2h) ‘/ Construct_L(t,temp->grid_pt,temp->conv,f1 ,temp- >multp,tf,tx,tn,tm,tl,bc); Construct_L(t,v->grid_pt,v->conv,f2,v- >multp,vf,vx,vn,vm,v1,bc1; /" f2 = L_h(u_h1- Sf _h ‘/ mat_sub(f2,f2,v->Sf,vn,vm,vl); /" f = I_h-2h(L_h(u_h1- Sf_h) */ Restrict(f,f2,tn,tm,tl1; /" Sf_2h = L_2h(u_2h1- l_h-2h(L_h(u__h1 - Sf_h) */ mat_sub(temp->Sf,fl,f,tn,trn,tl); /" Recursive call to FMV using the next coarser grid "/ FMV(t,temp,bc); /' f1 = u_2h - 1__h-2h(u_h1 ‘/ mat_sub(fl,tf,rf,tn,tm,t1); /* f2 = I2h-h(u_2h - [h-2h(u_h) ‘/ Interpolate(f2,f1,tn,tm,tl); f" u_h = u_h + 12h-h(u_2h - Ih-2h(u_h) */ mat_add(vf,f2,vf,vn,vm,vl); /" Relax u_h ‘/ find_next(t,v->grid_pt,v->conv,v->multp,vf,vx,v- >Tg,v->Xg,v->Sf,vn,vm,vl,nt2,bc); P Free allocated memory ‘/ free_3D(rf,tn,tm,tl1; free_3 D(rx,tn,tm,t11; free_3[)(f,tn,tm,tl); free_3D(fl,tn,tm,tl1; 125 free_31)(f2,vn,vm,vl); 1 l’ printf("Time :%1f; FMV (N = %d) ends here\n",t,vn1; */ return; 1 /" Fully Approximate Storage Algorithm */ FAS(t,v,bc1 double t; /" Time ‘/ struct varray *v; /" V Structure for current grid */ int he; /“ Nature of BC I"/ 1 extern int nt0; /" Iteration parameter for FAS "/ extern double delt,T_initia1; double thermal_cycle(); int i,j,k,cnt; struct varray ‘temp; double g; double stsf,sssn,tasm; double "“rf,err,trunc,norm(); double l"""‘mem_alloc_3 110; double "'"sub_T,***tf,”"tx,"*vf,""vx; int vn,vm,vl,tn,tm,tl; vn=v~>n; vm=v->m; v1=v->l; vf=v->T; vx=v->X; if(v->next==NULL) 1 /" print_source(v->Sf,v->n,v->rn,v->11; print _grid(v->grid__pt,v->n,v->m,v->l1; */ for(i=l;i<=nt0;i++) FMV(t,v,bc); 1 else 1 temp=v->next; tn=temp->n; tm=temp->m; t1=temp->l; tf=temp->T; tx=temp->X; Restrict(tf,vf,tn,tm,tl); Restrict(tx,vx,tn,tm,tl); f1=mem_alloc_3D(tn,tm,tl1; f=mem_alloc_3D(tn,tm,tl); rf=mem_alloc_31)(tn,tm,tl1; f2=mem_alloc_3D(vn,vm,vl1; mat_copy(rf,tf,tn,tm,tl); P f1 =L_2h(u_2h1 ‘/ Construct_L(t,temp->grid_pt,te mp->conv,f1 ,temp- >multp,tf,tx,tn,tm,tl,bc1; l" 12 = L_h(u_h1 ‘/ Construct_L(t,v->grid_pt,v->conv,f‘2,v- >multp,vf,vx,vn,vm,v1,bc1; l" f2 = L_h(u_h1- Sf_h ’/ mat_sub(fl,f2,v->Sf,vn,vm,vl1; /* f2 = l_h-2h(L_h(u_h1- Sf_h) */ Restrict(f,f2,tn,tm,tl); /“ Sf_2h = L_2h(u_2h1- I_h-2h(L_h(u_h1 - Sf_h) ‘/ mat_sub(temp->Sf,fl,f,tn,tm,tl1; /* Recursive call to FAS over next coarser grid */ FAS(t,temp,bc1; mat_sub(fl,tf,rf,tn,tm,tl); Interpolate(f2,fl,tn,tm,tl); P u_h = u_h + 12h-h(u_2h - Ih-2h(u_h1 "/ mat_add(vf,f2,vf,vn,vm,vl); err=1; cnt=1; 8=0; /" Perform FMV ntO times on u_h(corrected) */ while(cnt<=nt()) 1 if(cntl=1) g=err/v->trunc; else s=l; FMV(t,v,bc); Construct__L(t,v->grid_pt,v->conv,f2,v- >multp,vf,vx,vn,vm,v1,bc1; mat_sub(f2,v->Sf,f2,vn,vm,vl); err=norm(f2,vn,vm,v11; cnt++; 1 /‘ Free allocated memory */ free_31er,tn,tm,tl); free_3D(f,tn,tm,t11; free_SD(fl,tn,tm,tl); free_3D(f‘Z,vn,vm,v1); 1 return; 1 /‘ To allocate memory and read in parameters into appropriate multigrid structures (linked lists) */ arrange_struct(t,grid_pt,T,X,n,m,l,bc) double t; /" Time ‘/ struct point ‘grid_pt;/"' Shape parameters for finest grid ‘/ double *"T,"*X; /* Temperature and cure or dependent variables on the finest grid */ int n,m,l,bc; /* Mesh sizes and nature of boundary conditions *I 1 extern double delt; extern double "*transfer_coeff; extern struct varray top; struct point *arrange_transform(); struct varray *ptr,*temp; double jac,**"‘mem_alloc_3 D(1,"a; int iik; if(fabs(t)<=l.0e-8) 1 /"' Ifthe time is start time then create all structures and link them */ /" Mesh sizes ‘/ top.n=n; top.m=m; top.l=1; top.next=NULL; P Copy finest grid values */ top.T=mem_alloc_3D(n,m,l1; 126 top.X=mem_alloc_3D(n,m,l1; top.Xg=mem__alloc_3D(n,m,l); top.Tg=mem_alloc_3D(n,m,l); top.conv=mem_alloc_3D(n,m,l); /“' t=0 => Dependent variables are at initial values */ Initial_Condition(top.T,top.X,n,m,l); /’ Copy heat transfer coefficient values ‘/ mat_copy(top.conv,transfer_coeff,n,m,11; top.Sf=mem_alloc_3D(n,m,l1; top.multp=mem_alloc_3D(n,m,l); top.grid_pt=grid_pt; ptr=⊤ /‘ For all elements for all structures upto the coarsest grid ’/ while((ptr->n>31&&(ptr->m>31&&(ptr->1>311 1 temp=(struct varray 1")malloc(sizeof(struct varray11; /" Mesh Size */ temp->n=(ptr->n+l)/2; temp->m=(ptr->m+11/2; temp->l=(ptr->l+l)/2; l‘ Allocate memory for dependent variable values at current time step and previous time step "/ temp->T=mem_alloc_3D(temp->n,temp->m,temp->l); temp->X=mem_alloc_3D(temp->n,temp->m,temp->l); temp->Xg=mem_alloc_3D(temp->n,temp->m,temp- >11; temp->Tg=mem_alloc_3D(temp->n,temp->m,temp- >1); /* Allocate memory for transfer coeeficient matrix */ temp->00nv=mem_alloc_3D(temp->n,temp->m,temp- >1); /* Restrict transfer coefficient matrix to the coarser grid */ RestrictItemp->conv,ptr->conv,temp->n,temp- >m,temp->l); /" Allocate memory for time derivative coefiicients and source terms I'/ temp->Sf=mem_alloc_3D(temp->n,temp->m,temp->l); temp->multp=mem_alloc_3D(temp->n,temp- >m,temp->l); /“‘ Calculate shape factors (x,y,z and alpha,beta and local jacobians 8 functions of psi,eta and zei for the coarser grid */ temp->grid_pt=arrange_transform(ptr- >grid_pt,temp->n,temp->m,temp->l); temp->next=NULL; /’ Set current structure to point to next coarser grid */ ptr->next=temp; ptr=ptr~>next; 1 1 /"' rho‘cp=> coefficient of time derivatives evaluated at previous time */ mat_copy(top.Xg,top.X,n,m,l); mat_copy(top.Tg,top.T,n,m,l); evaluate_multiplier(top.multp,top.Tg,top.Xg,top.n,top .m,top.l); /"‘ Sf(source terms) evaluated at previous time ‘/ set_rhs(t,top.grid_pt,top.conv,top.Sf,top.multp,top.’l‘,t op.X,top.n,top.m,top.l,bc); ptr=⊤ l’ One time only calculations overall all multigrid passes "/ while(ptr->next) 1 temp=ptr->next; /" Calculate source terms, time derivative coefficients and values of dependent variable at previous time step for all the grid sizes ‘/ Restrict(temp->multp,ptr->multp,temp->n,temp- >m,temp->l); Restrict(temp->Tg,ptr->Tg,temp->n,temp->m,temp- >1); Restrict(temp->Xg,ptr->Xg,tempo>n,temp->m,temp- >1); ptr=temp; 1 return; 1 /* Setting up shape parameters for coarse grid, given the same for the finer grid */ struct point ‘arrange_transform(f,n,m,l1 struct point "f; /" Shape structure for fine grid */ int n,m,1; /"‘ Coarse grid sizes */ 1 double calc_alpha Jac( 1; struct point ‘ret,’set _grid( 1; double **"‘tempx,"‘**tempy,"‘**tempz,***c,*"c1,*"‘*mem_all oc_3D(); int o,p,q,r,s,t,i,j,k; /‘* Calculate fine grid sizes "‘/ p=2"‘n-l; q=2*m-l; r=2“l-1; c=mem_alloc_31)(p,q,r); cl=mem_alloc_3D(n,m,l); ret=set_grid(n,m,l1; /“‘ Restrict x,y,z values to coarse grid ’/ Restrict(reto>x,f->x,n,m,l); Restrict(ret->y,f->y,n,m,l1; Restrict(ret->z,f->z,n,m,l); l’ Restrict jacobian values to coarse grid */ /* Restrict alpha values over the coarse grid */ P for(i=0;iJ[i][j][kl=ca1c_alphaJac(mt->a[i][j][k].c,ret- >x,ret->y,ret->z,n,m,l,iJ,k); 1 1 */ update _grid(ret,f->x,f->y,f->z,n,m,l1; Restrict(ret->J,f->J,n,m,11; for(s=0;s<3;s++) 1 for(t=0;t<3;t++) 1 for(i=0;ia111[j llkldSHtl; 1 1 1 Restrict(cl,c,n,m,11; for(i=0;ia[ile 11kl.C181ltl=Cl[il[j11k1; 1 1 1 1 /" print _grid(ret,n,m,l1; */ free_3 D(c,p,q,r1; free_3D(c 1,n,m,l 1; return(ret); 1 l’ Calculate numerical derivative (example) of all the non-linear terms (include specific phenomena as callable routines or functions */ double nonlinear_derivative(T,X) double T,X; 1 extern double vol_f,rho_f,delh; extern double k lo,k20,k30,E 1,E2,E3,X _ge1,B; double kl,k2,k3; double ret,dT,den_dT,pow_dT,kin_dT,den,kinetics_rate(),M W_power(1,density(1; den=density(T,X); dT= 1.0e-6*T; den_dT=(density(T+dT,X1-density(T-dT,X)1/(2.0*dT); ifIXTs)) ret=Ts; /* If heat addition is negligible and convection is out of the body and the new temperature is less than ambient then set new temperature equal to ambient temperature */ if((fabs(heat_addition1<=1.0e- 51&&(Tg[ilU11k1>=Ts)&&(f=Ts)&&(f>Tglile][k111 ret=Tg1UL111k]; /" lf heat addition is +ve and convection is into the body and the new temperature is less than old then set new temperature equal to old temperature */ ifl(heat_addition>01&&(Tg1illj11k1<=Ts)&&(f=Ts1&&(f>Tg1 i IL) 11 111)) ret=Tg1ilLi 11k]; return(ret); 1 E ‘15 |.|. This code essenytially deals with the generalized boundary conditions associated with the process model and calculates using fictitious boundary point assumptions the various coefficients and source terms involving the interior points. These are then used to obtain the next iterates of the boundary point values in the Gauss-Siedel relaxation procedure. #inclu de #include #include #include "grid.h" /“ Collect all coefficients of the fictitous dependent variable in the finite diference representation at the boundary using fictitious boundaru points assumptions */ double collect_coefilt,c,grid_pt,T,X,Th,conv,i,j,k,n,m,l) double t; /* Time */ int c; /* Parameter to indicate the type of boundary 1. Psi = constant 2. Eta = constant 3. Zei = constant */ struct point ‘grid_pt; /‘ Shape structure */ double ""T,"*X,"*Th,”"‘conv; /"‘ Temperature and Cure variables Space derivative coefficients and transfer coefficient matrices */ int i,j,k,n,m,l; /“ Grid coordinates (i,j,k) and mesh size ‘/ 1 double ret__coeff,"a‘jac,bef,dcf,h,g,r,kavg; extern int N ,M,K; int bound(),m1,m2,m3; a=grid_pt->alillj11kl.c; jac=grid_pt->Jlilljllkl; h=(double)(N-11/(n-11; g=( double1(M- 11/(m- 1 1; r=(double)(K-11/(l-1 1; switch(c) 1 case 1: 1 /" Psi = constant surface */ m1=bound(i,n1; /* Boundary normal evaluation = -1 Psi =0; = 1 Psi =Psi_max*/ /" Calculate average space derivative coefficient */ ifl((j>0)&&(j01&&(k01&&(i01&&(k0)&&(i01&&(ia1 1 ID 11k l.c; jac=grid_pt->J[i1U11k1; h=(double)(N-11/(n-1); g=( double )(M- 11/(m- 1 1; r=(double1(K-1)/(l-1); switch(c) 1 /" Psi = constant surface */ case 1: { m1=bound(i,n1; dcf=ml*a[0]10V(2.0*h); ret_off=T1 i-ml 111 11k ]- 1.0/dcf'(a[0][ 1 1*f_eta(T,i,j,k,n,m,l,21+a10 11 2 ]*f_zei(T,i,j ,k.n.m.l,2)); break; 1 /* Eta = constant surface ‘/ case 2: 1 m2=bound(i,m); dcf=m2*a[111 11/(2.0*g); ret_off='1‘1il[j-m211kl- 1.0/dcf"‘(a[ 1][ 01*f_psi(T,i,j,k,n,m,l,21+a[ 1][ 2]*f_zei(T,i,j, k,n.m,l.2)); break; 1 /“ Zei = constant surface */ case 3: { m3=bound(k,11; dcf=m3"‘a121[2l/(2.0"r); ret_off=T1 1le ]1k-m3 ]- l.0/dcf“(a[2110]‘f_psi(T,i,j,k,n,m,l,2)+a[ 211 1]"‘f_eta(T,i,j ,k,n,m,l,2)); break; 1 1 return(ret_of0; 1 /“‘ Calculate the source terms (due to ambient temperature) in the BCs usin the fictitious point assumptions ‘/ double collect_source(t,c,grid_pt,T,X,Th,conv,i,j,k,n,m,l1 double t; /" Time ‘7 int c; /* Surface parameter */ struct point *grid_pt; /"‘ Shape structure *I double ¥**T’***x,***%,¥*$conv; int itiikmml; 1 extern double delt; double ret_source,*‘ajac,bcf,dcf,h.g,r.Ts,thermal_cycle(1,kav 8; extern int N,M,K; int bound(),m1,m2,m3; a=grid_pt->a1ill.111k].c; jac=grid_pt->J[ 1111 11k]; h=(double)(N-11/(n-l1; g=(double)(M- 11/(m- l); r=(double1(K-l1/(l- 1 1; Ts=thermal_cycle(t+delt1; switch(c) 1 /* Psi = constant surface */ case 1: 1 /" Evaluate normal to boundary */ m1=bound(i,n1; /" Calculate space derivative coeffn. average ”"/ ifl((j>01&&(j01&&(k<1-11)) kavg=0.25*(ThlilLil1k+1]+Th1i11jllk-1]+Th1i1Li- 111kl+TTrlilLi+lllk11; else kavg='l‘hlilljllk]; /* Calculate source denominator contributions */ bcf=m1“'kavg/(iac"‘sqrt(a[0][ 011); dcf=m1*a[01101/(2.0*h); break; 1 /" Eta = constant surface ‘/ case 2: [ m2=bound(j,m); if(((i>01&&(i01&&(k0)&&(i01&&(ja[ i ][j ][k ].c; jac=grid_pt->JIi]le[k]; h=(double)(N-11/(n-l); g=(double)(M-11/(m-1); r=(double)(K-11/(l-11; switch(c) 1 /" Psi = constant surface ‘/ case 1: 1 m1=bound(i,n1; /"‘ evaluate boundary normal */ /" Calculate surface average of space derivative coeffn. ‘/ kavg=0.25"‘(ThIillj+l 11k ]+Thl i ][j- 1 11kl+'l‘h[i11j 11k- 11+ThlilUllk+1 1); l‘ Collect appropriate coefficients for T(i,j,k1 from boundary conditions*/ coeff_i=collect_coefl(t,c,grid_pt,T,X,T‘h,conv,i,j,k,n,m,l 1; /"‘ Calculate averages for space derivatives ‘/ k_eta_avg=0. 50*(Th1 i 1Lj+l 11k1+Thlillj-1][k1); k_zei_avg=0.50*(Th[ilU11k+1]+'1‘h[illj11k-l l); l“ Collect coefficients from double derivatives and forcing terms ‘/ psi_psi_coeff=a[01101*ThlilUllkV(pow(jac‘h,2.0))*(- 2+coeff_i); eta_eta_coeff=-2.0"a[ 1 ]1 1 ]‘k_eta_avg/pow(iac*g,2 .01; zei_zei_coeff=-2.0‘a[2][2l*k_zei_avg/pow(iac"‘r,2.0); force_psi_coeff=m1*kavg"' P(i,j,k)/(2.0*h)*coeff_i; coeff=0.50"delt"(psi_psi_coeff+eta_eta_coeff+zei_zei_ coeff+force_psi_coeff); break; 1 case 2: l m2=bound(j,m); kavg=0.25*('1'hli+l 111 11k 1+Thli- 1 ]Li 11k 1+Thli ][j ][k- 11+T'111ilLi11k+11); /* Collect appropriate coefficients for T‘(iJ,k) from boundary conditions*/ coeff J=collect_coefi1t,c,grid_pt,T,X,'l‘h,conv,i,j,k,n,m,l 1; P Calculate averages for thermal conductivities ‘/ k_psi_avg=0.50‘(Th[i+1 ][j][k]+Th[i-11[jl[kl1; k_zei_avg=0.50"(ThIi][j][k+l]+Th1i][j][k-11); /* Collect coeficients from double derivatives and forcing terms */ psi_psi_coeff=-2.0"‘a[01101*k_psi_avg/pow(iac*h,2.0); eta_eta_coeff=a1 1 ][ 1]"‘Thli ]U ][k l/(pow(jac"g,2.0)1"‘(- 2+coeff_j); zei_zei_coeff=-2.0*a1 2 ][ 2 l‘k_zei_avg/pow(iac*r,2 .01; force_eta_coeff=kavg*m2*Q(iJ,k1/(2.0‘g1‘coefij; coeff=0.50*delt‘(psi_psi_coeff+eta_eta_coeff+zei_zei_ coeff+force_eta_coeff); break; 1 case 3: { m3=bound(k,l); kavg=0.25*('1'hli+1]Li11k1+Thli-1 1U11k1+Th1i11j- 111kl+ThIilLi+lllk111 /* Collect appropriate coefficients for 'l‘(i,j,k1 from boundary conditions*/ coefth=collect_coeff(t,c,grid_pt,T,X,Th,conv,i,j,k,n,m, l); /"‘ Calculate averages for thermal conductivities */ k_psi_avg=0.50*(Th1i+l 11j11k1+Th[i-1 lLil1k1); k_eta_avg=0.50*(Thlil[j+111kl+ThIilU-l111(1); l“ Collect coeficients from double derivatives and forcing terms ‘/ psi_psi_coeff=-2.0“a[ 0][ 0 1*k__psi_avg/ pow(j ac*h,2.01; eta_eta_coeff=-2 .0*a[ 1 ][ 1]‘k_eta_avg/pow(jac"‘g,2.01; zei__zei_coeff=al 2][2]*Th[i][j]1k 1/(pow(jac*r,2.0))*(- 2.0+coeff_k); force_zei_coeff=kavg*m3"‘R(i,j,k1/(2.0’r1*coeff_k; coeff=0.50*delt"‘(psi_psi_coeff+eta_eta_coeff+zei_zei_ coeff+force_zei_coefi1; break; 1 1 ret_var=-coe ff; /* printf("Coeff (%d,%d,%d1 : %6.31f \n",i,j,k,coeff); */ return(ret_var1; 1 double compute_boundary_off_diagonal(t,c,grid_pt,T,X,Th,co nv,mul,i,j,k,n,m,l) double t; int c; struct point *grid_pt; double ***T,***X,"""*Th,***conv,‘**mul; int i,j,k,n,m,l; 1 extern double delt; extern int N ,M,K; int bound(),m1,m2,m3; double f_psi(),f_eta(1,f_psi_psi(1,f_psi_eta(),f_psi_zei(),f_eta_e ta(),f_eta_zei(),f_zei_zei(),P().Q().R11; double h,g,r,"‘*a,jac; double kavg,k_psi_avg,k_eta_avg,k_zei_avg; double coeff_i_plus,coeff_i_minus,coeff_j_plus,coefLLminus, coeff_k_,plus,coeff_k_minus; double off_i,off _j,off_k,off_i_p1us,off_i_minus,of f J_plus,ofi J_ minus,off_k_plus,off_k_minus; double psi_psi_off,psi_eta_off,psi_zei__off,eta_eta_off,eta_zei_ off,zei_zei_off,force_psi__off,force_eta_off,force__zei_off; double off,ret_var; a=grid_pt->a[i]lj ][k ].c; jac=grid_pt->J[illj11k1; hddoubleXN- 11/(n- 1 ); g=( double )(M- 1 )/( m- 1 1; r=(double)(K- 11/(1- 1 ); switch(c) 1 case 1: l m1=bound(i,n1; /“ Collect appropriate coefficients from boundary conditions*/ coeff_j_plus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i,j+l .k.n,m.l); wefij_minus=mllect_coefflt,c,grid_pt,T,X,Th,conv,i,j - 1,k,n,m,l); coeff_k_plus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i,j,k +1,n,m,l); coeff_k_minus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i, i.k-1.n,m.1); l" Collect relevant off diagonal terms from boundary condition */ off_i=collect_ofi‘(t,c,grid_pt,T,X,Th,conv,i,j,k,n,m,l1; off 4' _plus=collect_off(t,c,grid_pt,T,X,Th,conv,i,j+1,k,n .mJ); off _j_minus=collect_off(t,c,grid _pt,T,X,Th,conv,i,j— 1.k.n,m.l); off_k_plus=collect_ofi(t,c,grid_pt,T,X,Th,conv,i,j ,k+ 1, n,m,l); off_k_minus=collect_off(t,c,grid_pt,T,X,Th,conv,iJ,k- 1,n,m,l1; /" printf("Coeff_j_plus : %6.31f \n",coeffj_plus); printf("CoeffJ_minus :%6.3lf \n",coeffJ_minus); printf("Coeff_k_plus : %6.31f \n",coeff_k_plus); printf("Coe1f_k_minus : %6.31f \n",coeff_k_minus); printf("l : %d; J : %d ; K : %d \n",i,j,k); printf("off_i: %6.31f \n”,off_i); printf("ofLLplus : %6.31f \n",off_j_plus1; printf("ofLLminus : %6.31f \n",offJ_minus1; printf("off_k_plus : %6.31f \n",off_k_plus1; printf("off_k_minus : %6.31f \n",off_k_minus); printf("Tml :%6.31f; Tj-l :%6.31f; Tk-l : %6.31f ; Tj+1 : %6.31f; Tk+1 : %6.3lf\n",T[i-m11U11k],T[i][j- 111k1.'1‘lilLillk-11,T111Li+111k1,'1‘1i]U11k+11); */ /* Calculate averages for thermal conductivities */ kavg=0.25"‘(Thli]Lj+ll[k]+ThlilLi-l][k]+Thli11j11k- 11+ThlilLi11k+1111 k_eta_avg=0.50"(Th[i 1Li+ 1][kl+Th[i]lj-1 ][k l); k_zei_avg=0.50*(ThIillj11k+11+Thlillj11k-l 1); /* Collect off diagonal terms from double derivatives and forcing terms */ psi_psi_off=a[ 0][0]"Th[i][j][k]/pow(iac*h,2.0)*(T1i- ml]U11k]+off_i); eta_eta_off=a( l][ 1 ]*k_eta_avg/pow0 ac‘g,2 .01*(T1 i le+ lllkl+Tlil[j-ll1kl); zei_zei_off=al2][2]*k_zei_avg/pOW(iac*r,2.01"(Tli][j 11k +1]+T[ilLi11k-11); psi_eta_off =2 *al 011 1 ]‘k_eta_avg/(4 .0“ pow(jac,2.01*h* g)*(ml*(coeff J_plus*’1‘[il[j+1l1kl+off J_plus1-m1*Tli- m111j+111k l+ml""l‘11-m1 ][j-111k ]- m1‘(coefij_nfinus*TTilU-l][kl+off_j_minus11; eta_zei_off=2.0"al 1][ 2 Vpoanc,2.0)* f_eta_zei(Th,T.iJ, k,m,l1; psi_zei__off=2 "al 011 2 1*k_zei__avg/(4.0*pow(jac,2.01*h*r )‘(m l*(coeff_k_plus“Tli ][j 11k+1 ]+off_k_plus)-m1*T[i- mlllj11k+l 1+ml‘T‘1i-m111jl1k-ll- m1"‘(coeff_k_minus*T[i ][j ][k- 1 ]+off_k_minus11; force_psi_off=kavg*P(i,j,k1/(2.0*h)*ml"(off_i-Tli- m1 1U11k1); force_eta_off=kavg*Q(iJ,k)‘f_eta(T,iJ,k,n,rn,l,2); force_zei_off=kavg*R(i,j,k1*f_zei(T,iJ,k,n,m,l,21; /"‘ Compute boundary temperature using source terms, off diagonals and coefficients ‘/ off=psi_psi_off+eta_eta_off+zei_zei_off+psi_eta_off+e ta_zei_off+psi_zei_off+force_psi_off+force_eta_off+for ce_zei_off; break; 1 case 2: l m2=bound(j,m1; /"‘ Collect appropriate coefficients from boundary conditions*/ coeff_i_plus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i+1,j ,k,n,m,l1; coeff_i_minus=collect_coefi‘(t,c,grid_pt,T,X,Th.conv,i- 1,j,k,n,m,l); coeff_k_plu s=collect_coeff(t,c,gri d_pt,T,X,’Ih,conv,i,J,k +1,n,m,l); coeff_k_minus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i, j,k-l,n,m,l1; /* Collect relevant off diagonal terms from boundary condition "/ off J=collect_ofi(t,c,grid_pt,T,X,Th,conv,i,j,k,n,m,l1; off_i_plus=collect_ofi(t,c,grid_pt,T,X,Th,conv,i+1,j,k,n .m.1); off_i_minu s=collect_ofl(t,c,grid_pt,T,X,Th,conv,i- l,j,k,n,m,l); off_k_plu s=collect_off(t,c,grid_pt,T,X,Th,conv,i,j,k +1, n,m,l); off_k_minu s=collect_off(t,c,grid_pt,T,X,Th,conv,i,j,k- l,n,m,11; /"‘ Calculate averages for thermal conductivities "‘/ kavg=0.25*(Th[i+lllj11k1+Th[i-11U11k1+Thlilljllk- 11+Thlillj 11k+1 l); k_psi_avg=0.50*('lh[i+11U11kl+Thli-l1L111k11; k_zei_avg=0.5()*(Thli11j11k+11+Thli 111 11k- 1 11; /" Collect off diagonal terms from double derivatives and forcing terms ‘/ psi_psi_off=a[0][0]*k_psi_avg/pow(jac‘h,2.0)*(T1i+1]Lj 11k1+T1i-111j1[k11; eta__eta_off=a[ 111 1]*Th[illj]1k VpowQac‘g,2.01*(T1i]Lj- m211kl+off_i1; zei_zei__off=a[ 2 ]12]‘k_zei_avg/pow(jac"r,2.0)*(T[ illj 11k +ll+'1‘1111j11k-1 11; psi_eta_off=2*a[ 1][01*k_psi_avg/(4.0‘pow(jac,2.01*h"g 1"‘(m2*(coeff_i_plus*T1i-i-1 llj]1k1+off_i_plus1- m2‘T1i+llLi-m211kl+m2*'l‘li-l ][j-m211k1- m2*(coeff_i_minus*'l‘1i—1 1L1 ][k]+off_i_minus)1; eta_zei_off=2"‘a[ l1121‘k_zei_avg/(4.0*pow(jac,2.01"‘g*r 1*(m2*(coeff_k_plus*T[i][j1[k+1 ]+off_k_plus1. m2*Tiil[j-m211k+1l+m2*T1illj-m211k-l1- m2*(coeff_k_minus"T[i][j11k-ll+off_k_minus11; psi_zei_off=2.0“a10I12VpowQ‘ac,2.01“f_psi_zei('lh,T,iJ, k,1,n); force_psi_off=kavg"‘P(i,j ,k1"f_psi(T,i;j,k,n,m,l,21; force_eta_off=kavg*Q(i,j,kll(2.0‘g1"‘m2*(offJ-T[ile- m211k1); force_zei_off=kavg*R(i,j,k)‘f_zei(T,i,j,k,n,m,l,2); l" Compute boundary temperature using source terms, off diagonals and coefficients "/ off=psi_psi_off+eta_eta_off+zei_zei_off+psi_eta_off+e ta_zei_off+psi_zei_off+force_psi_off+force_eta_off+for ce_zei_off; break; 1 case 3: l m3=bou nd(k,11; /"‘ Collect appropriate coefficients from boundary conditions*/ coeff_i_plus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i+1,j ,k,n,m,l); coeff_i dminus=collect_coefi1t,c,grid_pt,T,X,T‘h,conv,i- 1J,k,n,m,l); coeff_j_plus=collect_coeff(t,c,grid_pt,T,X,Th,conv,i,j+1 ,k,n,m,l1; weffj_minus=collect_coefi(t,c,grid_pt,T,X,Th,conv,i,j ' lrkrnrmvl)r /“ Collect relevant off diagonal terms from boundary condition */ off_k=collect_ofilt,c,grid_pt,T,X,Th,conv,i;j,k,n,m,l1; off_i_plus=collect_off(t,c,grid_pt,T,X,Th,conv,i+1,j,k,n ,m,l); off_i_minus=collect_ofi(t,c,grid_pt,T,X,Th,conv,i- 1,j,k,n,m,11; off _j_plus=collect_off(t,c,grid_pt,T,X,Th,conv,i,j+1,k,n .m.l); off _j__minus=collect_off(t,c,grid_pt,T,X,Th,conv,i,j- 1,k,n,m,l); l" Calculate averages for thermal conductivities ‘/ kavg=0.25"‘(Th1i+llLi]1kl+Thli-11U11kl+'1'hlillj- 111k 1+Thl i 1Li+l 111(1); k_psi_avg=0.50*(Th[i+l 1U11k1+Th1i- llU l1kl); k_eta_avg=0.50*(Th[ile+l ][k1+Th[i]U-1][k]1; /* Collect off diagonal terms from double derivatives and forcing terms ‘/ psi_psi_off=a[0][0]*k_psi_avg/pow(jac‘h,2.0)‘(T1i+1][j 11kl+'l‘li-lllj11kl); eta_eta_off=a[ 1 ][ 1]*k_eta_avg/pow(jac*g,2.01“(T[i lLi+ 1 11k1+Tlillj-1 111(1); zei_zei_off=a[ 2 ][ 2 1*ThIile 11k l/powfiac‘r,2.0)*(T1i1Lj 11k- m31+off_k); psi_eta_off=2.0*al 1][0l/powUac,2.01*f_psi_eta('lh,T,iJ, k,n,m); eta_zei_off=2"‘al l112l‘k_eta_avg/(4.0’powfiac,2.01*g*r )‘(m3‘(coeff _i_plus“‘T[ i][j+1 ][k]+off J_plus1- m3"’l‘1i1[j+1 11k-m3]+m3*’1‘[i]lj- 1 11k-m3]- m3*(coefij_minus‘T1ilU-1][kl+off_j_minus1); psi_zei_off=2*al 2 ][ 0 ]*k_psi_avg/(4. 0" powfi ac,2 .01*g*r 1*(m3"‘(coeff_i_plus"T[i+l 1U][k]+off_i_plus)- m3”’11i+11[j][k-m3]+m3*’l‘[i-llLillk-m31- m3‘(coeff_i_minus"‘Tli-1 llj]1k1+off_i_minus11; force_psi_off=kavg"‘ P(iJ,k1’f_psi(T,i,j,k,n,m,l,21; force_eta_off=kavg*Q(i,j,k1*f_eta(T,iJ,k,n,m,l,21; force_zei_off=kavg‘R(i,j.kV(2.0“‘r1"m3"‘(off_k-T[ i 1L1 ][k- m3l); /" Compute boundary temperature using source terms, off diagonals and coefficients */ off=psi_psi_off+eta_eta_off+zei_zei_off+psi_eta_off+e ta_zei_off+psi__zei_off+force_psi_off+force_eta_off +1hr ce_zei_off; break; 1 1 ret_var=-0.50"delt*off; P printf("Off (%d,%d,%d1 : %6.31f \n",i,j,k,ofi1; */ return(ret__var); 1 double combined_boundary_sou rce(t,c,grid_pt,T,X,TIi,conv,m ul,i,j,k,n,m,l) double t; int c; struct point *grid_pt; double *"T,***X,***Th,"’conv,*"‘*mul; int i,j,k,n,m,1; 1 extern double delt; extern int N,M,K; int bound(),m1,m2,m3; double P(),Q(),R(),collect_source( 1; double h,g,r,**a,jac; double kavg,k_psi_avg,k_eta_avg,k_zei_avg; double source_i,source_j,source_k,source__i_plus,source_i_mi nus,sourceJ_plus,source_jmminus,source_k_.plus,sour ce_k_minus; double psi_psi_source,psi_eta_source,psi_zei_source,eta_eta_ source,eta_zei_source,zei_zei_source,force_psi_source ,force_eta_source,force_zei_source; double source,ret_var; a=grid_pt-> a[ i 1L1 11k 1.c; jac=grid_pt->Jl i 111 11k 1; h=(double)(N- 11/(n- l 1; g=(doubleXM-1V(m1); r=(double)(K- 11/(1- 1 1; switch(c) 1 case 1: 1 m1=bound(i,n1; /* Collect relevant source diagonal terms from boundary condition ‘/ source_i=collect_source(t,c,grid_pt,T,X,'lh,conv,i,j,k,n ,m.1); source _j_plus=collect_source(t,c,grid_pt,T,X,Th,conv,i J+l.k,n.m,l); source J_minus=collect_source(t,c,grid_pt,T,X,Th,con Vrirj' l,k,n,m,l); source__k__plus=collect_source(t,c,grid_pt,T,X,Th,conv, i,j,k+1,n,m,l1; source_k_minus=collect_source(t,c,grid_pt,T,X,Th,co nv,i,j,k-1,n,m,l); /* Calculate averages for thermal conductivities */ kavg=0.25*('l‘hlilLi+111k1+Thlile-1llkl+Thlilellk- 1 1+Th1111j 11k+1 11; k_eta_avg=0.50*(Th[il[j+l][k]+Th[i]Li-111k]1; k_zei_avg=0.50*(’1‘h[ilL111k+1]+Th[i1U11k-1l); /“ Collect source diagonal terms from double derivatives and forcing terms "‘/ '43 'JJ psi_psi_source=a[ 011 0]*Th[i ][j ][k l/poanc*h,2.0)*sour ce 1; psi_eta_source=2*a[ 011 11*k_eta_avg/(4.0*pow(jac,2.01 ‘h*g1*ml‘(sourcej_plus-sourcej_minus); psi_zei_source=2*a10][21*k_zei_avg/(4.0*pow(jac,2.0)‘ h’r)*m1*(source_k_plusosource_k_minus); force_psi_source=kavg“‘ P(i,j,k1/(2.0*h1*m1*source_i; /* Compute boundary temperature using source terms, source diagonals and coefficients */ source:psi_psi_source+psi_eta_source+psi_zei_sourc e+force_psi_source; break; 1 case 2: { m2=bound(i,m1; /* Collect relevant source diagonal terms from boundary condition */ source _j=collect_source(t,c,grid_pt,T,X,Th,conv,i,J,k,n ,mJ); source_i_plus=collect__source(t,c,grid_pt,T,X,Th,conv,i +1J,k,n,m,l); source_i_minus=collect_source(t,c,grid__pt,T,X,Th,con vri' 1,j,k,n,m,1); source_k_plus=collect_source(t,c,grid_pt,T,X,Th,conv, i,j,k+1,n,m,11; source_k_minus=collect_source(t,c,grid_pt,T,X,Th,co nv,iJ,k- 1,n,m,l1; /" Calculate averages for thermal conductivities */ kavg=0.25*(ThIi+11[jl1k1+Th[i-11U11k1+ThlilL111k- 11+ThlilLi11k+ll); k_psi_avg=0.50*(Th[i+l 11.1 ][k]+Th[i- l 111 11k 1); k_zei_avg=0.50*(Th[i]Lj11k+ll+ThlilLi11k-11); /* Collect source diagonal terms from double derivatives and forcing terms */ eta_eta_source=a[ 111 l]*Th[ illj ][k Vpowfiac‘g,2.0)*sour 09.1; psi_eta_source=2*a[ 1][0]*k_psi_avg/(4.0*pow(iac,2.0) *h‘g)‘m2"'(source_i_plus-source_i_minus); eta_zei_source=2 ‘al 1 112 ]‘k__zei_avg/(4 .0*pow(jac,2.0) *g‘r)*m2*(source_k_plus-source_k_minus); force_eta_source=kavg*Q(i,j,k1/(2.0*g1*m2"‘source _j; /* Compute boundary temperature using source terms, source diagonals and coefficients ‘/ source=eta_eta_source+psi_eta_source+eta_zei_sourc e+force_eta_source; break; 1 case 3: 1 m3=bound(k,l); /* Collect relevant source diagonal terms from boundary condition */ source_k=collect_source(t,c,grid_pt,T,X,Th,conv,i,j,k, n,m,l); source_i__plus=collect_source(t,c,grid_pt,T,X,Th,conv,i +1,j,k,n,m,l); source_i_minus=collect_source(t,c,grid_pt,’l‘,X,Th,con v,i-1,j,k,n,m,l1; source _j_plus=collect_source(t,c,grid_pt,T,X,Th,conv,i ,j+1,k,n,m,l); source J_minu s=collect_source(t,c,grid_pt,T,X,Th,con v,i,j-1,k,n,m,l1; /* Calculate averages for thermal conductivities */ kavg=0.25*(Th[i+11U11k1+Thli-111J11k1+'1h[illj- 111kl+Thli1U+111kl); k_psi_avg=0.50*(’l‘h[i+1 1U11k1+Th1i-1 llj 111(1); k_eta_avg=0.50*(Thli]1J+1 11k1+Th1i1Li-111k]); /* Collect source diagonal terms from double derivatives and forcing terms ‘/ zei_zei_source=a[21121*Th1i11j11k l/pow(jac*r,2 .0)*sour ce k; eta_zei_source=2*a[ 211 11*k_eta_avg/(4.0‘pow0ac,2.01 ‘g‘rl‘mB‘(sourcequlus-sourcequinus); psi_zei_sou rce=2*a[ 21101‘k_psi_avg/(4.0*pow(jac,2.0)"‘ h‘r1‘m3"(source_i__plus-source_i_minus1; force_zei_source=kavg"‘R(i,j,k1/(2.0‘r1‘m3*source_k; /* Compute boundary temperature using source terms, source diagonals and coefficients “'l source=zei_zei_source+eta_zei_source+psi_zei_source +force_zei_source; break; 1 1 ret_var=-0.50*de1t*source; /* printf("Source (%d,%d,%d1 : %6.31f \n",i,j,k,souree); */ return(ret_var); 1 double boundary_operator(t,c,grid_pt,T,X,’lh,conv,mul,i,j,k,n .mJ) double t; int c; struct point *grid_pt; double ***T,***X,*"Th,***conv,***mul; int i,j,k,n,m,l; 1 double coeff,off,source,compute_boundary_coefficient(1,comp ute_boundary_off_diagonal(1,comhined_bou ndary_sou rce(),ret_operator; double rho,rxn_mw,kinetics__rate(),MW_power(),density(1; extern double rho_f,vol_f,delh; coeff=compute_boundary_coefficient(t,c,grid_pt,T,X,T h,conv,mul,i,j,k,n,m,l); off=compute_boundary_off_diagonal(t,c,grid_pt,T,X,T h,conv,mul,i,i,k,n,m,l1; ret_operator=( mu lli ][j ][k 1+coeff)*'l‘[ i ][j ][k ]+off; l“ printf("Operator(%d,%d,%d) : %6.31f \n",ij,k,ret_operator1; */ return(ret_operator1; 1 double calculate_boundary_variable(t,c,grid_pt,T,X,Th,conv, mul,Sf,i,j,k,n,m,l) double t; int c; struct point *grid__pt; double **"T,"‘**X,"**Th,***conv,***mul,***Sf; int iJ.k,n,m.l; 1 double coeff,off,source,compute_boundary_coefficient(1,comp ute_boundary_off_diagonal( ),combined_boundary_sou rce(1,var; double rho,rxn_mw,kinetics_rate(1,MW_power(),density(); extern double rho__f,vol_f,delh; coeff=mulli][j][k]+compute_boundary_coefficient(t,c,g rid_pt,T,X,Th,conv,mul,i,j,k,n,m,l1; off=compute_boundary_off_diagonal(t,c,grid_pt,T,X,T h,conv,mul,i,j,k,n,m,l); var=(Sfli]Lj][k]-ofl1/coeff; /‘ printf("T(%d,%d,%d) : %6.31f; X(%d,%d,%d1 : %6.31f \n",i,j,k,'1‘lilL1llkl,i,j,k,X[il[jl[kl); printf("Coefi(%d,%d,%d) : %6.31f\n",i,j,k,coeff1; printf("Off(%d,%d,%d1 : %6.31f \n",i,j,k,off1; printf("Sfl%d,%d,%d) : %6.31f \n",i,j,k,Sflilljl1k11; printfi"T(%d,%d,%d) : %6.31f \n",i,j,k,var); */ return(var); 1 This code contains all the routines associated with the calculation of the various properties of the material which are assumed to be functions of temperature and cure. The user who wishes to feed in his or her own functionalities into the routines can do so. Remember that the appropriate input prompts that the program currently requires from the input file (PROPERTYDAT) should also be modified. Further note that the the lines of the code prompting the program to read the various lines associted with the property inputs present in process_in_out.c should also be modified. #inclu de #include l" Calculate density using tempoerature and cure values at a given point "I double density(Temp,Cure) double Temp,Cure; 1 extern double vol_f,rho_f.rho_m; double rf,rm,ret; rf=rho_f; rm=rho_m; ret=( 1.0-vol_f1“rm+vol_l‘rf; f" printf("Density : %lf \n",ret1; ‘/ return(ret); 1 /* Calculate specific heat values using temperature and cure values at a given point ‘/ double specific_heat(Temp,Cure1 double Temp,Cure; 1 extern double rho _f,rho_m,vol_f ; double den,wf,cpm,cpf,ret,fiber_sp_heat(1,matrix_sp_heat( 1; cpm=matrix_sp_heat(Temp,Cure1; cpf=fiber_sp_heat(Temp,Cure1; den=density(Temp,Cure1; wf=rho_f"‘vol_f/den; ret=cpm‘(1.0-w0+wf*cpf; /"‘ printf("Sp. Heat : %lf \n",ret1; ‘/ return(ret); 1 /* These are specific functions depending upon the property to charactersic correspondence for the material and can be altered by the user depending upon his property predictor ‘/ double fiber_sp_heat(T,X1 double T,X; 1 extern double Cp_f; double ret; ret=Cp_f; return(ret); 1 'JI double matrix_sp_heat(T,X) double T,X; 1 extern double Cp_m; double ret; ret=Cp_m; return(ret); 1 double matrix_conductivity(T,X) double T,X; 1 extern double km; double ret; ret=km; return(ret); 1 double fiber_conductivity(T,X) double T,X; 1 extern double kf; double ret; ret=kf; return(ret); 1 double thermal_conductivity(Temp,Cure) double Temp,Cure; 1 extern double vol_f; double thm,thf,matrix_conductivity(1,fiber_conductivity(1,ret; thm=matrix_conductivity(Temp,Cure1; thf=fiber_conductivity(Temp,Cure1; ret=thm"( 1 .0.vol_f)+thf‘vol__f; /* printf("Th. Conductivity : %lf \n",ret1; */ return(ret); 1 fill_conductivity(Th,T,X,n,m,l) double sssmssspsssx; int n,m,1; 1 int i,j,k; extern double kf,km; double thermal__conductivity(1; for(i=0;i #include #include #include ”grid.h" I“ Read in RXN rate equation data */ read_rxn_data() 1 extern FILE ’“prop_f; extern double k10,](20,k30,E1,E2,E3,X_ge1,B; extern double delh; fscanflprop_f,"%lf',&k lo); fscanflprop_f,"%1f',&El); fscanflprop_f,"%lf',&k201; fscanflprop_f,"%lf',&E2); fscanflprop_f,"%lf’,&k30); fscanflprop_f,"%lf',&E3); fscanflprop_f,"%lf',&x_gel); fscanflprop_f,"%lf',&B); fscanflprop_f,"%lf',&delh); return; 1 /* Read in property data */ read_property_data() 1 extern FILE *prop_f; extern double vol_f; extern double Cp_f,Cp_m,rho_f,rho_m,kf,km; fscanflprop_f,"%lf',&kfi; fscanflprop_f,"%lf‘,&km); fscanflprop_f,"%lf',&Cp_.0; fscanf(prop_f,"%lf',&Cp_m); fscanflprop_f,"%lf',&vol_0; fscanflprop_f,"%1f',&rho_0; fscanflprop _f,"%lf',&rho_m1; return; 1 /* Read in transfer coefficient data if required */ read_transfer_coefl(n,m,l) int n,m,1; 1 extern FILE *proc_inp; extern double ***transfer_coeff; char be; double ret,""'menLalloc_3D(1; double val; int cnt,stp_i,stp_j,stp_k,ij,k; transfer_coeff=mem_alloc_3D(n,m,l); fscanflproc_inp,"%c",&bc); for(cnt=0;cnt<3;cnt++) 1 switch(cntl 1 case 0: [ stp_i=n-1; stPJ=1; stp_k=1; for(i=0;i #include double thermal_cyclc(t) double t; 1 extern double ”thermal_cycle_data; extern int no; double ret; int i,found; i=0; found=0; while((i=thermal_cycl e_data[i 1101)) 1 ret=(t- thermal_cycle_data[ i ][0 ])‘thermal_cycle_data[ 1+ 111 1 1 +thermal_cycle-data[i][2]; fou M: 1 ; 1 i++; 1 iflfou nd==0) ret=(t- thermaLcycle_data[ no 11 01)‘thermal_cycle_dat a[ no][ 1 1 +thermal_cycle_data[ no ][ 2 l; “Trek-=0) ret=thermal.cycle_datal no ][ 2 ]; return(ret); 1 InitiaLCondition(T.X,n.m.l) double “‘T,”“X; int n,m,1; 1 extern double T_initiaI,X.initial; int iJJr; for(i=0;i Sinclude I38 Cdefine gas__const 8.314 /" Function which calculates the rate of rxn using a given rate equation (should be replaced within this routine when using a different rate equation ) *I double kinetics_rate('I‘.X) double T,X; /‘ Temperature and cure ‘I 1 extern double klo.k20.k3o.El.EZ.E3.X_gel,B; double ret.k 1,k2,k3,k4; iKX>1.000) ret=0.00; else 1 klsk lo‘exp(-E 1/(gas_const"‘T)); k2-k20‘exp(.82/(gas_const"l‘)); k3=k3o‘exp(-Eal(gas_const"l‘)1; iKX1.0000) reb 1.0000; return(ret); l W The following set of routines are involved in the input and manipulation of surface coordinate information depending upon the availability of coordinate data (either surface or edge information). Iinelude Oinclude Iinclude Iinclude /’ Read in surface by surface information of the coordinates of the body '/ read_input(ch.fx,fy.fz) int ch; double ooo&.oo0fy.osofz; 1 extern int N,M,K; extern FILE ‘inp_f; extern double scale_x.scale_y.scale_z; int iJ,k,cnt; double "‘mem_alloc_3D(): iflch==2) 1 I‘ Read in surface values for x,y,z‘l l" Psi - constant surfaces in transformed coordinates‘l l‘ Psi = 0 ‘I forfi=Oj linclude linclude #include linclude linclude "engineh" l" Plot the shape of the object and additionally provide the color map using dependent variable values if available ‘/ plot_surflpxpyp,zp,c,n,m,l) int p; dOUb‘e ..‘xp'ttfiyp'COQZP’Ottc; int mm]; 1 Matrix ‘x.‘y,‘z,‘d."df; double ‘dr; extern Engine ‘eptr; double ‘rl.‘img; int iJ.k.s. p 1.total; double di.dj.dk,x_min.y_min.2,min,x_,max.y_max.z_max; double ‘tempc.‘tempx.‘tempy.‘tempz.find_max().find_min(): char ‘erg.str[60].str1[ 501.con1 501.err1 50].mat_comml 200].fi g_name{75]; totalzn‘m‘l; if((p==1)l l(p-=-=2)) 1 di=(double) n; dj=(double) m; dk=(double) l; tempz=(double ‘)malloc(total‘sizeofldouble»; tempx==(double ‘)malloc(total‘sizeofldouble)); tempy=(double ‘)malloc(total‘sizeof(double)); mat_set_fiill(xp.tempx.n.m.l); mat_set_firll(yp.tempy,n,m.l 1; mat_set_full(zp.tempz.n.m.l); x_min=find_min(xp.n.m.l); y_min=find_min(yp,n.m.l): z-min=find_min(zp.n.m.l); x.max=find_max(xp.n.m.l): y_max=find_max(yp.n.m.l); z_rnax=find_max(zp.n.m,l); engPutI-‘ull(eptr."n".1.l.&di.NU LL); engPutFull(eptr."m”. 1. 1.&dj.NULL); engPutFull(eptr.”l".1.l.&dk.NU LL); engPutFull(eptr."xmx".1.1.&x_max.NULL); engPutFull(eptr."ymx".l.l.&y_max.NULL); engPutFull(eptr."zmx”.1.l.&z_max.NU LL); engPutF‘ull(eptr.'xmn".1.1.&x_min.NULL); engPutFull(eptr.”ymn".1.1.&y_min.NULL); engPuthill(eptr."zmn”.l.1.&z_min.NULL); engPutFull(eptr."x”.total.1.tempx.NULL); engPutI-‘ull(eptr."y".total.1.tempy.NULL); engPutFull(eptr.”z".total.1.tempz.NULL); 1 iflp== 1) engEvalStri ng(eptr."figu re( 1);"); i111p==1)| l(p==2)) 1 iflpas2) engEvalStrIng(eptr."figure(find_pic(°Shape Modelling')):"): sprintflcon."subplot( 1.2.%d);”. p); engEvalString(eptr.con); engEvalString(eptr."hold off ;"); engEvalString(eptr."shape_plot(x.y.z.n.m.l.1);"); ifl -=1) 1 engEvalString(eptr."title('Initial Guess');"): sprint“mat_comm.”set(gcf.'name'.'Shape Modelling');"); engEvalString(eptr.mat_comm): 1 else engEvalString(eptr."title('Boundary Fitted Object’):”): free((char ‘)tempx): free((char ‘)tempy): free((char ")tempz): 1 else 142 1 tempe=(double ‘)malloc(n‘m‘l'sizeofidoubleD: mat_set_firll(c.tempc.n.m.l); engI’utIfirll(eptr."c".n‘m‘l.1.tempc.NULL); engEvalString(eptr."shape.plot(x.y.z.n.m,l.1.c);"); engEvalString(eptr."shading interp;"); free((char l")tempc): 1 return; 1 P Set up the three dimensional structures into matlab vectors for plotting ‘/ mat_set_1i.ill(p.temp.n.m.l) double ”“p; double ‘temp; int n,m,1; 1 int r.iJ.k; r20; for(i=0;i : "): scanfl"%d".&ch); while(ch1=4) 1 if(count== 1) 1 engEvalString(eptr."figure;"); engEvalString(eptr."hi=gcf ;"1; you want to view ifics=(double “‘)NU LL) 1 if(b== 1) sprintflmat_comm."set(gcf.'name'.'&r rf ace Plots-1');"): else sprintflmat.comm."set(gcf.'name'.'Su rface Plots-2');"): 1 else 1 iflb==1) sprintfimat_comm."set(gcf.'name'.'Surface Temperature Map');"): else sprint“mat_comm.”set(gcf.'name'.'Surface Map'):"): 1 engEvalString(eptr.mat_comm); count-H; 1 switch(ch) 1 case 1:1 printfl”Psi - "): scanfl"%d".&sur0: engEvalString(eptr."figure(hi):"): engEvalStr-ing(eptr."subplot( 1.3.1 ):"): engEvalStri ng(eptr."hold 011;”); plot_su rface( 1.x.y.z.e.n.m.l.sur0; engEvalString(eptr.”view(90.20):"): engEvalString1eptr.”axis om"); sprintfipsi_ar."title('u = %d plane');".surf); engEvaIString(eptr.psi-ar); engEvalString(eptr.”end;"); engEvalString(eptr."pause;"): cnti++; break; 1 case 2: 1 printf("Eta = "): scanfi”%d”.&sur1): engEvalString(eptr."figur8(hi):"): engEvaIString(eptr."subplot( 1.3,2):"): engEvalStr-ing(eptr."hold 011;"); plot_surface(2.x.y.z.c.n.m.l.sur1): engEvalString(eptr.”view(100.30);"): engEvalString(eptr."axis 011;"): sprintfieta_ar."title('v = %d plane');".surfi; engEvalString(eptr,eta_ar); engEvalString(eptr.”end:"): engEvalString(eptr."pause:"): cntj++z break; 1 case 3: 1 printfl”Zei = ”1: scanfl"%d”.&surfi: engEvalString : "): scanfl"%d”.&ch): 1 mat_lab_options(): return; 1 you want to view I‘ Plot 8 single surface/coordinate-section of the shape ‘/ plot_surfacdp.x.y.z.c,n.m.l.i) I‘ To plot the section requested ‘/ int p; /" Psi. Eta or Zei surface index (1. 2 or 3 resptly) ‘/ double Temperature matrices ‘/ int n.m.l.i; /‘ Size of plane and corresponding index for 3rd coordinate ‘/ 1 extern Engine ‘eptr; /" Matlab Workspace ‘/ double dp.di.‘tempx.‘tempy.‘tempz.‘tempc; int j.k; dp=(double)p; di=(double)(i+1): engPutIfirll(eptr.”index".1.1.&dp.NULL): engPutFull(eptr."i".1.1.&di.NULL); if(c==(double ""l NULL) engEvalString1eptr."surface_plot(x.y.z.n.m.l.i.index);" ): else 1 tempc:(double "')malloc(n‘m‘l‘sizeofidouble)1: mat_set_full(c.tempc,n.m.l); engPutFull(eptr."c",n‘m‘l,1.tempc.NULL); engEvalString(eptr."surface_plot(x.y.z.n.m.l,i.index,c) z"): engEvalString(eptr,”shading interp;"); free((char ‘)tempc): 1 engEvalString(eptr.”clear c;”); return; 1 cssx'csoy,oooz.osoc; [0 and Shape /" Matlab options routine. Acts as conduit between the user input at terminal to the matlab software workspace ‘/ mat_lab_options() 1 extern Engine ‘eptr; char yn13l.‘err; err=(char ‘)malloc(5‘sizeofichar)): printf("Do you want to execute Matlab commands (y/n) :"): scanfl"%s".err); while((stremp(err."quit")1=0)&&(strcmp(err."n")!=0)& &(stremp(err.”quit;")!=0)) 1 free((char ‘)err): em(chu ')malloc(75‘sizeof1char)): printf(”Command : ")z scanfi"%s".err): ifl(strcmp(err."quit")1=0)&&(strcmp(err."quit;")!=0)) 1 engEvalString(eptr.err); engEvalString(eptr.”figure(gcf);"); 1 1 free((char ‘)err): return; 1 print_plot_options() 1 char yn13].‘mat_comm1; int fig_no; mat_comm1=(char ‘)malloc(60’sizeofichar)); printf("Do you want any plots to be saved ? :"): scanfi”%s".yn); while(strcmp(yn."y”)==0) 1 printf("Figure No. for which hard copy is required ? ”)z scanf("%d".&fig_no); sprintflmat_comm1."figure(%d);print Figure%d”.fig_no.fig_no); engEvalString(eptr.mat_comm1); printf("Do you want any other figure to be saved ? :"): scanfl"%s”.yn); 1 free( mat.comm 1 ); return; 1 WW #include flinclude linclude linclude linclude ”grid.h" Oinclude ”engine.h" P Graphical display of the required variables as 3D plots along two directions. color maps indicating trends and time plots at different points as specified in the 'init_display.m' matlab file ‘/ Display_variables(t.grid_pt.T.X.n.m.l) double t; /‘ Time ‘/ struct point ‘grid_pt; l" Grid point structure containing physical coordinate information ‘I double “‘T."”X; l‘ Temperature and Cure variables ‘/ int n,m,1; P Mesh size ‘/ 1 extern int psi.eta.zei: int no_fig.cnt; extern char plot_choicel 10]; extern char map151.ttx151: 144 extern Engine ‘eptr; extern double t_display; double ‘tempT.’teme.‘"x."‘y."‘z; “11:0; 111 psil=o 1) cnt++; ifleta1=-1) cnt++; iflzei !=- 1) cnt++; no_fig=cnt; /’ Plot Temperature maps and profiles ‘/ i11(no_fig1=0)&&((strcmp(plot_choice."l‘emp")==0)l I ( strcmp(plot_choice."Both")==0))) 1 if(stremp(map."yes")==0) 1 ifit<=1.0eo8) 1 engEvalString1eptr.”figure;set(gcf.'name'."l‘emperatu re Profiles');"); plot _graph(t.no_fig,grid_pt.T.n.m.l): engEvalString(eptr."figure;set(gcf.'name'.'Temperatu re Mash"): plot_map(t.grid_pt.T.n.m.l.1); 1 else 1 engEvalString(eptr."figurc(find_pic('Temperature Profiles'));"): plot_graph(t.no_fig.grid_pt.T.n.m.l): engEvalString(eptr."figure(find_pic('l‘emperature Map')):"): plot_map(t.grid_pt.T,n.m.l.1); 1 l /" Update time-temperature plot ‘/ iflstrcmp(ttx."yes")==0) update_time_plot(t.T.n.m.l."T"); l /‘ Plot Cure maps and profiles */ ifi(no_fig1=0)&&((strcmp(plot_choice."Cure")==0) I l (s trcmp(plot_choice."Both")==0))) 1 iflstrcmp(map.”yes")==0) 1 ifit<=1.0e-8) 1 engEvalString(eptr.”figure(gcf+ 1 );set(gcf.'name‘.'Cure Profiles');"): plot _graph(t.no_fig.grid_pt.X.n.m.l): engEvalString(eptr."figure;set(gcf.'name'.'Cure Map'):"): plot_map(t.grid_pt.X.n.m.l.2): 1 else 1 engEvalString(eptr.”figure(find__pic('Cure Profiles'));"); plot _graph(t.no_fig.grid_pt.X.n.m,l); engEvalString(eptr."figure(find_pic('Cure Map'));"); plot_map(t.grid_pt.X.n.m.l.2); 1 1 /‘ Update time-cure plots ‘/ ifistrcmp(ttx.”yes”)==0) update_time_plot(t.X.n.m.l."X"): 1 return: 1 /‘ Graphical output of Temperature or Cure Vs (Psi.Eta).(Eta.Zei) and (Zei.Psi) directions ‘/ plot_graph(t.no_fig.grid_pt.f.n.m.l) double t; int no_fig; struct point ‘grid_pt: double “‘1‘; int n.m.l; 1 extern int psi.eta.zei; extern Engine ‘eptr; double ‘temp; char stri501: int}; j=1: if((psi1=-l)) 1 sprintflstr."subplot(1.%d.%d);”.no_fig.j); engEvalString(eptr.str); plot_profitJ.psi.f.n.m.l); engEvalString(eptr."pause;"): j++; 1 metals-1) 1 sprintflstr."subplot(1.%d,%d);".no_fig.j); engEvalString(eptr.str); plot_profitj,eta.f.n.m.l); engEvalString(eptr.”pause;"): j++; 1 ifizei1=-1) 1 sprintfistr.”subplot( 1.%d.%d);".no_figJ ): engEvalString(eptr.str); plot_profltj.zei.f.n,m.l); engEvalString(eptr."pause:"); 1 return; 1 plot_proflt.ch.i,f.n.m,l) double t; int ch.i; double ""f; int n.m.l; 1 extern Engine ‘eptr; int sJ.k; char stri501.con'11501: double 'temp.‘”mem_alloc_3D().di; di=(double) i; temp=(double ")mallodn‘m‘l‘sizeofidouble»: mat_set_fi.rll(f.temp.n.m.l); engPutNll(eptr.”P”.n‘m’l.1.temp.NU LL); I45 engPutFull(eptr."I".1.1.&di.NULL): switch(ch) 1 case 1: 1 engEvalString(eptr."plot_xy(P.1.I.n.m.l);"); sprintf(str."title('Psi=%d;Time:%6.llf units');”.i,t); engEvalString(eptr.”xlabel('Eta');ylabel('Zei');zlabel('T '):"): break; 1 case 2: 1 engEvalString(eptr."plot_xy(P.2.l,n,m,l):"): sprintflstr.”title('Eta=%d;Time:%6. llf units'):".i.t); engEvalString(eptr."xlabel('Psi');ylabel('Zei');zlabe1('T' 1:”); break; 1 case 3: 1 engEvalStri ng(eptr."plot_xy(P.3.I.n.m.l);"); sprintflstr."titlc('Zei=%d;Time:%6.11funits'):".i.t); engEvalString(eptr.”xlabel('Psi');ylabel('Eta');zlabel('T '1:"): break; 1 1 engEvalString(eptr.str); engEvalString(eptr.”clear P;"); free((char *)tcmp); return; 1 l‘ Plot the color map of Temperature and or Cure for entire shape ‘/ plot_map(t.grid_pt.T.n.m.l.b) double t; struct point ‘grid_pt; double "‘1‘; int n.m.l.b; 1 extern Engine ‘eptr; char str1100]; engEvalString(eptr."hold off;"); engEvalString(eptr.”set(gca.'Visible°.'off);"); engEvalString(eptr."view(3);"); engEvalStri ng(eptr.”hold on;”); sprintfistr."title(Temperature Map: units);".t); engEvalString(eptr.str); plot_surf(3.grid_pt->x.grid_pt->y.grid_pt->z.T.n.m.l); engEvalString(eptr."hold om"); Time=%6. llf l‘ view_options(grid_pt->x.grid__pt->y.grid_pt. >z.T.n,m.l.b); ‘/ return; 1 I‘ Update the time-temperature and/or time-cure plots ‘/ update_time_plct(t.f.n.m.l.str) double t.‘"f; int n,m,1; char ‘str; 1 double ‘temp; extern Engine ‘eptr; temp=(double ‘)mallodn‘m‘l‘sizeofidouble»: mat_set_full(f,temp.n.m,l); engPutFull(eptr."t”.1.1.&t.NULL); engPutFull(eptr."c".n‘m‘l.1,temp.NULL); ifit<=1.0e-8) 1 engEvalString(eptr."figure:"); iflstremp(str.”'l‘")==0) 1 engEvalString(eptr."set(gcf.'name'.'Time- Temperature Plots');"): engEvalString(eptr.”time_tplots(t,c.n.m.l);"): engEvalString(eptr.”title("l‘emperature Vs Ti me');"): engEvalString(eptr.”xlabel("I‘ime ---->'):"): engEvalString(eptr."ylabel(‘Temperature ----- >');"): 1 else 1 engEvalString(eptr.”set(gcf.'name'.'Time.Cure Plots');”); engEvalString(eptr.”time_xplots(t.c.n.m.l):"): engEvalString(eptr.”title('Cure Vs Time');"); engEvalString(eptr.”xlabel('Time -~-->');"); engEvalString(eptr.”ylabel('Cure ----- >'):"): engEvalString(eptr."axis(10 t 0 11):"): 1 1 else 1 ifistrcmp(str."'l'")==0) 1 engEvalString(eptr.”figure(find_pic('Time- Temperature Plots'));"); engEvalString(eptr."time_tplots(t.c.n.m.l):"): engEvalString(eptr."axis(10 t Tin max(c)1):"): 1 else 1 engEvalString(eptr.”figure(find_pic("l‘ime-Cure Plots')):"): engEvalString(eptr."time_xplots(t.c.n.m.l):"): engEvalString(eptr."axis(10 t 0 11):"): 1 1 engEvalString(eptr."clear cz"); free(temp); return; 1 initialize_workspace(grid_pt,n.m.l) struct point ‘grid__pt; int n.m.l; 1 extern Engine ‘eptr; extern double final_time.t_display; extern double T.initial,X_initial; double ‘temp; engEvalString(eptr.”init_display;"); temp=(double ‘)malloc(n‘m'l’sizeofidouble»: mat_set_full(grid_pt->x.temp.n.m.l); engPutNll(eptr.”x".n‘m‘l.1.temp.NULL); I46 mat_set_full(grid_pt->y.temp.n.m.l); engPuthJll(eptr."y".n‘m“l.1.temp.NU LL); mat_set_full(grid_pt->z.temp.n.m.l); engPutFull(eptr."z”.n‘m‘l.1.temp.NU LL); engPutNll(eptr.”tfinal".1.1.&final_time.NULL); engPutFull(eptr.”tdisp".1.1.&t_display.NULL): engPutFull(eptr."'I‘in".1.1.&T_initial.NU LL): engPutFull(eptr."Xin".1.1.&X_initial.NULL); free(temp); return; 1 The following set of programs contain general purpose routines on memory management. matrix operations and the various derivative operations performed numerically in space w.r.t natural coordinates. We #include #include #include #include "grid.h" /"' 'I‘ransfinite Interpolation Routine */ transfinite_30(f,n,m.l) double "*f; int n,m,1; 1 double "*F1,"*mem_alloc_3[)(),int_plt(); double ”*F‘2; double "*FB; int iJ,k; F1=mem_alloc__3[)(l,n,m); /"‘ Correction along j(eta) direction */ for(k=0;k #include mat_add(c.a,b.n,m,l) double ttsc’sssa,**sb; int n,m,1; 1 int ickk; for(i=0;i=fli1[j 111(1) ret=fli 11.111111; 1 1 1 return(ret); 1 double norm(f,n,m,l) double "*f; int n,m,1; 1 double ret,find_abs_rnax(); int i,j,k; ret=find_abs_max(f,n,m,l); return(ret); 1 PM atrix Solution*/ matrix_solve(a,m) double "a; int m; 1 int iJ,k,l,r; double temp,t,s; r=0; while(r1.0e-4) 1 for(l=0;l<=m;l++) 1 temp=a[r][ll; 31r1111=81k1111; 81k1111=temp; 1 l 1 1 for(k=0;k #include double f_psi(f,iJ,k,n.m,l,bc) /* Psi derivative */ double "’"f; int i,j,k,n,m,l; int be; 1 extern int N; double fp,h; h4doubleXN-l.0)l(n-1.0); if1(i>0)&&(i0)&&1i0)&&(k