// Simple 1D Mass Transfer (parameters a function of moisture) // Type: Fickian Diffusion // Shape: Slab // BC: Fixed Moisture (infinite moisture sinks, "Dirichlet") // Params: Variable! Diffusivity & density can be a function of moisture, time, and position!!! // Method: Explicit (small memory version) // By: Wayne Pafko // Created: June 25th, 2001 // Updated: June 26th, 2001 #include // The standard C library #include // How fast is this program? #define SUCCESS 1 #define FAILURE 0 #define SMALL 1e-6 // A small number #define ASCII_SPACE 32 // ASCII decimal value for a space #define ASCII_TAB 9 // ASCII decimal value for a tab // Allows the use of "excel" style formulas to define D & DENSITY (just cut and paste from excel to formula's below) #define IF(A,B,C) ((A)?(B):(C)) // Solution specific definitions (change these & recompile) #define NODES 10 // How many nodes do we have (no decimal point please) #define STEPS 1000 // How many time steps do we have (no decimal point please) #define TIME 30.0 // (sec) Total time #define WIDTH 0.4 // (cm) Total width of simulation // Functional relationships (change these & recompile) Excel friendly (paste back and forth)!!! // In general, each parameter is a function of concentration - c_100 (% of total), time - t_s (s), and position - x_cm (cm) // Concentration range: 0 < c < 100 // Time range: 0 s < t < "TIME" s // Position range: 0 cm < x < "WIDTH" cm. (Center node is at "WIDTH"/2 * cm) #define D(c_100, t_s, x_cm) 3E-5 // for Potato from Corrigan's Model // (cm2/s) Diffusivity #define DENSITY(c_100,t_s,x_cm) 0.8 // (g/cm3) Density // note: We use macros instead of "functions" because they are much faster! // Important! Leave extra parenthesis around the individual variables, otherwise complex expressions won't be inserted into equation correctly!!! // Problem specific definitions (change these & recompile) #define INITIAL_CONC 40.0 // (% of total mass) Initial concentration across slab #define BC_LEFT 2.0 // (% of total mass) Left BC (0.0 to 100.0) #define BC_RIGHT 2.0 // (% of total mass) Right BC (0.0 to 100.0) #define K_LEFT 0.5 // (dimless) Partition Coefficient on left side of slab (0.0 to 1.0) #define K_RIGHT 0.5 // (dimless) Partition Coefficient on right side of slab (0.0 to 1.0) //#define C_ADSORB 1.0 // (% of total mass) Adsorbed water (tightly bound, does not participate in mass transfer) double Cmin_predict = 0.0; // (% of total mass) Predicted minimum concentration. Should solution converge? Also effects parameter output range. double Cmax_predict = 100.0; // (% of total mass) Predicted maximum concentration. Should solution converge? Also effects parameter output range. // Unit conversion factors #define FOURIER_FACTOR 1 // [dimless: (cm2/s) * (s) / (cm2)] Conversion constant for equation k/cp/density*dt/dx/dx //#define M_OVER_CM 0.01 // [dimless: (1 m) / (100 cm)] Conversion constant from meters to cm // Parameters that effect how the results are saved (but not the solution method) #define FILE_NAME "data.txt" // Save results in this file #define PRINTSTEP 10 // How many time steps before a result is printed (if "1" see data at every time step, if "2" see every other step) #define PRINTNODE 1 // How many nodes before a result is printed (if "1" see all the nodes, if "2" see half) #define PRINTNODESTART 1 // What node should I start printing at? (usually 1, always >=1) // Derived values based on above definitions double dt = (double)TIME/STEPS; // (sec) Time step double dx = (double)WIDTH/NODES; // (cm) Distance between nodes double dtdx2 = (double)TIME/STEPS * NODES/WIDTH * NODES/WIDTH; // So don't have to keep recalculating (dt/dx/dx) at each time step // Memory stuff double storage_a[NODES+1]; // Storage for algorithm; "NODES+1" for safety reasons (don't want to fly past alotted space accidently) double storage_b[NODES+1]; // More storage for algoritm // Other useful numbers that will be calculated (leave these alone) double * Conc_now; // (% of total mass) Current concentration of moisture (this time step) (will point to storage_a or storage_b) double * Conc_past; // (% of total mass) Past concentration of moisture (one time step ago) (will point to storage_a or storage_b) double aveConc[STEPS]; // (% of total mass) Average concentration of moisture in the slab double fourierDtDx2; // (dimless) Fourier number for a single step & node (must be small for convergence) double CBC_left[STEPS]; // (% of total mass) Left BC (may change with time) double CBC_right[STEPS]; // (% of total mass) Right BC (may change with time) //double Cmin; // (% of total mass) Minimum concentration found during solution (Same as CMIN_PREDICT?) //double Cmax; // (% of total mass) Maximum concentration found during solution (Same as CMAX_PREDICT?) // Time the algorithm time_t startTime, endTime; // "time_t" structure defined in time.h // Sets initial conditions and boundary conditions int initial_conditions(void) { long i, t; // Loop variable // Set up pointers (to storage space) // Establish where Temp_now & Temp_past point to... Conc_now = & storage_a[0]; Conc_past = & storage_b[0]; // Check expected memory & file usage printf("Expected sizes: memory %.2f mb, file %.2f mb\n",1.1*(sizeof(storage_a)+sizeof(storage_b))/1E6, NODES/PRINTNODE*STEPS/PRINTSTEP*8*8*2/1E6); // Initial Conditions for interior nodes for(i=0;iconverge) printf("\t\tbut it may oscillate if above %.3f!\n",converge); printf("\t\tnote: Truncation error minimized below %.3f\n",minError); } else { printf("\t\tERROR! The solution will not converge!\n\t\tIncrease # of time steps (STEPS) and recompile...\n"); return(FAILURE); } // Print material properties printf("\n---Material Properties---\n"); printf(" Density:\t%.4f g/cm3\n",(double)DENSITY(Cmax_predict,(float)TIME * dt,0.0)); printf(" Diffusivity:\t%.6g cm2/s\n",diffusivity); //SAVE TO FILE STUFF///// printf("\nDoing math and saving in %s...\t ",FILE_NAME); // Open the file (check for failure) if( (dataFile = fopen(FILE_NAME,"w"))==FAILURE ) { printf("[failure]\nIs %s open in another program?\n",FILE_NAME); return(FAILURE); } // Start printing to the file fprintf(dataFile, "1D Mass Transfer\nby Wayne Pafko\n"); //Print solution properties fprintf(dataFile, "\nNodal Representation\n"); fprintf(dataFile, "Nodes:\t%d\n",(long)NODES); fprintf(dataFile, "Slab Width:\t%.4g\tcm\n",(double)WIDTH); fprintf(dataFile, "Spacing:\t%.6g\tcm\n",dx); fprintf(dataFile, "Steps:\t%d\n",(long)STEPS); fprintf(dataFile, "Total Time:\t%.2f\ts\n",(double)TIME); fprintf(dataFile, "Time Step:\t%.6g\ts\n",dt); // Calculate & print material properties fprintf(dataFile, "\nMaterial Properties\n"); // Print table of density fprintf(dataFile, "Density (g/cm3)\n"); fprintf(dataFile, "(%%)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx); for(Cloop=Cmin_predict;Cloop<=Cmax_predict+SMALL;Cloop+=(Cmax_predict-Cmin_predict)/4.0) { fprintf(dataFile, "%.2f\t", Cloop); for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx) fprintf(dataFile, "%.4f\t", DENSITY(Cloop,0.0,xloop) ); fprintf(dataFile, "\n"); } fprintf(dataFile, "%.2f\t", Cmax_predict); for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx) fprintf(dataFile, "%.4f\t", DENSITY(Cmax_predict,(double)TIME,xloop) ); fprintf(dataFile, "%.2f (s)\n", (double)TIME); // Print table of diffusivity fprintf(dataFile, "Diffusivity (cm2/s)\n"); fprintf(dataFile, "(%%)\t%.4f\t%.4f\t%.4f\t:(cm)\n",0.5*dx,(float)WIDTH/2.0,(float)WIDTH-0.5*dx); for(Cloop=Cmin_predict;Cloop<=Cmax_predict+SMALL;Cloop+=(Cmax_predict-Cmin_predict)/4.0) { fprintf(dataFile, "%.2f\t", Cloop); for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx) { diffusivity=D(Cloop,0.0,xloop); fprintf(dataFile, "%.6g\t", diffusivity ); } fprintf(dataFile, "\n"); } fprintf(dataFile, "%.2f\t", Cmax_predict); for(xloop=0.5*dx; xloop<=((double)WIDTH-0.5*dx+SMALL); xloop+=(float)WIDTH/2.0-0.5*dx) { diffusivity=D(Cmax_predict,(double)TIME,xloop); fprintf(dataFile, "%.6g\t", diffusivity ); } fprintf(dataFile, "%.2f (s)\n", (double)TIME); // Print Boundary condition information fprintf(dataFile, "\nBoundary Conditions\n"); fprintf(dataFile, "\t[left]\t[right]\n"); fprintf(dataFile, "Partition Coefficient:\t%.4f\t%.4f\n",(double)K_LEFT,(double)K_RIGHT); fprintf(dataFile, "Exterior Moisture:\t%.4f\t%.4f\t(%%)\n",(double)BC_LEFT,(double)BC_RIGHT); // Print headings for main temperature table fprintf(dataFile, "\nConcentration Results (%% of total mass)\n"); fprintf(dataFile, "\t\tx (cm):\t%.4f\t",0.0); for(i=PRINTNODESTART-1;i to continue..."); do{ presskey=getchar(); //putchar(presskey); } while (presskey!='\n'); //scanf("%c",&presskey); //Another way of doing the same thing } // The "main" program (start here) int main (void) { printf("1D Mass Transfer\nby Wayne Pafko\n\n"); // Run appropriate subroutines (if they fail, don't continue) if( initial_conditions() !=SUCCESS) {press_enter(); return(FAILURE);} if ( solve_and_save_problem() !=SUCCESS) {press_enter(); return(FAILURE);} //if( save_results() !=SUCCESS) {press_enter(); return(FAILURE);} append_source_file(); // No big deal if it fails (couldn't find the source file) // All done, ready to leave press_enter(); return(SUCCESS); } // Created and compiled with Microsoft Visual C++ 4.0