PDF

Pyrolysis of Wood with Time-Integrated Objective
Introduction
The process of pyrolyzing wood to produce tar and charcoal has been important since ancient times. Tar was used to impregnate wood for ships, and char was essential for iron smelting. Material for both of these applications have later been replaced by fossil sources, but environmental concerns has since kindled the interest for products produced by pyrolyzing wood.
Pyrolysis is thermal decomposition in the absence of an oxidizing atmosphere. In other words; heat something up without air until it decomposes. Historically, tar and charcoal were produced in piles, or pits, where wood was covered with for example dirt to prevent air from reaching the inside of the pile. The wood was lit on fire and allowed to smolder, but not burn (combust). During smoldering (pyrolysis), volatile species, water, and light decomposition products, leave the solid, resulting in charcoal. In modern times, steel reactors with an inert atmosphere are used.
The products that result from pyrolysis will depend on the feedstock type and particle size, the heating rate, the final temperature, and the duration of the process. Due to the complexity of the reaction mechanism, so called lumped-reaction models are often used. The reaction products are lumped into pseudospecies based on their phase, which gives a simplified reaction scheme that can be used for engineering purposes. One such reaction scheme is seen in Figure 1.
Figure 1: The reaction scheme used in this model consists of four pseudospecies, namely gas (g), tar (t), intermediate solid (is), and char (c).
The scheme was proposed by Park and others (Ref. 1) to describe the pyrolysis of a wooden sphere, approximately 1 cm in diameter, inserted into a hot furnace. The experimental setup is illustrated in Figure 2.
Figure 2: The experimental system setup consists of an isothermal furnace with inert atmosphere, a sample holder attached to a scale, and thermocouples (TC) measuring the temperature of the system. The sample radius is approximately 1 cm.
The experimental system consists of an isothermal furnace with inert atmosphere. The temperature of the furnace, measured by thermocouples (TC), is kept constant, and the inert atmosphere is achieved by nitrogen flowing through the furnace chamber. For each experiment, the sample is inserted into the isothermal furnace, and the sample temperature and the sample mass are recorded during the pyrolysis process. The temperature gradient within the sample is significant and the temperature is thus measured at three positions within the sample; at the surface, mid and center position. Wood is a porous, anisotropic material and in this study the temperature was measured along the fibers in the horizontal direction.
The reaction scheme in Figure 1 describes both primary and secondary pyrolysis reactions. The primary decomposition steps will convert wood into the pseudospecies gas, tar, and intermediate solid. Gaseous species are those that do not condense at room temperature, for example carbon monoxide. Tar species are all the condensable volatiles, for example water, carboxylic acids and phenols. The gases and the tars leaving the particle result in mass loss. On its way out from the porous particle, the tar may decompose to form gas or char. The intermediate solid further converts into char.
Experimental results from Ref. 1 showing development of the sample mass and the temperature are shown in Figure 3.
Figure 3: Experimental data from Ref. 1. The mass has been normalized by the initial sample mass.
This example model, based on Ref. 1 and Ref. 2, consists of two parts. The first part demonstrates how to set up a model describing the pyrolysis process, heat, and momentum transfer in an anisotropic wood sphere. In the second part, parameter estimation is used to optimize the model using the experimental data in Figure 3. The parameters to be estimated are one Arrhenius constant, two reaction heats, and one external heat transfer coefficient (see Table 1).
Model Definition
The pyrolysis of a centimeter-sized wood particle presents a fully coupled multiphysics problem with mass transfer, fluid flow and heat transfer. In this example, both the conductive heat transfer and the permeability of the solid are anisotropic.
Mass Transport
The pyrolysis reaction scheme in Figure 1 consists of gaseous and solid species. The reaction rate expressions (kg/(m3·s)), for the solid species wood (w), intermediate solid (is), and char (c) are expressed in terms of their respective density in the manner of
(1)
(2)
and
(3)
The density is defined as the bulk density of the solid, and thus includes the porosity ε of the solid domain. ki is the Arrhenius rate constant (1/s), as indicated in the reaction scheme (Figure 1)
(4)
No transport terms are needed for the solid species, and Equation 1Equation 3 are thus sufficient to conserve the mass of the solid species.
The mass conservation equation for the gas species i includes diffusion, convection, and the reaction rate terms. The gas mixture inside the particle consists solely of gas, tar, and inert gas, the mass balance is expressed in terms of the respective mass fractions ω as
(5)
Here, ρ is the density of the fluid in the pores, derived using the ideal gas law, and ε is the porosity of the porous domain:
The initial wood porosity, εw0, is 0.4.
In Equation 5, ji is the diffusional flux as described by Fick’s law, with diffusion coefficients Di, and a Millington and Quirk model to derive the effective diffusivity:
The mass averaged velocity of the mixture in the pores, u, is derived using Darcy’s law:
(6)
Here, μ is the viscosity (kg/(m·s)), p is the pressure (Pa), and κ is the effective permeability (m2):
where j indicates across or along the fiber direction.
The reaction rate expressions for the gas phase species tar (t) and gas (g) are
and
Mass transfer through the exterior boundary of the particle is dominated by convection. A boundary condition assuming no diffusive flux is thus applicable:
Here, n denotes the outward pointing normal of the exterior boundary.
Momentum Transport
The fluid flow is defined with the continuity equations
with a Darcian flow (see Equation 6), and a mass source term (evolution of fluid species) defined as
(7)
At the exterior boundary, a zero relative pressure with respect to the reference pressure (pref = 1 atm) is prescribed:
Heat Transport
The energy balance equation applied to the fluid (f) in the pores, and the solid bulk (b) of the wood sample assumes local thermal equilibrium, and considers heat transfer through convection, radiation, and conduction:
(8)
Here, Q is the heat of reaction
(9)
In Equation 8 (ρCp)eff is defined as
The heat capacities at constant pressures for the fluid and bulk phases are
and
The dry bulk density ρb is (1 − ε)ρw,0 and the fluid density ρf is derived from the ideal gas law.
The effective thermal conductivity is the weighted sum of the conductivity of the fluid, kf, of the solid bulk, kb, as well as a contribution from the radiation in the pores:
where σ is the Boltzmann constant (W/(m2·K4)), e is the emissivity, and deff is the effective pore diameter (m) defined as the weighted sum of the pore diameter in the wood and the char:
The degree of pyrolysis, η is
The conductivity of the solid bulk is anisotropic
where j indicates across or along the fiber direction.
The external boundary of the particle has a heat flux boundary condition:
where the heat flux q0 is the sum of convective and radiative heat flux:
(10)
Here, hconv is the heat transfer coefficient in the gas surrounding the particle, Tgas is the gas temperature in the reactor, σ is the Stefan–Boltzmann constant, es is the surface emissivity, Treactor is the reactor temperature, and T is the surface temperature of the sample. All temperatures are expressed in K.
Parameter Estimation
Parameter estimation problems consist of three components: (i) experimental data; (ii) a forward model that represents the physics of the experiments; and (iii) an optimization algorithm that compares the two and updates the model parameters to minimize the difference. This can be formulated mathematically as a nonlinear least-squares minimization problem, but in this model we will use an Interpolation function to convert the experimental data to a time varying function, so that the objective can be expressed as
(11)
with
(12)
Herein, q is the vector of control parameters (ξ) that we want to estimate, N is the number of experiments, is the interpolation function of experiment n, and Pn(u(q), q, t) denotes the corresponding model prediction given the PDE solution u.
In this example, we consider N = 4 experimental data sets from Ref. 1 (see Figure 3), for which the measured quantity Pn is either the temperature (at one of the thermocouple positions, see Figure 2) or the normalized solid mass, and u is the solution to the multiphysics model set up to describe the system.
The normalized solid mass Y is defined as
(13)
In this model, the control parameters to estimate are q = (AisΔHtΔHchconv), where Ais is the Arrhenius frequency factor for the primary pyrolysis step where wood turns to intermediate solid (Equation 4). ΔHt and ΔHc are the heat of reaction for formation of tar from wood, and formation of char from intermediate solid (Equation 9). In this model the primary reactions all share the same heat of reactions, namely ΔHt = ΔHg = ΔHis. The final control parameter hconv is the convective heat transfer coefficient external to the particle (Equation 10). The parameters along with the initial guess of their values are provided in Table 1.
Results and Discussion
The results from the forward model, using the initial guess of the parameter values in Table 1, are shown in Figure 4 and Figure 5. The model describes the trends in the temperatures and solid mass quite well, especially for the middle temperature (Figure 5). Both the timings and the absolute values of the peak temperatures for each position are lower than in experiments, and the experimental final solid mass is not captured by the model at all.
The model predictions after parameter estimation are illustrated in Figure 4 and Figure 5. For comparison, the results from the forward model are also included in those figures. It is clear that the optimized model better predicts the experimental results. The timing of the center temperature peak and its value are now captured, as is the final solid mass.
Figure 4: Surface and middle temperatures from the forward and the optimized model.
Figure 5: Center temperature and normalized solid mass from the forward and optimized model.
Figure 6 illustrates the changes in solid composition in the particle at three different times. Early in the process, there is mainly wood. This wood is converted into gases and the solid intermediate species (is). Late in the process, the wood is fully converted, and most of the particle consists of char.
Figure 6: Normalized densities of wood ρw/ρw,0, intermediate solid ρis/ρw,0, and char ρc/ρw,0, at three different times.
Figure 4 and Figure 5 above show the temperatures at the thermocouple positions in the particle. Figure 7Figure 9 illustrate the temperature in the solid domain, together with both the mass source (Qm in Equation 7) and the heat source (Q in Equation 9). As expected, at an early stage (Figure 7), a positive mass source is accompanied with a negative heat source, since the primary pyrolysis reactions, producing gaseous species, are endothermic.
Figure 7: Temperature, mass source, and heat source in the modeled geometry, 150 s into the pyrolysis process.
As the process progresses (Figure 8), the positive mass source moves inward, followed by a positive heat source, indicating the formation of char. This is also seen in Figure 6.
Figure 8: Temperature, mass source, and heat source in the modeled geometry, 270 s into the pyrolysis process.
At the last stage of the mass loss process, depicted in Figure 9, there are practically no gas forming reactions, and only exothermic processes, as seen by a dominating negative mass source and a positive heat source.
Figure 9: Temperature, mass source, and heat source in the modeled geometry, 433 s into the pyrolysis process. This time corresponds to the temperature peak seen in Figure 5.
One thing seen in Figure 6-Figure 9 is the anisotropic particle properties. Figure 10 illustrates this further by showing the relative pressure in the particle, p/pref, the total Darcy velocity magnitude, U, the Darcy velocity vector, u, the porosity, ε, and the solid mass fraction, Y (see Equation 13).
Figure 10: The relative pressure in the particle, p/pref, the total Darcy velocity magnitude, U, the Darcy velocity, u, the porosity, ε, and the solid mass fraction, Y.
Notes About the COMSOL Implementation
The solid phase reactions in Figure 1 are implemented in the Domain ODEs and DAEs interface.
In parameter estimation problems, it is good practice to first set up and test the forward model before solving the inverse problem. In this example, the normalized solid mass, Y, is tracked during computation using a Domain Probe. The data associated with this probe is then used in postprocessing to illustrate the data. For the optimization study, create a new domain probe and disable the first one, as otherwise the data from the forward study will be overwritten with the data from the optimization run.
This models does not use the Parameter Estimation functionality available in COMSOL Multiphysics. Instead, an objective based on integration is constructed. COMSOL evaluates transient objectives at the final time, so this requires the use of an ODE. Interpolation functions are used to estimate the transient data between measurements.
For most least-squares problems, the Levenberg–Marquardt algorithm with a finite difference approximation of the Jacobian is a robust and efficient choice of optimization solver, but it requires the evaluation of the gradient for every measurement. If this is too time consuming, a gradient-free algorithm like the BOBYQA can be used or gradient based method can be used together with an aggregated objective — as demonstrated here. To increase the stability of the optimization process, the logarithm of the control parameters can be optimized. This gives scales of 1, and initial values equal to zero. This strategy is used in this example model together with the gradient based IPOPT optimization algorithm.
References
1. W.C. Park, A. Atreya, and H.R. Baum, “Experimental and theoretical investigation of heat and mass transfer processes during wood pyrolysis,” Combustion and Flame, vol. 157, pp. 481–494, 2010.
2. X. Shi, F. Ronsse, and J.G. Pieters, “Finite element modeling of intraparticle heterogeneous tar conversion during pyrolysis of woody biomass particles,” Fuel Processing Technology, vol. 148, pp. 302-316, 2016.
Application Library path: Optimization_Module/Parameter_Estimation/pyrolysis_wood_odeobj
Modeling Instructions
This model consists of two parts: setting up the forward model and performing parameter estimation to calibrate that model with experimental data. Start by setting up the forward model.
Use the Model Wizard to add a Component (2D axisymmetric due to the anisotropic sphere); the physics Darcy’s Law, Transport of Concentrated Species in Porous Media, Heat Transfer in Porous Media, and Domain ODEs and DAEs (for the solid reactions); and a Time Dependent study to follow the decomposition progress.
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Chemical Species Transport > Transport of Concentrated Species in Porous Media (tcs).
3
Click Add.
4
In the Select Physics tree, select Fluid Flow > Porous Media and Subsurface Flow > Darcy’s Law (dl).
5
Click Add.
6
In the Select Physics tree, select Heat Transfer > Porous Media > Heat Transfer in Porous Media (ht).
7
Click Add.
8
In the Select Physics tree, select Mathematics > ODE and DAE Interfaces > Domain ODEs and DAEs (dode).
9
Click Add.
10
In the Select Physics tree, select Mathematics > ODE and DAE Interfaces > Global ODEs and DAEs (ge).
11
Click Add.
12
Click  Study.
13
In the Select Study tree, select General Studies > Time Dependent.
14
Before setting up the geometry, load all the parameters and variables for this model from files. Since we have not yet defined the physics, some of the expressions in the variable files will be undefined at this point.
Global Definitions
Sample Properties
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Sample Properties in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file pyrolysis_wood_odeobj_sample_properties_parameters.txt.
Experimental Conditions
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Experimental Conditions in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file pyrolysis_wood_odeobj_experimental_conditions_parameters.txt.
Reaction Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Reaction Parameters in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Optimization Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Optimization Parameters in the Label text field.
3
Locate the Parameters section. Click  Load from File.
4
Definitions
Solid Species Variables
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, type Solid Species Variables in the Label text field.
3
Locate the Variables section. Click  Load from File.
4
Reaction Variables
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, type Reaction Variables in the Label text field.
3
Locate the Variables section. Click  Load from File.
4
Fluid Species Variables
1
Right-click Definitions and choose Variables.
2
In the Settings window for Variables, type Fluid Species Variables in the Label text field.
3
Locate the Variables section. Click  Load from File.
4
External Boundary Variables
1
Right-click Definitions and choose Variables.
2
In the Settings window for Variables, type External Boundary Variables in the Label text field.
3
Locate the Variables section. Click the Load button. From the menu, choose Load from File.
4
Set up the geometry. It consists of a quarter of a Circle, with a Circular Arc used only for meshing, and a Point at the location of the middle thermocouple.
Geometry 1
Circle 1 (c1)
1
In the Geometry toolbar, click  Circle.
2
In the Settings window for Circle, locate the Size and Shape section.
3
In the Radius text field, type r_sample.
4
In the Sector angle text field, type 90.
Circular Arc 1 (ca1)
1
In the Geometry toolbar, click  More Primitives and choose Circular Arc.
2
In the Settings window for Circular Arc, locate the Radius section.
3
In the Radius text field, type r_sample/3.
4
Click  Build All Objects.
5
Click the  Zoom Extents button in the Graphics toolbar.
Middle Along
1
In the Geometry toolbar, click  Point.
2
In the Settings window for Point, type Middle Along in the Label text field.
3
Locate the Point section. In the r text field, type r_sample/2.
4
Click  Build All Objects.
Mesh Control Edges 1 (mce1)
1
In the Geometry toolbar, click  Virtual Operations and choose Mesh Control Edges.
2
On the object fin, select Boundary 7 only.
3
In the Geometry toolbar, click  Build All.
Domain ODEs and DAEs (dode)
In the Domain ODEs and DAEs interface, the solid species reaction rates are added. Using this interface allows choosing the unit for the dependent variables, namely the densities.
1
In the Model Builder window, under Component 1 (comp1) click Domain ODEs and DAEs (dode).
2
In the Settings window for Domain ODEs and DAEs, locate the Units section.
3
Click  Select Dependent Variable Quantity.
4
In the Physical Quantity dialog, type density in the text field.
5
In the tree, select General > Density (kg/m^3).
6
7
In the Settings window for Domain ODEs and DAEs, locate the Units section.
8
In the Source term quantity table, enter the following settings:
9
Click to expand the Dependent Variables section. In the Field name (kg/m³) text field, type rho.
10
In the Number of dependent variables text field, type 3.
11
In the Dependent variables (kg/m³) table, enter the following settings:
Distributed ODE 1
In the Source Term fields, add the decomposition rate expressions for the solid species.
1
In the Model Builder window, under Component 1 (comp1) > Domain ODEs and DAEs (dode) click Distributed ODE 1.
2
In the Settings window for Distributed ODE, locate the Source Term section.
3
In the f text-field array, type -(k_t+k_g+k_is)*rho_w on the first row.
4
In the f text-field array, type k_is*rho_w-k_c*rho_is on the second row.
5
In the f text-field array, type k_c*rho_is+k_c2*tcs.rho*w_t on the third row. The dependent variable w_t is not yet defined. It will be added in the Transport of Concentrated Species in Porous Media Interface.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the rhow text field, type rho_w_init.
Transport of Concentrated Species in Porous Media (tcs)
Continue by setting up the mass transfer equations and reactions. The workflow is to first go through the settings for the default features, and then add the features that are needed.
1
In the Model Builder window, under Component 1 (comp1) click Transport of Concentrated Species in Porous Media (tcs).
2
In the Settings window for Transport of Concentrated Species in Porous Media, locate the Transport Mechanisms section.
3
From the Diffusion model list, choose Fick’s law.
4
Click to expand the Dependent Variables section. In the Number of species text field, type 3.
5
In the Mass fractions (1) table, enter the following settings:
6
Locate the Species section. From the From mass constraint list, choose w_N2.
Species Molar Masses 1
1
In the Model Builder window, under Component 1 (comp1) > Transport of Concentrated Species in Porous Media (tcs) click Species Molar Masses 1.
2
In the Settings window for Species Molar Masses, locate the Molar Mass section.
3
In the Mwt text field, type Mw_t.
4
In the Mwg text field, type Mw_g.
5
In the MwN2 text field, type Mw_N2. The parameters for the molar masses were loaded from file and can be found in the Sample Properties node under Global Definitions.
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Transport of Concentrated Species in Porous Media (tcs) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Convection section.
3
From the u list, choose Total Darcy velocity field (dl/porous1).
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
In the εp text field, type epsilon. This parameter was also loaded from file.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Transport of Concentrated Species in Porous Media (tcs) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the ω0,wt text field, type 0.
4
In the ω0,wg text field, type 0.
Since the geometry we want to model is actually a sphere, add a Symmetry boundary condition.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Finally, define the reactions. Add two Reaction Sources features: one with the reactions that involve mass transfer to other phases, and one with only gas phase species.
Reaction Sources with Phase Transfer
1
In the Physics toolbar, click  Domains and choose Reaction Sources.
2
In the Settings window for Reaction Sources, type Reaction Sources with Phase Transfer in the Label text field.
3
4
Locate the Reactions section. Select the Mass transfer to other phases checkbox.
5
In the Rwt text field, type k_t* rho_w-k_c2*w_t*tcs.rho.
6
In the Rwg text field, type k_g*rho_w.
Reaction Sources Gas to Gas
1
In the Physics toolbar, click  Domains and choose Reaction Sources.
2
In the Settings window for Reaction Sources, type Reaction Sources Gas to Gas in the Label text field.
3
4
Locate the Reactions section. In the Rwt text field, type -k_g2*w_t*tcs.rho.
5
In the Rwg text field, type k_g2*w_t*tcs.rho.
Darcy’s Law (dl)
Now, define the fluid flow in the system.
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the ρ list, choose Density (tcs/porous1/fluid1).
4
From the μ list, choose User defined. In the associated text field, type viscosity.
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type epsilon.
4
From the κ list, choose User defined. From the list, choose Diagonal.
5
Specify the κ matrix as
Initial Values 1
Now we have gone through the settings for the default features. Next we need to add a Mass Source feature, a Pressure boundary condition feature for the external surface, and a Symmetry feature.
Mass Source 1
1
In the Physics toolbar, click  Domains and choose Mass Source.
2
3
In the Settings window for Mass Source, locate the Mass Source section.
4
Click in the Qm text field, then press Ctrl+Space. From the menu, choose Component 1 (comp1) > Transport of Concentrated Species in Porous Media > tcs.Qmass - Net mass source - kg/(m³·s).
Pressure 1
1
In the Physics toolbar, click  Boundaries and choose Pressure.
2
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Heat Transfer in Porous Media (ht)
The last physics interface to set up is the heat transfer interface. Use the same workflow: go through the default features and then add the additional features that you need. In this case, add the features Symmetry, Heat Flux, and Heat Source.
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Heat Convection section.
3
From the u list, choose Total Darcy velocity field (dl/porous1).
4
Locate the Model Input section. From the pA list, choose Absolute pressure (dl).
5
Locate the Heat Conduction, Fluid section. From the kf list, choose User defined. In the associated text field, type k_f.
6
Locate the Thermodynamics, Fluid section. From the ρf list, choose Density (tcs/porous1/fluid1).
7
From the Cp,f list, choose User defined. In the associated text field, type cp_f.
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type epsilon.
4
Locate the Heat Conduction, Porous Matrix section. From the kb list, choose User defined. From the list, choose Diagonal.
5
Specify the kb matrix as
6
Locate the Thermodynamics, Porous Matrix section. From the ρb list, choose User defined. In the associated text field, type rho_b.
7
From the Cp,b list, choose User defined. In the associated text field, type cp_b.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type T0.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Heat Flux 1
1
In the Physics toolbar, click  Boundaries and choose Heat Flux.
2
3
In the Settings window for Heat Flux, locate the Heat Flux section.
4
In the q0 text field, type q0.
Heat Source 1
1
In the Physics toolbar, click  Domains and choose Heat Source.
2
3
In the Settings window for Heat Source, locate the Heat Source section.
4
In the Q0 text field, type Q.
Before we can compute the study, we need to define the mesh that discretizes the modeling domain into finite elements.
Mesh 1
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Sequence Type section.
3
From the list, choose User-controlled mesh.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Finer.
Mapped 1
1
In the Mesh toolbar, click  Mapped.
2
3
In the Settings window for Mapped, locate the Domain Selection section.
4
From the Geometric entity level list, choose Domain.
5
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type 25.
Distribution 2
1
In the Model Builder window, right-click Mapped 1 and choose Distribution.
2
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 15.
Size 1
1
In the Model Builder window, right-click Free Triangular 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Point.
4
5
Locate the Element Size section. From the Predefined list, choose Finer.
6
In the Model Builder window, right-click Mesh 1 and choose Build All.
Study 1: Forward Model (Initial-Value Based)
Set up the study for the forward model. Add the experimental data to compare it to the results from the forward model. We also add a probe to derive Y (the normalized solid mass) at the same time as we compute the study.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Study 1: Forward Model (Initial-Value Based) in the Label text field.
3
Locate the Study Settings section. Clear the Generate default plots checkbox.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.01*tmax,tmax).
4
Click to expand the Results While Solving section. From the Probes list, choose None.
5
Locate the Physics and Variables Selection section. In the Solve for column of the table, under Component 1 (comp1), clear the checkbox for Global ODEs and DAEs (ge).
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node. The choice Show Default Solver creates the node Solver Configurations where we can edit the solver settings. Since we know the scales for the dependent variables, we will enter them and not use the default values. If a scale is too high (orders higher than the value of the dependent variable), then we will not get an accurate solution for that variable. If instead the scales are too low, the solver will take more time steps than necessary, giving high accuracy but increasing the computation time.
3
In the Model Builder window, expand the Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node, then click Pressure (comp1.p).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Dependent Variable Rho_c (comp1.rho_c).
7
In the Settings window for Field, locate the Scaling section.
8
From the Method list, choose Manual.
9
In the Scale text field, type rho_w_init.
10
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Dependent Variable Rho_is (comp1.rho_is).
11
In the Settings window for Field, locate the Scaling section.
12
From the Method list, choose Manual.
13
In the Scale text field, type rho_w_init.
14
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Dependent Variable Rho_w (comp1.rho_w).
15
In the Settings window for Field, locate the Scaling section.
16
From the Method list, choose Initial-value based.
17
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Temperature (comp1.T).
18
In the Settings window for Field, locate the Scaling section.
19
From the Method list, choose Initial-value based.
20
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Mass Fraction (comp1.w_g).
21
In the Settings window for Field, locate the Scaling section.
22
From the Method list, choose Manual.
23
In the Scale text field, type 0.1.
24
In the Model Builder window, under Study 1: Forward Model (Initial-Value Based) > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 click Mass Fraction (comp1.w_t).
25
In the Settings window for Field, locate the Scaling section.
26
From the Method list, choose Manual.
27
In the Scale text field, type 0.1.
Since we want to compare the forward model with the experimental data, add the experimental data. We should also add a probe to derive Y (the normalized solid mass) while solving the model.
Results
Experimental Data
1
In the Model Builder window, expand the Results node.
2
Right-click Results > Tables and choose Node Group.
3
In the Settings window for Group, type Experimental Data in the Label text field.
Experimental data: Y
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, locate the Data section.
3
Click  Import.
4
5
In the Label text field, type Experimental data: Y.
6
Locate the Column Headers section. In the table, enter the following settings:
Experimental data: T_surface
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Experimental data: T_surface in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file pyrolysis_wood_odeobj_experimental_data_T_surface.txt.
5
Locate the Column Headers section. In the table, enter the following settings:
Experimental data: T_middle
1
Right-click Experimental data: T_surface and choose Duplicate.
2
In the Settings window for Table, type Experimental data: T_middle in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file pyrolysis_wood_odeobj_experimental_data_T_middle.txt.
5
Locate the Column Headers section. In the table, enter the following settings:
Experimental data: T_center
1
Right-click Experimental data: T_middle and choose Duplicate.
2
In the Settings window for Table, type Experimental data: T_center in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file pyrolysis_wood_odeobj_experimental_data_T_center.txt.
5
Locate the Column Headers section. In the table, enter the following settings:
Plot the experimental data. Follow these instructions to generate Figure 3 in the model documentation.
Experimental Data
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Experimental Data in the Label text field.
3
Locate the Data section. From the Dataset list, choose None.
4
Click to expand the Title section. From the Title type list, choose None.
Normalized Solid Mass
1
Right-click Experimental Data and choose Table Graph.
2
In the Settings window for Table Graph, locate the Coloring and Style section.
3
Find the Line style subsection. From the Line list, choose None.
4
From the Color list, choose Red.
5
Find the Line markers subsection. From the Marker list, choose Plus sign.
6
Click to expand the Legends section. Select the Show legends checkbox.
7
Find the Include subsection. Select the Label checkbox.
8
Clear the Headers checkbox.
9
In the Label text field, type Normalized Solid Mass.
Surface Temperature
1
Right-click Normalized Solid Mass and choose Duplicate.
2
In the Settings window for Table Graph, type Surface Temperature in the Label text field.
3
Locate the Data section. From the Table list, choose Experimental data: T_surface.
4
Locate the Coloring and Style section. From the Color list, choose Magenta.
5
Find the Line markers subsection. From the Marker list, choose Asterisk.
Middle Temperature
1
Right-click Surface Temperature and choose Duplicate.
2
In the Settings window for Table Graph, type Middle Temperature in the Label text field.
3
Locate the Data section. From the Table list, choose Experimental data: T_middle.
4
Locate the Coloring and Style section. From the Color list, choose Black.
5
Find the Line markers subsection. From the Marker list, choose Circle.
Center Temperature
1
Right-click Middle Temperature and choose Duplicate.
2
In the Settings window for Table Graph, type Center Temperature in the Label text field.
3
Locate the Data section. From the Table list, choose Experimental data: T_center.
4
Locate the Coloring and Style section. From the Color list, choose Blue.
5
Find the Line markers subsection. From the Marker list, choose Point.
Experimental Data
1
In the Model Builder window, click Experimental Data.
2
In the Settings window for 1D Plot Group, locate the Plot Settings section.
3
Select the Two y-axes checkbox.
4
Select the x-axis label checkbox. In the associated text field, type Time (s).
5
Select the y-axis label checkbox. In the associated text field, type Normalized Solid Mass (-).
6
Select the Secondary y-axis label checkbox. In the associated text field, type Temperature (K).
7
In the table, select the Plot on secondary y-axis checkboxes for Surface Temperature, Middle Temperature, and Center Temperature.
8
Locate the Legend section. From the Position list, choose Middle right.
9
In the Experimental Data toolbar, click  Plot.
10
Click the  Zoom Extents button in the Graphics toolbar.
We need to enter a model expression for the surface temperature. Add a probe feature for this purpose. Also, add the probe to derive Y during computations.
Definitions
Probes for Parameter Estimation
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Node Group.
2
In the Settings window for Group, type Probes for Parameter Estimation in the Label text field.
Point Probe Surface
1
In the Definitions toolbar, click  Probes and choose Point Probe.
2
In the Settings window for Point Probe, type Point Probe Surface in the Label text field.
3
In the Variable name text field, type T_surface.
4
5
Locate the Expression section. In the Expression text field, type T.
6
Select the Description checkbox.
Point Probe Middle
1
Right-click Point Probe Surface and choose Duplicate.
2
In the Settings window for Point Probe, type Point Probe Middle in the Label text field.
3
In the Variable name text field, type T_middle.
4
Point Probe Center
1
Right-click Point Probe Middle and choose Duplicate.
2
In the Settings window for Point Probe, type Point Probe Center in the Label text field.
3
In the Variable name text field, type T_center.
4
Domain Probe Y
1
In the Definitions toolbar, click  Probes and choose Domain Probe.
2
In the Settings window for Domain Probe, type Domain Probe Y in the Label text field.
3
In the Variable name text field, type domY.
4
Locate the Expression section. In the Expression text field, type Y.
5
Select the Description checkbox.
Study 1: Forward Model (Initial-Value Based)
Now we are ready to compute the forward study.
1
In the Study toolbar, click  Compute.
Before starting the optimization study, verify that the model has conservation of mass. In other words, check that the mass in the system at any time equals the initial mass in the system.
Results
Mass Conservation Check
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, locate the Transformation section.
3
From the Transformation type list, choose General.
4
Select the Keep child nodes checkbox.
5
In the Column header text field, type m(t)/m(0) = 1.
6
In the Label text field, type Mass Conservation Check.
Gas and Tar Inside Sample
1
Right-click Mass Conservation Check and choose Integration > Surface Integration.
2
In the Settings window for Surface Integration, type Gas and Tar Inside Sample in the Label text field.
3
4
Locate the Expressions section. Click  Clear Table.
5
Gas and Tar Leaving Sample
1
In the Model Builder window, right-click Mass Conservation Check and choose Integration > Line Integration.
2
In the Settings window for Line Integration, type Gas and Tar Leaving Sample in the Label text field.
3
Locate the Expressions section. Click  Clear Table.
4
5
6
Locate the Data Series Operation section. From the Transformation list, choose Integral.
7
Select the Cumulative checkbox.
Intermediate + Char
1
Right-click Mass Conservation Check and choose Integration > Surface Integration.
2
3
In the Settings window for Surface Integration, type Intermediate + Char in the Label text field.
4
Locate the Expressions section. Click  Clear Table.
5
Wood
1
Right-click Mass Conservation Check and choose Integration > Surface Integration.
2
In the Settings window for Surface Integration, type Wood in the Label text field.
3
Locate the Expressions section. Click  Clear Table.
4
5
Mass Conservation Check
1
In the Model Builder window, click Mass Conservation Check.
2
In the Settings window for Evaluation Group, locate the Transformation section.
3
In the Expression text field, type (int1+int2+int3+int4).
4
In the Mass Conservation Check toolbar, click  Evaluate.
Mass Conservation Check
1
Go to the Mass Conservation Check window.
Normalize with the mass of wood at time 0 s.
2
In the Settings window for Evaluation Group, locate the Transformation section.
3
In the Expression text field, type (int1+int2+int3+int4)/0.0026954.
4
In the Mass Conservation Check toolbar, click  Evaluate.
5
Go to the Mass Conservation Check window.
By looking at the values in column 2 in the table of the Mass Conservation Check, evaluation group, we see that the mass is conserved with a precision of at least three decimals.
Global Definitions
Interpolation 1 (int1)
1
In the Home toolbar, click  Functions and choose Global > Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
From the Data source list, choose Result table.
4
Locate the Data Column Settings section. In the table, click to select the cell at row number 1 and column number 3.
5
In the Unit text field, type s.
Interpolation 2 (int2)
1
Right-click Interpolation 1 (int1) and choose Duplicate.
2
In the Settings window for Interpolation, locate the Definition section.
3
From the Table from list, choose Experimental data: T_surface.
4
Locate the Data Column Settings section. In the Unit text field, type K.
Interpolation 3 (int3)
1
Right-click Interpolation 2 (int2) and choose Duplicate.
2
In the Settings window for Interpolation, locate the Definition section.
3
From the Table from list, choose Experimental data: T_middle.
Interpolation 4 (int4)
1
Right-click Interpolation 3 (int3) and choose Duplicate.
2
In the Settings window for Interpolation, locate the Definition section.
3
From the Table from list, choose Experimental data: T_center.
Definitions
Objective
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, type Objective in the Label text field.
3
Locate the Variables section. In the table, enter the following settings:
Global ODEs and DAEs (ge)
Global Equations 1 (ODE1)
1
In the Model Builder window, under Component 1 (comp1) > Global ODEs and DAEs (ge) click Global Equations 1 (ODE1).
2
In the Settings window for Global Equations, locate the Global Equations section.
3
Now add the Parameter Estimation study.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Time Dependent.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2: Parameter Estimation
1
In the Settings window for Study, type Study 2: Parameter Estimation in the Label text field.
2
Locate the Study Settings section. Clear the Generate default plots checkbox.
General Optimization
1
In the Study toolbar, click  Optimization and choose General Optimization.
2
In the Settings window for General Optimization, locate the Optimization Solver section.
3
From the Method list, choose IPOPT.
4
In the Optimality tolerance text field, type 0.01.
5
Click Add Expression in the upper-right corner of the Objective Function section. From the menu, choose Component 1 (comp1) > Global ODEs and DAEs > comp1.obj - State variable obj - 1.
6
Locate the Control Variables and Parameters section. Click  Add four times.
7
8
Click to expand the Solver Settings section. Find the Objective settings subsection. From the Objective scaling list, choose Initial solution based.
9
Click to expand the Output section. From the Probes list, choose None.
Step 1: Time Dependent
1
In the Model Builder window, click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,tmax/120,tmax).
Edit the scales for the dependent variables in the same way as for the Forward Study.
Solution 2 (sol2)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 2 (sol2) node.
3
In the Model Builder window, expand the Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 node, then click Pressure (comp1.p).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Dependent Variable Rho_c (comp1.rho_c).
7
In the Settings window for Field, locate the Scaling section.
8
From the Method list, choose Manual.
9
In the Scale text field, type rho_w_init.
10
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Dependent Variable Rho_is (comp1.rho_is).
11
In the Settings window for Field, locate the Scaling section.
12
From the Method list, choose Manual.
13
In the Scale text field, type rho_w_init.
14
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Dependent Variable Rho_w (comp1.rho_w).
15
In the Settings window for Field, locate the Scaling section.
16
From the Method list, choose Initial-value based.
17
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Temperature (comp1.T).
18
In the Settings window for Field, locate the Scaling section.
19
From the Method list, choose Initial-value based.
20
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Mass Fraction (comp1.w_g).
21
In the Settings window for Field, locate the Scaling section.
22
From the Method list, choose Manual.
23
In the Scale text field, type 0.1.
24
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Dependent Variables 1 click Mass Fraction (comp1.w_t).
25
In the Settings window for Field, locate the Scaling section.
26
From the Method list, choose Manual.
27
In the Scale text field, type 0.1.
28
In the Model Builder window, expand the Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Optimization Solver 1 node, then click Time-Dependent Solver 1.
29
In the Settings window for Time-Dependent Solver, locate the General section.
30
From the Times to store list, choose Steps taken by solver.
31
In the Model Builder window, expand the Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Optimization Solver 1 > Time-Dependent Solver 1 node, then click Advanced.
32
In the Settings window for Advanced, locate the General section.
33
From the Solver log list, choose Minimal.
34
In the Model Builder window, under Study 2: Parameter Estimation > Solver Configurations > Solution 2 (sol2) > Optimization Solver 1 > Time-Dependent Solver 1 click Direct, pressure (dl) (Merged).
35
In the Settings window for Direct, locate the General section.
36
From the Solver list, choose PARDISO.
Follow these instructions to generate Figure 4.
Results
Experimental Data 1
In the Model Builder window, right-click Experimental Data and choose Duplicate.
Experimental Data 2
Right-click Experimental Data and choose Duplicate.
Optimized, Forward Model, and Experimental Data: T_surface and T_middle
1
In the Settings window for 1D Plot Group, type Optimized, Forward Model, and Experimental Data: T_surface and T_middle in the Label text field.
2
In the Model Builder window, expand the Optimized, Forward Model, and Experimental Data: T_surface and T_middle node.
Center Temperature, Normalized Solid Mass
1
In the Model Builder window, under Results > Optimized, Forward Model, and Experimental Data: T_surface and T_middle, Ctrl-click to select Normalized Solid Mass and Center Temperature.
2
Surface Temp. (experimental)
In the Settings window for Table Graph, type Surface Temp. (experimental) in the Label text field.
Middle Temp. (experimental)
1
In the Model Builder window, under Results > Optimized, Forward Model, and Experimental Data: T_surface and T_middle click Middle Temperature.
2
In the Settings window for Table Graph, type Middle Temp. (experimental) in the Label text field.
Surface Temp. (forward model)
1
In the Model Builder window, right-click Optimized, Forward Model, and Experimental Data: T_surface and T_middle and choose Point Graph.
2
3
In the Settings window for Point Graph, type Surface Temp. (forward model) in the Label text field.
4
Locate the Data section. From the Dataset list, choose Study 1: Forward Model (Initial-Value Based)/Solution 1 (sol1).
5
6
Locate the y-Axis Data section. In the Expression text field, type T.
7
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dashed.
8
From the Color list, choose Magenta.
9
Click to expand the Legends section. Select the Show legends checkbox.
10
Find the Include subsection. Select the Label checkbox.
11
Clear the Point checkbox.
12
Clear the Solution checkbox.
Surface Temp. (optimized model)
1
Right-click Surface Temp. (forward model) and choose Duplicate.
2
In the Settings window for Point Graph, type Surface Temp. (optimized model) in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
4
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Solid.
Middle Temp. (forward model)
1
In the Model Builder window, right-click Surface Temp. (forward model) and choose Duplicate.
2
In the Settings window for Point Graph, type Middle Temp. (forward model) in the Label text field.
3
Locate the Selection section. Click to select the  Activate Selection toggle button.
4
5
Locate the Coloring and Style section. From the Color list, choose Black.
Middle Temp. (optimized model)
1
Right-click Middle Temp. (forward model) and choose Duplicate.
2
In the Settings window for Point Graph, locate the Data section.
3
From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
4
In the Label text field, type Middle Temp. (optimized model).
5
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Solid.
Optimized, Forward Model, and Experimental Data: T_surface and T_middle
1
In the Model Builder window, click Optimized, Forward Model, and Experimental Data: T_surface and T_middle.
2
In the Settings window for 1D Plot Group, locate the Plot Settings section.
3
Clear the Two y-axes checkbox.
4
In the y-axis label text field, type Temperature (K).
5
Locate the Legend section. From the Layout list, choose Inside graph axis area.
6
From the Position list, choose Lower right.
7
In the Number of columns text field, type 2.
8
In the Maximum relative width text field, type 1.
Middle Temp. (experimental)
1
In the Model Builder window, click Middle Temp. (experimental).
2
Drag and drop below Surface Temp. (optimized model).
Experimental Data 2
In the Model Builder window, expand the Results > Experimental Data 2 node.
Middle Temperature, Surface Temperature
1
In the Model Builder window, under Results > Experimental Data 2, Ctrl-click to select Surface Temperature and Middle Temperature.
2
Optimized, Forward Model, and Experimental Data: Y and T_center
1
In the Model Builder window, under Results click Experimental Data 2.
2
In the Settings window for 1D Plot Group, type Optimized, Forward Model, and Experimental Data: Y and T_center in the Label text field.
3
Locate the Legend section. From the Layout list, choose Outside graph axis area.
4
From the Position list, choose Bottom.
5
In the Number of rows text field, type 2.
Normalized Mass (experimental)
1
In the Model Builder window, under Results > Optimized, Forward Model, and Experimental Data: Y and T_center click Normalized Solid Mass.
2
In the Settings window for Table Graph, type Normalized Mass (experimental) in the Label text field.
Center Temp. (experimental)
1
In the Model Builder window, under Results > Optimized, Forward Model, and Experimental Data: Y and T_center click Center Temperature.
2
In the Settings window for Table Graph, type Center Temp. (experimental) in the Label text field.
Optimized, Forward Model, and Experimental Data: Y and T_center
Right-click Results > Optimized, Forward Model, and Experimental Data: Y and T_center > Center Temp. (experimental) and choose Global.
Normalized Mass (forward model)
1
In the Settings window for Global, type Normalized Mass (forward model) in the Label text field.
2
Locate the Data section. From the Dataset list, choose Study 1: Forward Model (Initial-Value Based)/Solution 1 (sol1).
3
Click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Definitions > domY - Domain Probe Y - 1.
4
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dashed.
5
From the Color list, choose Red.
6
Click to expand the Legends section. Find the Include subsection. Select the Label checkbox.
7
Clear the Solution checkbox.
8
Clear the Description checkbox.
Normalized Mass (optimized model)
1
Right-click Normalized Mass (forward model) and choose Duplicate.
2
In the Settings window for Global, type Normalized Mass (optimized model) in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
4
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Solid.
Center Temp. (forward model)
1
In the Model Builder window, right-click Optimized, Forward Model, and Experimental Data: Y and T_center and choose Point Graph.
2
In the Settings window for Point Graph, type Center Temp. (forward model) in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1: Forward Model (Initial-Value Based)/Solution 1 (sol1).
4
Locate the y-Axis Data section. In the Expression text field, type T.
5
Locate the y-Axis section. Select the Plot on secondary y-axis checkbox.
6
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dashed.
7
From the Color list, choose Blue.
8
Locate the Legends section. Select the Show legends checkbox.
9
Find the Include subsection. Select the Label checkbox.
10
Clear the Point checkbox.
11
Clear the Solution checkbox.
12
Center Temp. (optimized model)
1
Right-click Center Temp. (forward model) and choose Duplicate.
2
In the Settings window for Point Graph, type Center Temp. (optimized model) in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
4
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Solid.
5
In the Optimized, Forward Model, and Experimental Data: Y and T_center toolbar, click  Plot.
Turn on manual time stepping to avoid that interpolation errors affect the objective function. Using manual time stepping will increase the computational time significantly. With automatic time stepping the computational error may vary between iterations, affecting the value of the objective function, and thus the final results. If the time list is instead fixed, as with manual time stepping, the interpolation error does not change between iterations, and the objective function is evaluated without this variation. This means that the optimized values will deviate less between runs. When using a manual time step the relative and absolute tolerances are not used. A step convergence study should thus be performed. Such a study consists of re-computing the model with different manual time steps until the results from the different runs do not vary significantly. For this model, a time step of 1 s will give accurate results.
During optimization, a table will be available to inspect the values for the control parameters and objective function for each model evaluation.
Study 2: Parameter Estimation
General Optimization
1
In the Model Builder window, under Study 2: Parameter Estimation click General Optimization.
2
In the Settings window for General Optimization, locate the Output section.
3
Select the Plot checkbox.
4
Solve the optimization problem.
5
In the Study toolbar, click  Compute.
Results
Optimized, Forward Model, and Experimental Data: T_surface and T_middle
1
In the Model Builder window, under Results click Optimized, Forward Model, and Experimental Data: T_surface and T_middle.
2
In the Optimized, Forward Model, and Experimental Data: T_surface and T_middle toolbar, click  Plot.
3
Click the  Zoom Extents button in the Graphics toolbar.
Optimized, Forward Model, and Experimental Data: Y and T_center
1
In the Model Builder window, click Optimized, Forward Model, and Experimental Data: Y and T_center.
2
In the Optimized, Forward Model, and Experimental Data: Y and T_center toolbar, click  Plot.
3
Click the  Zoom Extents button in the Graphics toolbar.
The optimization has finished. The optimized values for the physical parameters can be derived using a Global Evaluation feature.
Optimized parameter values
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, type Optimized parameter values in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
4
From the Time selection list, choose Last.
Global Evaluation 1
1
Right-click Optimized parameter values and choose Global Evaluation.
2
In the Settings window for Global Evaluation, click Add Expression in the upper-right corner of the Expressions section. From the menu, choose Global definitions > Parameters > A_is - Frequency factor w -> is (intermediate solid) - 1/s.
3
Click Add Expression in the upper-right corner of the Expressions section. From the menu, choose Global definitions > Parameters > DH_c - Heat of reaction intermediate solid -> char - J/kg.
4
Click Add Expression in the upper-right corner of the Expressions section. From the menu, choose Global definitions > Parameters > DH_t - Heat of reaction wood -> tar - J/kg.
5
Click Add Expression in the upper-right corner of the Expressions section. From the menu, choose Global definitions > Parameters > hconv - External convective heat transfer coefficient - W/(m²·K).
6
In the Optimized parameter values toolbar, click  Evaluate.
Objectives
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, type Objectives in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2: Parameter Estimation/Solution 2 (sol2).
Global Evaluation 1
1
Right-click Objectives and choose Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Expressions section.
3
4
Locate the Data Series Operation section. From the Transformation list, choose Integral.
5
In the Objectives toolbar, click  Evaluate.
Objectives
Go to the Objectives window.
One can see that the objective related to the normalized solid mass is the largest, which explains why this property matches the data particularly well. A lower value for Tdelta can be used to change this balance.
Results
Set up a 3D plot to illustrate how the amounts of the three solid species change during the pyrolysis process. When done, we will have set up Figure 6 in the model documentation.
Solid Species
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Solid Species in the Label text field.
3
Locate the Color Legend section. From the Position list, choose Right double.
4
Click to expand the Plot Array section. Select the Enable checkbox.
5
From the Array shape list, choose Square.
6
From the Array plane list, choose xz.
Mirror 2D 1
Prepare a dataset that illustrates a half sphere.
1
In the Results toolbar, click  More Datasets and choose Mirror 2D.
2
In the Settings window for Mirror 2D, locate the Axis Data section.
3
In row Point 2, set R to 1.
4
In row Point 2, set Z to 0.
5
Half Sphere
1
In the Results toolbar, click  More Datasets and choose Revolution 2D.
2
In the Settings window for Revolution 2D, type Half Sphere in the Label text field.
3
Locate the Data section. From the Dataset list, choose Mirror 2D 1.
4
Click to expand the Revolution Layers section. In the Revolution angle text field, type 180.
5
150s Wood
1
In the Model Builder window, right-click Solid Species and choose Surface.
2
In the Settings window for Surface, type 150s Wood in the Label text field.
3
Locate the Data section. From the Dataset list, choose Half Sphere.
4
From the Time (s) list, choose Interpolation.
5
In the Time text field, type 150.
6
Locate the Expression section. In the Expression text field, type rho_w/rho_w_init.
There are many color tables to choose from. There is also the possibility to create your own. Follow these steps to create three new color tables that will simplify the interpretation of the plot.
7
Click the  Show More Options button in the Model Builder toolbar.
8
In the Show More Options dialog, select Results > Color Tables in the tree.
9
10
Wood
1
In the Model Builder window, under Results right-click Color Tables and choose Color Table.
2
In the Settings window for Color Table, type Wood in the Label text field.
3
Locate the Definition section. In the table, enter the following settings:
Intermediate
1
Right-click Color Tables and choose Color Table.
2
In the Settings window for Color Table, type Intermediate in the Label text field.
3
Locate the Definition section. In the table, enter the following settings:
Char
1
Right-click Color Tables and choose Color Table.
2
In the Settings window for Color Table, type Char in the Label text field.
3
Locate the Definition section. In the table, enter the following settings:
150s Wood
1
In the Model Builder window, under Results > Solid Species click 150s Wood.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Wood.
4
From the Color table transformation list, choose Reverse.
5
From the Color table type list, choose Discrete.
150s Intermediate Solid
1
Right-click 150s Wood and choose Duplicate.
2
In the Settings window for Surface, type 150s Intermediate Solid in the Label text field.
3
Locate the Expression section. In the Expression text field, type rho_is/rho_w_init.
4
Locate the Coloring and Style section. From the Color table list, choose Intermediate.
150s Char
1
Right-click 150s Intermediate Solid and choose Duplicate.
2
In the Settings window for Surface, type 150s Char in the Label text field.
3
Locate the Expression section. In the Expression text field, type rho_c/rho_w_init.
4
Locate the Coloring and Style section. From the Color table list, choose Char.
5
Click to expand the Plot Array section. Select the Manual indexing checkbox.
6
In the Column index text field, type 2.
270s Wood
1
In the Model Builder window, right-click 150s Wood and choose Duplicate.
2
In the Settings window for Surface, type 270s Wood in the Label text field.
3
Locate the Data section. In the Time text field, type 270.
4
Click to expand the Inherit Style section. From the Plot list, choose 150s Wood.
5
Locate the Plot Array section. Select the Manual indexing checkbox.
6
In the Row index text field, type 1.
270s Intermediate Solid
1
In the Model Builder window, right-click 150s Intermediate Solid and choose Duplicate.
2
In the Settings window for Surface, type 270s Intermediate Solid in the Label text field.
3
Locate the Data section. In the Time text field, type 270.
4
Locate the Inherit Style section. From the Plot list, choose 150s Intermediate Solid.
5
Locate the Plot Array section. Select the Manual indexing checkbox.
6
In the Row index text field, type 1.
7
In the Column index text field, type 1.
270s Char
1
In the Model Builder window, right-click 150s Char and choose Duplicate.
2
In the Settings window for Surface, type 270s Char in the Label text field.
3
Locate the Data section. In the Time text field, type 270.
4
Locate the Inherit Style section. From the Plot list, choose 150s Char.
5
Locate the Plot Array section. In the Row index text field, type 1.
410 s Wood
1
In the Model Builder window, right-click 270s Wood and choose Duplicate.
2
In the Settings window for Surface, type 410 s Wood in the Label text field.
3
Locate the Data section. In the Time text field, type 410.
4
Locate the Plot Array section. In the Row index text field, type 2.
410 s Intermediate Solid
1
In the Model Builder window, right-click 270s Intermediate Solid and choose Duplicate.
2
In the Settings window for Surface, type 410 s Intermediate Solid in the Label text field.
3
Locate the Data section. In the Time text field, type 410.
4
Locate the Plot Array section. In the Row index text field, type 2.
410 s Char
1
In the Model Builder window, right-click 270s Char and choose Duplicate.
2
In the Settings window for Surface, type 410 s Char in the Label text field.
3
Locate the Data section. In the Time text field, type 410.
4
Locate the Plot Array section. In the Row index text field, type 2.
Solid Species
1
In the Model Builder window, click Solid Species.
2
In the Settings window for 3D Plot Group, click to expand the Title section.
3
From the Title type list, choose None.
Wood
1
Right-click Solid Species and choose Annotation.
2
In the Settings window for Annotation, type Wood in the Label text field.
3
Locate the Data section. From the Dataset list, choose Half Sphere.
4
Locate the Annotation section. In the Text text field, type $\frac{\rho_{\omega}}{\rho_{\omega,0}}$.
5
Select the LaTeX markup checkbox.
6
Locate the Position section. In the z text field, type -0.02.
7
Locate the Coloring and Style section. Clear the Show point checkbox.
8
From the Anchor point list, choose Center.
9
Click to expand the Plot Array section. Select the Manual indexing checkbox.
Intermediate Solid
1
Right-click Wood and choose Duplicate.
2
In the Settings window for Annotation, type Intermediate Solid in the Label text field.
3
Locate the Annotation section. In the Text text field, type $\frac{\rho_{is}}{\rho_{\omega,0}}$.
4
Locate the Plot Array section. In the Column index text field, type 1.
Char
1
Right-click Intermediate Solid and choose Duplicate.
2
In the Settings window for Annotation, type Char in the Label text field.
3
Locate the Annotation section. In the Text text field, type $\frac{\rho_{c}}{\rho_{\omega,0}}$.
4
Locate the Plot Array section. In the Column index text field, type 2.
150 s
1
In the Model Builder window, right-click Wood and choose Duplicate.
2
In the Settings window for Annotation, type 150 s in the Label text field.
3
Locate the Annotation section. In the Text text field, type 150 s.
4
Locate the Position section. In the z text field, type 0.
5
In the x text field, type -0.025.
270 s
1
Right-click 150 s and choose Duplicate.
2
In the Settings window for Annotation, type 270 s in the Label text field.
3
Locate the Annotation section. In the Text text field, type 270 s.
4
Locate the Plot Array section. In the Row index text field, type 1.
410 s
1
Right-click 270 s and choose Duplicate.
2
In the Settings window for Annotation, type 410 s in the Label text field.
3
Locate the Annotation section. In the Text text field, type 410 s.
4
Locate the Plot Array section. In the Row index text field, type 2.
Solid Species
1
In the Model Builder window, click Solid Species.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
From the View list, choose View 3D 2.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
In the Solid Species toolbar, click  Plot.
T, Qmass and Q at 150 s
Having illustrated the progress of the solid reactions, look at the temperature, mass source, and heat source. Follow the steps below to set up Figure 7Figure 9 in the model documentation.
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type T, Qmass and Q at 150 s in the Label text field.
3
Locate the Data section. From the Dataset list, choose Probe Solution 3 (sol1).
4
From the Time (s) list, choose Interpolation.
5
In the Time text field, type 150.
6
Click to expand the Title section. From the Title type list, choose None.
7
Locate the Color Legend section. Select the Show units checkbox.
8
From the Position list, choose Right double.
9
Click to expand the Plot Array section. Select the Enable checkbox.
10
From the Array shape list, choose Square.
11
From the Order list, choose Column-major.
Temperature
1
Right-click T, Qmass and Q at 150 s and choose Annotation.
2
In the Settings window for Annotation, type Temperature in the Label text field.
3
Locate the Annotation section. In the Text text field, type $T$.
4
Select the LaTeX markup checkbox.
5
Locate the Coloring and Style section. Clear the Show point checkbox.
6
From the Anchor point list, choose Center.
7
Click to expand the Plot Array section. Select the Manual indexing checkbox.
8
Locate the Position section. In the R text field, type 0.005.
9
In the Z text field, type 0.0135.
Mass Source
1
Right-click Temperature and choose Duplicate.
2
In the Settings window for Annotation, type Mass Source in the Label text field.
3
Locate the Annotation section. In the Text text field, type $Q_{mass}$.
4
Locate the Plot Array section. In the Row index text field, type 1.
Heat Source
1
Right-click Mass Source and choose Duplicate.
2
In the Settings window for Annotation, type Heat Source in the Label text field.
3
Locate the Annotation section. In the Text text field, type $Q$.
4
Locate the Plot Array section. In the Row index text field, type 2.
T
1
In the Model Builder window, right-click T, Qmass and Q at 150 s and choose Surface.
2
In the Settings window for Surface, type T in the Label text field.
3
Locate the Expression section. In the Expression text field, type T.
4
Locate the Coloring and Style section. From the Color table list, choose HeatCameraLight.
5
Click to expand the Plot Array section. Select the Manual indexing checkbox.
dl.Qm
1
Right-click T and choose Duplicate.
2
In the Settings window for Surface, type dl.Qm in the Label text field.
3
Locate the Expression section. In the Expression text field, type dl.Qm.
4
Locate the Coloring and Style section. From the Color table list, choose Crocus.
5
From the Scale list, choose Linear symmetric.
6
Locate the Plot Array section. In the Row index text field, type 1.
Q
1
Right-click dl.Qm and choose Duplicate.
2
In the Settings window for Surface, type Q in the Label text field.
3
Locate the Expression section. In the Expression text field, type Q.
4
Locate the Coloring and Style section. From the Color table list, choose Avicularia.
5
Locate the Plot Array section. In the Row index text field, type 2.
T, Qmass and Q at 150 s
1
Click the  Show Grid button in the Graphics toolbar.
2
In the Model Builder window, click T, Qmass and Q at 150 s.
3
In the T, Qmass and Q at 150 s toolbar, click  Plot.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
In the T, Qmass and Q at 150 s toolbar, click  Plot.
T, Qmass and Q at 270 s
1
Right-click T, Qmass and Q at 150 s and choose Duplicate.
2
In the Settings window for 2D Plot Group, type T, Qmass and Q at 270 s in the Label text field.
3
Locate the Data section. In the Time text field, type 270.
4
In the T, Qmass and Q at 270 s toolbar, click  Plot.
5
Click the  Zoom Extents button in the Graphics toolbar.
6
In the T, Qmass and Q at 270 s toolbar, click  Plot.
T, Qmass and Q at 433 s
1
Right-click T, Qmass and Q at 270 s and choose Duplicate.
2
In the Settings window for 2D Plot Group, type T, Qmass and Q at 433 s in the Label text field.
3
Locate the Data section. In the Time text field, type 433.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
In the T, Qmass and Q at 433 s toolbar, click  Plot.
Pressure, velocity, porosity, and normalized solid mass at 270 s
Now plot the relative pressure, the total Darcy velocity magnitude, porosity, normalized solid mass, and the total Darcy velocity field. This gives Figure 10 in the model documentation.
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Pressure, velocity, porosity, and normalized solid mass at 270 s in the Label text field.
3
Locate the Data section. From the Dataset list, choose Probe Solution 3 (sol1).
4
From the Time (s) list, choose Interpolation.
5
In the Time text field, type 270.
6
Locate the Title section. From the Title type list, choose None.
7
Locate the Color Legend section. Select the Show units checkbox.
8
From the Position list, choose Right double.
9
Locate the Plot Array section. Select the Enable checkbox.
10
From the Array shape list, choose Square.
Relative Pressure
1
Right-click Pressure, velocity, porosity, and normalized solid mass at 270 s and choose Annotation.
2
In the Settings window for Annotation, type Relative Pressure in the Label text field.
3
Locate the Annotation section. In the Text text field, type $p/p_{ref}$.
4
Select the LaTeX markup checkbox.
5
Locate the Position section. In the R text field, type 0.005.
6
In the Z text field, type 0.0135.
7
Locate the Coloring and Style section. Clear the Show point checkbox.
8
From the Anchor point list, choose Center.
9
Locate the Plot Array section. Select the Manual indexing checkbox.
Total Darcy Velocity Magnitude
1
Right-click Relative Pressure and choose Duplicate.
2
In the Settings window for Annotation, type Total Darcy Velocity Magnitude in the Label text field.
3
Locate the Annotation section. In the Text text field, type $U$.
4
Locate the Plot Array section. In the Column index text field, type 1.
Porosity
1
Right-click Total Darcy Velocity Magnitude and choose Duplicate.
2
In the Settings window for Annotation, type Porosity in the Label text field.
3
Locate the Annotation section. In the Text text field, type $\epsilon$.
4
Locate the Plot Array section. In the Row index text field, type 1.
5
In the Column index text field, type 0.
Normalized Solid Mass
1
Right-click Porosity and choose Duplicate.
2
In the Settings window for Annotation, type Normalized Solid Mass in the Label text field.
3
Locate the Annotation section. In the Text text field, type $Y$.
4
Locate the Plot Array section. In the Column index text field, type 1.
dl.pA/dl.pref
1
In the Model Builder window, right-click Pressure, velocity, porosity, and normalized solid mass at 270 s and choose Surface.
2
In the Settings window for Surface, type dl.pA/dl.pref in the Label text field.
3
Locate the Expression section. In the Expression text field, type dl.pA/dl.pref.
4
Locate the Coloring and Style section. From the Color table list, choose Iodinea.
5
Locate the Plot Array section. Select the Manual indexing checkbox.
dl.U
1
Right-click dl.pA/dl.pref and choose Duplicate.
2
In the Settings window for Surface, type dl.U in the Label text field.
3
Locate the Expression section. In the Expression text field, type dl.U.
4
Locate the Coloring and Style section. From the Color table list, choose Acanthaster.
5
Locate the Plot Array section. In the Column index text field, type 1.
dl.u
1
In the Model Builder window, right-click Pressure, velocity, porosity, and normalized solid mass at 270 s and choose Arrow Surface.
2
In the Settings window for Arrow Surface, type dl.u in the Label text field.
3
Locate the Expression section. In the R-component text field, type dl.u.
4
In the Z-component text field, type dl.w.
5
Locate the Arrow Positioning section. Find the R grid points subsection. In the Points text field, type 10.
6
Find the Z grid points subsection. In the Points text field, type 10.
7
Locate the Coloring and Style section. From the Arrow type list, choose Cone.
8
From the Arrow length list, choose Normalized.
9
From the Arrow base list, choose Center.
10
Select the Scale factor checkbox. In the associated text field, type 0.14.
11
From the Color list, choose Black.
12
Click to expand the Plot Array section. Select the Manual indexing checkbox.
13
In the Column index text field, type 1.
epsilon
1
In the Model Builder window, right-click dl.pA/dl.pref and choose Duplicate.
2
In the Settings window for Surface, type epsilon in the Label text field.
3
Locate the Expression section. In the Expression text field, type epsilon.
4
Locate the Coloring and Style section. From the Color table list, choose Metasepia.
5
From the Color table transformation list, choose Reverse.
6
Locate the Plot Array section. In the Row index text field, type 1.
Y
1
Right-click epsilon and choose Duplicate.
2
In the Settings window for Surface, type Y in the Label text field.
3
Locate the Expression section. In the Expression text field, type Y.
4
Locate the Coloring and Style section. From the Color table list, choose Cynanthus.
5
From the Color table transformation list, choose None.
6
Locate the Plot Array section. In the Column index text field, type 1.
Pressure, velocity, porosity, and normalized solid mass at 270 s
1
Click the  Zoom Extents button in the Graphics toolbar.
2
In the Model Builder window, click Pressure, velocity, porosity, and normalized solid mass at 270 s.
3
In the Pressure, velocity, porosity, and normalized solid mass at 270 s toolbar, click  Plot.