This shows you the differences between two versions of the page.
| Previous revisionNext revision | |||
| — | tutorials:t23 [2019/06/25 15:20] – [Complementary files] pwarczok | ||
|---|---|---|---|
| Line 1: | Line 1: | ||
| + | ===== T23: Introduction to cell simulations: | ||
| + | //This tutorial was tested on\\ | ||
| + | MatCalc version 6.01 rel 1.003\\ | ||
| + | license: free\\ | ||
| + | database: mc_fe.tdb; mc_fe.ddb// | ||
| + | |||
| + | ==== Complementary files ==== | ||
| + | |||
| + | Click {{: | ||
| + | |||
| + | ==== Contents ==== | ||
| + | |||
| + | * setting up a simulation grid | ||
| + | * assigning properties to cells | ||
| + | * setting up display of cell simulation results | ||
| + | * setting up simulation parameters | ||
| + | * boundary conditions | ||
| + | |||
| + | MatCalc can be used to carry out simulations of long-range diffusion problems, and problems in which diffusion is coupled with precipitation, | ||
| + | |||
| + | The functionality for setting up and running diffusion simulations is only available through the scripting language and not through the GUI. This tutorial therefore also serves as a refresher for some of the basic scripting commands that we have already encountered, | ||
| + | |||
| + | ===== Description of the problem ===== | ||
| + | |||
| + | Fick's Second Law of Diffusion is a partial differential equation of the form:\\ | ||
| + | |||
| + | \[ \left({\partial C \over \partial t}\right) = D . \left({\partial^2 C \over \partial x^2}\right) \] | ||
| + | |||
| + | where $D$ is the diffusion coefficient, | ||
| + | |||
| + | One of these is the thin-film solution, in which an very thin layer containing a finite amount of the diffusing species is sandwiched between two semi-infinite bars of material. The boundary conditions are that: | ||
| + | |||
| + | * initially, the concentration of the solute in the bar is zero; this is expressed mathematically as $ C\{x=0, | ||
| + | |||
| + | * the total amount of solute is finite and constant throughout the time period considered; this is expressed mathematically as | ||
| + | |||
| + | \[\int_{-\infty}^{\infty} C \{x,t\} \mathrm{d}x = B,\] where $B$ represents the amount of solute in the system. | ||
| + | |||
| + | {{: | ||
| + | |||
| + | The standard solution when the boundary conditions are of this type is of Gaussian form: | ||
| + | |||
| + | \[ C\{x,t\} = {B \over 2 \sqrt{\pi D t}} . \mathrm{exp}\left({-x^2 | ||
| + | |||
| + | In this tutorial, we will set up a cell-based numerical simulation of a situation that fulfils these conditions to a reasonably good approximation, | ||
| + | ===== Setting up the system ===== | ||
| + | |||
| + | |||
| + | Start MatCalc and create a new script file as explained in [[tutorials: | ||
| + | ==== a. Thermodynamic and kinetic setup ==== | ||
| + | |||
| + | The first task that the script must perform is to create a new MatCalc workspace. Type the following command into your script window: | ||
| + | |||
| + | < | ||
| + | use-module core | ||
| + | </ | ||
| + | |||
| + | This tells MatCalc to use the //core module//, which it selects at the begining by default anyway; later in this tutorial we will see the use of the // | ||
| + | |||
| + | < | ||
| + | new-workspace | ||
| + | </ | ||
| + | |||
| + | This creates a new workspace. If an existing unsaved workspace is open when the script is run, you will be prompted to save it. To disable this prompt, append '' | ||
| + | |||
| + | Some workspace information can be added as shown: | ||
| + | |||
| + | < | ||
| + | @$************************************ | ||
| + | $ enter workspace info | ||
| + | @$************************************ | ||
| + | @ set-workspace-info Script T23 | ||
| + | @ set-workspace-info +Simulation of long-range diffusion of C | ||
| + | @ set-workspace-info +out of thin film | ||
| + | @ set-workspace-info +with matrix phase FCC_A1. | ||
| + | </ | ||
| + | |||
| + | The next step is to set up the thermodynamic and kinetic data. Enter the following commands into your script window and save the script. | ||
| + | |||
| + | < | ||
| + | open-thermodyn-database mc_fe.tdb | ||
| + | select-elements fe c | ||
| + | select-phases FCC_A1 | ||
| + | read-thermodyn-database | ||
| + | read-mobility-database mc_fe.ddb | ||
| + | </ | ||
| + | |||
| + | Following this, define the variables that will specify the initial carbon content in the bar ('' | ||
| + | |||
| + | < | ||
| + | set-variable-value c_comp_bar_wp 1e-10 | ||
| + | set-variable-value c_comp_layer_wp 0.8 | ||
| + | </ | ||
| + | |||
| + | In the cell simulation, compositions must be entered in molar notation '' | ||
| + | |||
| + | < | ||
| + | enter-composition wp c=c_comp_bar_wp | ||
| + | set-variable-value c_comp_bar un$c | ||
| + | enter-composition wp c=c_comp_layer_wp | ||
| + | set-variable-value c_comp_layer un$c | ||
| + | </ | ||
| + | |||
| + | Then set a temperature and calculate an equilibrium to initialise the system and check that everything is working OK. | ||
| + | |||
| + | < | ||
| + | set-temperature-celsius 1000 | ||
| + | calculate-equilibrium | ||
| + | </ | ||
| + | |||
| + | At this stage, also specify the simulation time (in seconds) for subsequent use: | ||
| + | |||
| + | < | ||
| + | set-variable-value sim_time 1e7 | ||
| + | </ | ||
| + | |||
| + | Save the script. It is recommended to test it out at this stage to catch any typing errors before going further. | ||
| + | |||
| + | ==== b. Precipitation domain setup ==== | ||
| + | |||
| + | Before setting up the cell simulation, you must define a // | ||
| + | |||
| + | < | ||
| + | create-precipitation-domain aus_domain | ||
| + | set-precipitation-parameter aus_domain X FCC_A1 | ||
| + | </ | ||
| + | |||
| + | |||
| + | |||
| + | ==== c. Creating the simulation grid ==== | ||
| + | |||
| + | The next step is to set up the cell simulation. First change to the // | ||
| + | |||
| + | < | ||
| + | use-module simulation | ||
| + | </ | ||
| + | |||
| + | Define a variable representing the number of cells: | ||
| + | |||
| + | < | ||
| + | set-variable-value num_cells_x 99 | ||
| + | </ | ||
| + | |||
| + | Then create a simulation grid using the command: | ||
| + | |||
| + | < | ||
| + | create-simulation-grid num_cells_x 1 1 | ||
| + | </ | ||
| + | |||
| + | The three numbers here represent the number of cells in the x-, y- and z-directions, | ||
| + | |||
| + | < | ||
| + | set-variable-value barlength 0.099 | ||
| + | </ | ||
| + | |||
| + | Specify the dimensions in these directions (in m) using the command: | ||
| + | |||
| + | < | ||
| + | set-grid-coordinates barlength 1 1 | ||
| + | </ | ||
| + | |||
| + | and the geometry of the grid, using the command: | ||
| + | |||
| + | < | ||
| + | set-grid-geometry p | ||
| + | </ | ||
| + | |||
| + | The three geometry options are: | ||
| + | * p - planar | ||
| + | * c - cylindrical | ||
| + | * s - spherical | ||
| + | |||
| + | To visualise this grid, create a new window of the type g5 (Paint grid: 2D cells): | ||
| + | |||
| + | < | ||
| + | new-gui-window g5 | ||
| + | </ | ||
| + | |||
| + | Give this window a name so that it can be referred to easily later. The variable '' | ||
| + | |||
| + | < | ||
| + | set-variable-value paint_window_id active_frame_id | ||
| + | </ | ||
| + | |||
| + | Modify the position, size and zoom level using the following two commands: | ||
| + | |||
| + | < | ||
| + | move-gui-window paint_window_id 0 100 750 160 $ change position and size | ||
| + | set-gui-window-property paint_window_id i Z 200 $ zoom 200% | ||
| + | </ | ||
| + | |||
| + | The grid should now look as shown: | ||
| + | |||
| + | {{: | ||
| + | |||
| + | In the next section, properties will be assigned to the cells in the grid. | ||
| + | |||
| + | ==== d. Assigning properties to the cells ==== | ||
| + | |||
| + | For the diffusion simulation, it is necessary to create a // | ||
| + | |||
| + | < | ||
| + | create-material aus_material | ||
| + | set-material-property aus_material G D aus_domain | ||
| + | $ [g]eneral properties: [d]omain | ||
| + | </ | ||
| + | |||
| + | Then set the diffusion coefficient. Here, a user-specified value is used, instead of one calculated using information from the database, to better facilitate comparison with the analytical solution. The diffusion coefficient is stored in a variable for later use. The syntax below means [d]iffusion [c]oefficient of [c]arbon, [f]unction or expression. | ||
| + | |||
| + | < | ||
| + | set-variable-value diffCoeff 1e-12 | ||
| + | set-material-property aus_material D C C F diffCoeff | ||
| + | $[d]iffusion [c]oefficient of [c]arbon, [f]unction or expression | ||
| + | </ | ||
| + | |||
| + | Following this, assign properties to the cells. Using the syntax below, the material '' | ||
| + | |||
| + | < | ||
| + | set-cell-property a m aus_material $[a]ll cells, [m]aterial | ||
| + | </ | ||
| + | |||
| + | Set the temperature of all cells to 1000˚C, by first specifying that the temperatures should be in ˚C: | ||
| + | |||
| + | < | ||
| + | set-simulation-parameter N Y | ||
| + | </ | ||
| + | |||
| + | specifying the temperature in the form of a variable: | ||
| + | |||
| + | < | ||
| + | set-variable-value sim_temp 1000 | ||
| + | </ | ||
| + | |||
| + | and then attaching this temperature property to all cells: | ||
| + | |||
| + | < | ||
| + | set-cell-property a V T sim_temp $ [a]ll cells, [v]ariable, [t]emperature | ||
| + | </ | ||
| + | |||
| + | Assign the lower carbon content, '' | ||
| + | |||
| + | < | ||
| + | set-cell-property a V C C c_comp_bar | ||
| + | $ [a]ll cells, [v]ariable, [c]omposition [c]arbon | ||
| + | </ | ||
| + | |||
| + | The next stage is to select the cell in the middle of the bar using its cell index. First, we note that the first cell has the index 0, so the highest cell index in the array obtained using: | ||
| + | |||
| + | < | ||
| + | set-variable-value max_cell_index num_cells_x-1 | ||
| + | </ | ||
| + | |||
| + | and the index of the cell in the middle is given by: | ||
| + | |||
| + | < | ||
| + | set-variable-value middle_cell_index max_cell_index/ | ||
| + | </ | ||
| + | |||
| + | Then select the cell in the centre by its c[e]ll index: | ||
| + | |||
| + | < | ||
| + | add-cell-selection E middle_cell_index $ select the cell in the centre. | ||
| + | </ | ||
| + | |||
| + | Assign the higher carbon content, '' | ||
| + | |||
| + | < | ||
| + | set-cell-property * V C C c_comp_layer | ||
| + | </ | ||
| + | |||
| + | The g5 window shows, by default, the selected cell(s) in blue: | ||
| + | |||
| + | {{: | ||
| + | |||
| + | |||
| + | |||
| + | |||
| + | ==== e. Display ==== | ||
| + | |||
| + | The next stage is to set up the graphical display to look at the results. Firstly, assign a property to the existing paint window created above, '' | ||
| + | |||
| + | < | ||
| + | set-gui-window-property paint_window_id I W V $ display cell variable | ||
| + | set-gui-window-property paint_window_id I V _cwp$c $ set which variable | ||
| + | set-gui-window-property paint_window_id I R 0..1 | ||
| + | </ | ||
| + | |||
| + | In the first line in the code above, the syntax represents s[i]mulation results, display [w]hat, [v]ariable. In the second, the [v]ariable is specified as '' | ||
| + | |||
| + | The cell array is shown schematically below, where blue represents a lower carbon content, and yellow a higher. | ||
| + | |||
| + | {{: | ||
| + | |||
| + | In addition, make a profile plot of the carbon content. Create a window of type **g1: plot grid: 1D profile**. | ||
| + | |||
| + | < | ||
| + | new-gui-window g1 | ||
| + | </ | ||
| + | |||
| + | Assign a name to both the plot window and the plot itself: | ||
| + | |||
| + | < | ||
| + | set-variable-value profile_window_id active_frame_id | ||
| + | set-variable-value profile_plot_id last_plot_id | ||
| + | </ | ||
| + | |||
| + | Add a [s]eries, [n]ew of the [s]imulation variable '' | ||
| + | |||
| + | < | ||
| + | set-plot-option profile_plot_id S N S _cwp$c | ||
| + | </ | ||
| + | |||
| + | Specify that the default x-axis should be used, and give it a title: | ||
| + | |||
| + | < | ||
| + | set-gui-window-property profile_window_id S U Y | ||
| + | set-gui-window-property profile_window_id S T Position [m] | ||
| + | </ | ||
| + | |||
| + | Label the y-axis and specify its range: | ||
| + | |||
| + | < | ||
| + | set-plot-option profile_plot_id A Y 1 T Carbon content [wt.%] | ||
| + | set-plot-option profile_plot_id A Y 1 S 0..c_comp_layer_wp | ||
| + | </ | ||
| + | |||
| + | Give the plotted series the name ' | ||
| + | |||
| + | < | ||
| + | set-plot-option profile_plot_id S M 0 Numerical | ||
| + | </ | ||
| + | |||
| + | Specify the range of cells over which the profile is to be taken. In the syntax here, the first line gives the index of the s[t]arting cell (here 0) and the second line gives the index of the st[o]p cell (here max_cell_index). | ||
| + | |||
| + | < | ||
| + | set-gui-window-property profile_window_id I t 0 | ||
| + | set-gui-window-property profile_window_id I o max_cell_index | ||
| + | </ | ||
| + | |||
| + | Give the window an appropriate size: | ||
| + | |||
| + | < | ||
| + | move-gui-window profile_window_id 20 20 600 400 | ||
| + | </ | ||
| + | |||
| + | The plot should look as shown: | ||
| + | |||
| + | {{: | ||
| + | |||
| + | We will now create a third window to display the time-dependent evolution of the carbon concentration in selected cells. The required window type is **g2: Plot grid: cell history**. Create the window and assign a name to the window and plot. | ||
| + | |||
| + | < | ||
| + | new-gui-window g2 | ||
| + | set-variable-value history_window_id active_frame_id | ||
| + | set-variable-value history_plot_id last_plot_id | ||
| + | </ | ||
| + | |||
| + | Use the following commands to specify that the default x-axis should be used and to assign it a title and a range. | ||
| + | |||
| + | < | ||
| + | set-gui-window-property history_window_id S U Y $use default x-axis | ||
| + | set-gui-window-property history_window_id S T Time [s] | ||
| + | set-gui-window-property history_window_id S S 0..sim_time | ||
| + | </ | ||
| + | |||
| + | Label the y-axis and specify its range: | ||
| + | |||
| + | < | ||
| + | set-plot-option history_plot_id A Y 1 T Carbon content [wt.%] | ||
| + | set-plot-option history_plot_id A Y 1 S 0..c_comp_layer_wp | ||
| + | </ | ||
| + | |||
| + | Now add two series, firstly the composition of the middle cell, and then that of the cell adjacent to it. :!: | ||
| + | |||
| + | < | ||
| + | set-plot-option history_plot_id S N S _cwp$c{middle_cell_index} | ||
| + | set-plot-option history_plot_id S N S _cwp$c{middle_cell_index-1} | ||
| + | </ | ||
| + | |||
| + | Give these series titles. Note that the quotation marks are needed for a multi-word title, as in the case of the series 'next to centre' | ||
| + | |||
| + | < | ||
| + | set-plot-option history_plot_id S M 0 centre | ||
| + | set-plot-option history_plot_id S M 1 'next to centre' | ||
| + | </ | ||
| + | |||
| + | Finally, resize and move the history plot. | ||
| + | |||
| + | < | ||
| + | move-gui-window history_window_id 800 20 600 400 | ||
| + | </ | ||
| + | |||
| + | |||
| + | ==== f. Analytical solution ==== | ||
| + | |||
| + | In the final stage before setting up and running the simulation, the analytical solution will be added to the existing profile plot. Since the expression is rather complicated, | ||
| + | |||
| + | '' | ||
| + | |||
| + | where '' | ||
| + | |||
| + | \[B = \int_{-\infty}^{\infty} C \{x,t\} \mathrm{d}x, | ||
| + | |||
| + | '' | ||
| + | |||
| + | \[ {1 \over 2 \sqrt{\pi D t}} \] | ||
| + | |||
| + | and '' | ||
| + | |||
| + | \[ \mathrm{exp}\left({-x^2 | ||
| + | |||
| + | The value of the term $B$ can be estimated by assuming that the initial concentration of carbon in the bar is small enough to be neglected, such that initially all the carbon is in the layer in the middle of the bar. The value of $B$ is then simply the carbon concentration in the layer (in the units in which the curve is to be plotted, here weight percent,'' | ||
| + | |||
| + | < | ||
| + | set-variable-value Bval c_comp_layer_wp*(barlength/ | ||
| + | </ | ||
| + | |||
| + | The cell variables (prefixed with **'' | ||
| + | |||
| + | < | ||
| + | set-function-expression term_1 1/ | ||
| + | </ | ||
| + | |||
| + | It should be noted that the term //x// in the mathematical formulation of the analytical expression is measured from the centre of the range, whereas the cell variable corresponding to the x-position of the centre of the cells, '' | ||
| + | |||
| + | < | ||
| + | set-variable-value midpt barlength/2 | ||
| + | set-function-expression x_shifted _ccenter_x-midpt | ||
| + | set-function-expression term_2 exp((-1*x_shifted**2)/ | ||
| + | </ | ||
| + | |||
| + | The analytical expression is then given by: | ||
| + | |||
| + | < | ||
| + | set-function-expression analytical_expr Bval*term_1*term_2 | ||
| + | </ | ||
| + | |||
| + | and this series is plotted onto the existing profile plot using the command: | ||
| + | |||
| + | < | ||
| + | set-plot-option profile_plot_id S N F analytical_expr 0..barlength | ||
| + | </ | ||
| + | |||
| + | The '' | ||
| + | |||
| + | Note that when the analytical expression is first plotted, it looks slightly strange. This is because the time is initially zero, leading to a value of zero in the denominator of the first term, '' | ||
| + | |||
| + | ==== g. Setting up and running the simulation ==== | ||
| + | |||
| + | The final set of commands deals with the simulation time and temperature, | ||
| + | |||
| + | < | ||
| + | set-simulation-parameter e sim_time $ end | ||
| + | </ | ||
| + | |||
| + | Specify the update interval: | ||
| + | |||
| + | < | ||
| + | set-simulation-parameter u 10 $ update interval | ||
| + | </ | ||
| + | |||
| + | It is now necessary to tell the simulation module what type of problem should be simulated, i.e. how temperature and diffusion should be treated. In this case, we wish to simulate a **diffusion-field problem**, in which the elements are allowed to diffuse through the system: | ||
| + | |||
| + | < | ||
| + | set-simulation-parameter d f $ diffusion-field simulation | ||
| + | </ | ||
| + | |||
| + | The temperature control, by contrast, is set to be isothermal and the hold temperature specified in terms of the variable '' | ||
| + | |||
| + | < | ||
| + | set-simulation-parameter t i sim_temp $ temperature: | ||
| + | </ | ||
| + | |||
| + | The final step is to set up the default boundary conditions, i.e. the boundary conditions that apply in the absence of any other boundary conditions being set for specific cells or groups of cells. The type of boundary condition required is [g]eometric, | ||
| + | |||
| + | < | ||
| + | set-default-boundary-cond G O $ open to environment. | ||
| + | </ | ||
| + | |||
| + | Save the script. | ||
| + | |||
| + | To run the script, click on **Script > Run Script**, press **Shift + F2** or click on the {{: | ||
| + | |||
| + | When the script has been correctly executed, select **Simulation > Start simulation** or press **Ctrl + J** to start the simulation. | ||
| + | |||
| + | < | ||
| + | ##################################################################### | ||
| + | |||
| + | | ||
| + | |||
| + | ##################################################################### | ||
| + | </ | ||
| + | |||
| + | and finishing with: | ||
| + | |||
| + | < | ||
| + | #################################################################################################### | ||
| + | |||
| + | |||
| + | microstructure simulation finished in 59.07 secs (tot 59.07 secs) on Fri, 2016-06-10, 17: | ||
| + | |||
| + | |||
| + | #################################################################################################### | ||
| + | </ | ||
| + | |||
| + | |||
| + | ===== Interpreting the results ===== | ||
| + | |||
| + | You will notice that, at first, the peak value of the analytical expression is much higher than the maximum value of 0.8 wt. % in the numerical simulation. The reason for this is that the analytical solution is based on the assumption that all the carbon is initially in an infinitesimally thin layer, whereas in the numerical solution, we have put the carbon into a cell of finite thickness. As the simulation proceeds, the two curves start to approach each other until they eventually almost exactly coincide. By default, the profile plot is set to show a ' | ||
| + | |||
| + | When the simulation is finished, you can examine the profiles at different times by looking at **Global > Buffers > Edit Buffer States**. You can obtain a set of profiles as shown below by loading up selected buffer states and using ‘Duplicate and Lock all Series’ on the plot right-click menu. | ||
| + | |||
| + | {{: | ||
| + | |||
| + | On the cell history plot, you should be able to see how the carbon concentration of the central cell decreases, while that of the adjacent cell decreases, until at longer exposure times, they converge. | ||
| + | |||
| + | |||
| + | {{: | ||
| + | |||
| + | |||
| + | ===== Subsequent articles ===== | ||
| + | |||
| + | Go to the [[: | ||