Department of Science and Technology
Institutionen för teknik och naturvetenskap
Linköpings Universitet
Linköpings Universitet
Examensarbete
LITH-ITN-MT-EX--05/015--SE
Monte Carlo Simulation of Light
Scattering in Paper
Ronnie Dahlgren
LITH-ITN-MT-EX--05/015--SE
Monte Carlo Simulation of Light
Scattering in Paper
Examensarbete utfört i medieteknik
vid Linköpings Tekniska Högskola, Campus
Norrköping
Ronnie Dahlgren
Handledare Ludovic Coppel
Examinator Li Yang
Norrköping 2005-02-14
Rapporttyp Report category Examensarbete B-uppsats C-uppsats D-uppsats _ ________________ Språk Language Svenska/Swedish Engelska/English _ ________________ Titel Title Författare Author Sammanfattning Abstract ISBN _____________________________________________________ ISRN _________________________________________________________________
Serietitel och serienummer ISSN
Title of series, numbering ___________________________________
Nyckelord
Datum
Date
URL för elektronisk version
Avdelning, Institution
Division, Department
Institutionen för teknik och naturvetenskap Department of Science and Technology
2005-02-14
x
x
LITH-ITN-MT-EX--05/015--SE
Monte Carlo Simulation of Light Scattering in Paper
Ronnie Dahlgren
Paper is a very complex optical material. Analytical models explaining some of the optical properties of paper exist, but they often rely on bold simplifications. Monte Carlo simulation models are less
constrained and allow for a greater degree of complexity. Grace is a three-dimensional light scattering simulation tool for paper, previously implemented in Matlab. During this project, the basesheet model was implemented in C++. This model simulates a layer containing a network of wood fibers and filler material. The new implementation makes simulations much faster. In addition, some new features and enhancements were developed. Wavelength dependency of parameters and fluorescence were
implemented and tested. A problem with Grace was that there was no simple way to calculate the grammage of the simulated paper. A method to analytically determine grammage was developed so that the user has complete control over the grammage of the fiber network. Modifications were made to improve light scattering at the pore boundaries. A new feature was also added to study how the fibers’ and fillers’ geometry affects the light scattering at the paper surface.
Upphovsrätt
Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –
under en längre tid från publiceringsdatum under förutsättning att inga
extra-ordinära omständigheter uppstår.
Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,
skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för
ickekommersiell forskning och för undervisning. Överföring av upphovsrätten
vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av
dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,
säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ
art.
Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i
den omfattning som god sed kräver vid användning av dokumentet på ovan
beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan
form eller i sådant sammanhang som är kränkande för upphovsmannens litterära
eller konstnärliga anseende eller egenart.
För ytterligare information om Linköping University Electronic Press se
förlagets hemsida
http://www.ep.liu.se/
Copyright
The publishers will keep this document online on the Internet - or its possible
replacement - for a considerable time from the date of publication barring
exceptional circumstances.
The online availability of the document implies a permanent permission for
anyone to read, to download, to print out single copies for your own use and to
use it unchanged for any non-commercial research and educational purpose.
Subsequent transfers of copyright cannot revoke this permission. All other uses
of the document are conditional on the consent of the copyright owner. The
publisher has taken technical and administrative measures to assure authenticity,
security and accessibility.
According to intellectual property law the author has the right to be
mentioned when his/her work is accessed as described above and to be protected
against infringement.
For additional information about the Linköping University Electronic Press
and its procedures for publication and for assurance of document integrity,
please refer to its WWW home page:
http://www.ep.liu.se/
Abstract
Paper is a very complex optical material. Analytical models explaining some of the optical properties of paper exist, but they often rely on bold simplifica-tions. Monte Carlo simulation models are less constrained and allow for a greater degree of complexity. Grace is a three-dimensional light scattering simulation tool for paper, previously implemented in Matlab. During this project, the base-sheet model was implemented in C++. This model simulates a layer containing a network of wood fibers and filler material. The new implementation makes simulations much faster. In addition, some new features and enhancements were developed. Wavelength dependency of parameters and fluorescence were imple-mented and tested. A problem with Grace was that there was no simple way to calculate the grammage of the simulated paper. A method to analytically determine grammage was developed so that the user has complete control over the grammage of the fiber network. Modifications were made to improve light scattering at the pore boundaries. A new feature was also added to study how the fibers’ and fillers’ geometry affects the light scattering at the paper surface.
Acknowledgements
I wish to thank my tutor Ludovic Coppel for the opportunity to do my master’s project at Acreo and for his superb guidance.
My gratitude also to Dr. Li Yang at Link¨oping University for being my examiner.
My thanks as well to all friends who helped me with linguistic revisions and valuable comments on my report.
Contents
1 Introduction 1 1.1 Background . . . 1 1.2 Project Objective . . . 1 1.3 Method . . . 2 1.4 Thesis Overview . . . 2 2 Paper 3 2.1 Raw Material . . . 3 2.2 Papermaking . . . 5 2.3 Optical Properties . . . 53 Models in Paper Optics 7 3.1 Kubelka-Munk . . . 7
3.2 DORT . . . 7
3.3 Monte Carlo . . . 8
4 Grace 11 4.1 Overview . . . 11
4.2 On the Notation used in the Thesis . . . 11
4.3 Light Description . . . 12
4.4 Surface Scattering . . . 14
4.4.1 Specular Deflection . . . 14
4.4.2 Reflection vs. Transmission . . . 15
4.4.3 Specular vs. Diffuse Deflection . . . 15
4.4.4 Diffuse Deflection . . . 16
4.5 Scattering in Homogeneous Layers . . . 16
4.6 Light Sources . . . 17
4.7 Detectors . . . 18
4.8 The Paper Structure . . . 19
4.9 The Basesheet Layer . . . 19
4.9.1 Fibers . . . 19
4.9.2 Pores and Fillers . . . 20
4.10 Boundary Conditions . . . 20
5 Implementation 23 5.1 Statistical Distributions . . . 23
CONTENTS
6 Improvements 27
6.1 Grammage . . . 27
6.1.1 Evalutation of Grammage by Simulation . . . 27
6.1.2 Analytical Determination of Grammage . . . 28
6.2 Fluorescence . . . 29
6.3 A New Pore Representation . . . 31
6.4 A Different Approach to Surface Scattering . . . 31
7 Results 33 7.1 Basesheet Simulation Example . . . 33
7.2 Simulation of Real Basesheets Made of Eucalyptus Fibers . . . . 33
7.3 Using Fiber Geometry for Surface Reflection . . . 36
7.4 New Pore Representation . . . 36
7.5 Fluorescence . . . 39
7.6 Speed and Memory Usage . . . 41
8 Conclusions and Future Work 43 8.1 Conclusion . . . 43
8.2 Suggestions for Future Work . . . 43
8.2.1 Performance . . . 44
8.2.2 Other Representations of Paper Components . . . 44
8.2.3 Static Fiber Network . . . 44
8.2.4 Semi-dynamic Fiber Network . . . 44 A Generating Random Numbers 49 B Path Length 51 C Grammage Calculation 53 D Pore Size Distribution 57 E Diffuse Scattering 61
Chapter 1
Introduction
1.1
Background
Paper is a very complex optical material. Analytical models explaining some of the optical properties of paper exist, but often oversimplify the problem. Monte Carlo simulation models are less constrained and allow for a greater degree of complexity.
Grace is the name of a computer software tool for simulation of light scat-tering in printed paper. It is being developed by Acreo AB for four paper companies in joint projects. The model was first presented to the public in 1999 at a conference [11]. It is a Monte Carlo simulation model that describes paper as a layered 3D structure and simulates the interaction of light with paper and ink. This simulation tool is also the main subject of this master’s thesis.
1.2
Project Objective
The Grace model was previously implemented in Matlab. Within the Paper Optics and Perception project at Acreo, Grace is now being ported to C++. The transition to C++ will improve speed and allow a more efficient use of the computer memory. Being an object oriented programming language, C++ also makes it easier to upgrade the code and continously make modifications and improvements. The first goal of the project was to implement what remained to be done in C++, namely the basesheet layer. This is the layer containing the fiber network with fillers and pores. The second goal of the project was to identify and implement some improvements in order to make the model more consistent with real paper.
The scope of this master’s project is limited to the basesheet layer. Only minor modifications and adaptations needed for the implementation were done in other parts of Grace. Existing functions were used for the actual light scat-tering and the optics theory behind them is covered in this thesis but there were no attempts to improve or modify these functions.
CHAPTER 1. INTRODUCTION
1.3
Method
First, literature and manuals covering the theories and principles of Grace were studied, followed by a study of the existing program. The existing model was then implemented and adapted to C++. During this work several ideas on how to improve the model emerged and some of these ideas were then implemented and tested.
1.4
Thesis Overview
Chapter 2 introduces some paper properties and discusses the model represen-tation. In chapter 3 some different models used in paper optics are presented. Chapter 4 describes Grace and the theory involved in detail. Chapter 5 discusses briefly the C++ implementation from a programmer’s point of view. Chapter 6 addresses the modifications that were implemented and chapter 7 presents the results obtained throughout the project. Chapter 8 concludes the thesis and suggests some future improvements. The different appendices provide the mathematical theory for the physical models developed in the thesis.
Chapter 2
Paper
This chapter will briefly discuss the properties of papers and explain some of the assumptions used in the model representing paper in Grace.
2.1
Raw Material
Fibers are long tubes made of cellulose. They are the main building component of trees so they need to be very strong to support the whole weight of the tree. At the same time they are hollow and act as a channel for water. The cavity in the middle of the fiber is called the lumen. The fiber wall size and the total fiber diameter vary with the season and this can be seen as annual rings in the wood.
In the paper making process, the fibers are damaged in different ways. The fibers are cut into shorter lengths and studies show that the distribution of fiber lengths in paper can be approximated as log-normal [4], see figure 2.1. Some fibers also tend to collapse [6] due to, for example, beating of the pulp in the paper-making process or when the paper is put under pressure. Figure 2.2 shows some different possibilities for fiber collapse. In Grace, some simplifications are done. The fibers are modeled as cylinders with an elliptic cross-section, where the ellipticity represents different degrees of collapse. Even though the surface of an ellipse is smooth, the rough surfaces, as seen on the right in figure 2.2, are achieved by introducing a surface microroughness, which adds a degree of randomness to the light scattering. Chapter 4 will discuss how these different properties of the fiber are represented in Grace.
The bleaching is not always sufficient to remove the yellowish tint of the paper so paper often contains an amount of filler material with pigments and fluorescent material that increase the brightness and whiteness. Filler normally gathers as clusters in the paper. In Grace these clusters are modeled as ellipsoids and different parameters control the scattering of light at the surfaces and within the material. Fillers added to paper are often fluorescent, absorbing light at one wavelength and emitting at another. How this behavior can be simulated in Grace is discussed in chapter 6.
CHAPTER 2. PAPER
fiber length (µm)
probability density
Log normal distribution
0 2000
Figure 2.1: The log normal distribution is often used to estimate different geo-metrical parameters, such as fiber and pore sizes. [4]
Figure 2.2: Different types of fiber collapse. [6]
CHAPTER 2. PAPER
2.2
Papermaking
First, the fibers in the wood need to be extracted and separated from each other and this is done in the pulping process, either chemically or mechanically but both processes damage some of the fibers. Chemical processing tends to produce longer fibers but also dissolve some of the fiber wall.
The resulting fiber mass, called pulp, can be further refined and often mate-rials can be added. The pulp is then formed into a paper sheet which is pressed and dried. During this process, some fibers collapse, as seen in the previous section.
The paper can also be calendered, a process in which it is pressed very hard. This gives the paper a flatter and glossier surface. Since the fibers are pressed together and thus come more in contact with each other, the total amount of optical surfaces where scattering can occur will decrease and this can sometimes make the paper more transparent [20].
2.3
Optical Properties
There are a number of standard measures that quantify the optical properties of paper.
Brightness is defined as the reflectance from an infinitely thick paper sheet with light at a wavelength of 457 nm.
Whiteness is a measure of how white the paper is perceived to be. Brightness tells us how much light a surface emits for one single wavelength, while whiteness is a color value that represents how much all wavelengths are reflected by the paper.
Opacity is a measure for the paper’s ability to resist passage of light. It is defined as the ratio between the reflectance of a single sheet (R1) and that of
an infinite number of sheets (R∞) at wavelength 557 nm:
Opacity = R1 R∞
(2.1) The bi-directional scattering distribution function (BSDF) is a way of mea-suring the scattering as a function of angle, both in transmittance and in reflec-tion. It is defined as:
BSDF = dPs/dw P0cos θs
(2.2) where dPs is the power of the light scattered into a solid angle dw, P0 is the
incident power and θs is the polar scattering angle.
Total reflectance is the ratio of the incident intensity of light that is reflected from the paper. The total transmittance is the proportion of the incident in-tensity of light that is transmitted through the paper.
The total reflectance alone does not tell us much about how the paper is perceived by the eye. In this thesis the total reflectance will however be used for comparisons with real paper. This is because it is very easy to measure in simulations.
CHAPTER 2. PAPER
In this thesis, some simulations will be compared to measurements done with the Perkin-Elmer spectrophotometer ([2] and [3]). It is a commonly used instrument that measures the total spectral reflectance or the transmittance with light at an 8◦ incident angle.
One of the most commonly used instruments is the Lorentzen & Wettre El-repho spectrophotometer which measures the reflectance from paper with diffuse incoming light. This instrument can also be simulated but is more sensitive to noise in simulations. It requires long simulation times, which is why it won’t be used in this thesis.
Chapter 3
Models in Paper Optics
Since paper is such a complex material, some simplifications have to be made when studying paper from an optical point of view. Even if we knew the ex-act geometry of the components that make up the paper sheet it would not be possible to find a solution to the Maxwell equations describing the light interac-tion with the material. Until the end of the 19th century when Lord Rayleigh published his classic paper explaining why the sky is blue [22], there were no established theories for light scattering. Rayleigh explained light scattering at a molecular level and later, in 1908, Mie explained light scattering by larger particles [17]. However, these theories only deal with scattering from one single particle so another theory that deals with multiple scattering was needed. One of the theories dealing with multiple scattering is known as radiative transfer theory. The Kubelka-Munk and Discrete Ordinate Radiative Transfer (DORT) theories presented below are special cases of radiative transfer theory. Another approach for multiple scattering analysis is the use of so called Monte Carlo methods. This technique is old, but has gained in status only during the past decades, stimulated by the progress in modern computers.
3.1
Kubelka-Munk
The Kubelka-Munk theory was developed in the 1930s [15] for the study of light propagation in colorant layers. It is a two-flux approximation used to calculate the total diffuse reflection from a material and it is still predominant in the field of paper optics today. It is based on some assumptions, such as the light is equally scattered in all directions and it only considers the average of all the fluxes towards the up- and down-hemispheres. These are quite bold simplifications, but the theory has proven to be useful and is often used because of its simplicity.
3.2
DORT
The Kubelka-Munk theory has often been extended and adapted to different situations and one of these extensions is Discrete Ordinate Radiative Transfer (DORT). It was applied on paper by Mudgett and Richards [18] in 1971 and it takes into account the angular distribution of the light. The two surrounding
CHAPTER 3. MODELS IN PAPER OPTICS
hemispheres are divided into several channels according to figure 3.1. Applied to paper, the Kubelka-Munk only allows us to study the total flux reflected from or transmitted through the paper, while DORT lets us study in which proportions the light will be scattered in different directions.
A problem with DORT models is that they often do not take into account the contribution of the paper surface to the light scattering. Currently, a model for surface reflections is being implemented in DORT2002 [5]. However, it still assumes a homogeneous bulk without any thickness variations.
x y z 0° 20° 40° 60° 80°
Figure 3.1: DORT models. The two hemispheres are divided into several chan-nels and the amount of light propagating in each channel is calculated.
3.3
Monte Carlo
The Monte Carlo approach to multiple scattering is what this thesis will focus on.
The city of Monte Carlo is well known for its casinos and this has given name to a computation method which imitates the random behavior of the games of chance in a casino. In this method, variables that are uncertain but have a known range of values with a probability distribution are generated.
For instance, the value of π can be approximated with a Monte Carlo ap-proach. Consider figure 3.2 below. If we generate uniform random values of x and y and count the points (x, y) that fall within the circle, multiply this num-ber by 4 and divide by the total numnum-ber of points we will approximate π. It is actually the surface of the circle that is approximated, from which π is deduced and the accuracy increases with the number of points generated.
To return to the case with light scattering, we can in the same way use this kind of method to approximate the properties of a multiple-scattering medium. This is done by making assumptions about the random behavior of light and by performing simulations. Now, instead of generating random points on a surface, random wave packets representing light are generated. It is simulated how they propagate through and stochastically interact with the paper. Probability dis-tributions are specified for the stochastic variables which determine the outcome for the different processes. The angle of a diffuse scattering is an example of 8
CHAPTER 3. MODELS IN PAPER OPTICS 0 1 0 1 x y
Figure 3.2: Estimating π with a Monte Carlo method. Points on the surface are generated randomly. The number of points within the circle will be proportional to the circle area, which can be used to approximate π.
such a variable. If enough wave packets are simulated we can approximate with increasing accuracy the general optical properties of the paper.
The disadvantage of the Monte Carlo approach is that it is necessary to do time-consuming simulations. During this master’s project, simulation times were in the range of several hours, the longest being the simulations to make the images in figure 7.4 on page 37 which took about 10 hours per image on a Pentium III 450 MHz. One advantage is however that there is no limit to the complexity of the model. If we for example want to know how much of the incoming light is reflected back from a paper sheet, the Kubelka-Munk approach would be to describe the paper by a mathematical model and analytically cal-culate the ratio reflection/transmission. In order to do this the model needs to be simple enough to allow the equations to be solved exactly. When using a Monte Carlo approach there is no such limitation on the complexity of the model, it is just a matter of having enough computer power to simulate enough wave packets to sufficiently minimize the statistical error.
Chapter 4
Grace
Grace is a computer simulation tool developed at Acreo AB1 and this chapter
will deal with its main principles. For the moment Grace is unfortunately not commercially available to the public, and most material that this thesis is built on is confidential. The theory explained in this chapter is covered in [12], which is a non-public report.
4.1
Overview
The user defines light sources, a paper structure and detectors. Using a Monte Carlo approach, the light is represented as wave packets that enter the paper structure and bounce around until they leave it in some direction or are absorbed somewhere. Different detectors then collect the scattered wave packets. Figure 4.1 shows an example of a paper structure with two layers in Grace.
The paper is assumed to be composed of different layers that are always sep-arated by surfaces. The layers can represent ink films, coatings and basesheets containing fibers, fillers and pores. Each layer has different parameters such as refractive indices, scattering- and absorption coefficients. At the top and the bottom of the paper structure are two semi-infinite layers that contain the light source and the detectors.
The light is modeled as many wave packets and the overall effect of a mul-titude of wave packets simulated in this manner will approximate the outcome of a real physical experiment. When different processes occur, the outcome is randomly determined with a probability corresponding to the physical rules. If for instance 90 % of the light should be refracted and 10 % reflected according to the Fresnel equations [21], then there will be a probability of 0.9 for refraction and 0.1 for reflection.
4.2
On the Notation used in the Thesis
In order to avoid misinterpretations some terms as they are used in this thesis are explained here.
A surface is the interface separating two layers in the paper structure.
CHAPTER 4. GRACE
Figure 4.1: Example of a paper model in Grace. The paper is composed of different layers and each layer is always contained between two surfaces.
When referring to a fiber, the fiber wall and the lumen are included. If only the fiber wall made of cellulose is intended, the term fiber wall is used.
The basesheet is made up of three different components: Fibers, fillers and pores. In this thesis they are referred to as paper components or simply compo-nents.
Sometimes the phrase paper object is used. This is the umbrella term for surfaces, layers, light sources and detectors.
Angles called θ and ϕ with different subscripts are used in the thesis. Figure 4.2 illustrates these. θ is the polar angle from the z-axis and ϕ is the azimuthal angle from the x-axis.
We sometimes encounter the phrase plane of incidence. This is the plane perpendicular to the surface and parallel to the direction vector of the wave packet, see figure 4.3.
4.3
Light Description
The wave packet is the representation of light in Grace. It can sometimes be seen as a particle and sometimes as a wave, depending on the process it is submitted to. Every wave packet contains information about its current state. This information is updated as the wave packet propagates through the paper structure. The state of a wave packet includes:
• Position vector • Direction vector • Polarization vector 12
CHAPTER 4. GRACE
x
y z
Figure 4.2: Visualization of the angles often used throughout the thesis.
n
Plane of incidence
Surface
CHAPTER 4. GRACE • Wavelength
The direction and polarization vectors are always normalized and orthogonal to each other. Each wave packet is linearly polarized in one direction. By generating a multitude of wave packets with random polarization directions the effect of unpolarized light can still be achieved.
4.4
Surface Scattering
This section discusses scattering on surfaces and the discussion is equally true for scattering on the surfaces of the paper components in the basesheet. Figure 4.4 illustrates the different scattering processes.
Specularly deflected wave packet
Diffusely scattered wave packet
In-material scattering Surface
Fiber
Figure 4.4: The light scattering processes. The wave packet can be specularly or diffusely deflected at surfaces. When propagating through a material it can be scattered or absorbed.
When a wave packet intercepts a surface, the two possible outcomes are: • A specular deflection occurs. The wave packet is reflected or refracted
according to the classic laws of optics as if the surface was perfectly flat. • A diffuse deflection occurs. The wave packet is reflected or refracted
in a random direction according to the Lambertian distribution. This simulates the effect of the microroughness of the surface.
4.4.1
Specular Deflection
If we assume that the surface is locally flat, the angle of specularly reflected and transmitted light can be determined with the well-known law of reflection (equation 4.1) and Snell’s law (equation 4.2), illustrated in figure 4.5, where θi, θrand θtare the angles of incidence, reflection and transmission, and ni and
ntare the refractive indices of the incident and transmitting media.
θi= θr (4.1)
nisin θi= ntsin θt (4.2)
CHAPTER 4. GRACE A specularly deflected wave packet keeps the polarization state relative to its direction vector, i.e. the polarization is rotated to the new traveling direction.
n
in
tFigure 4.5: An illustration of Snell’s law and the law of reflection (refractive indices ni< nt).
4.4.2
Reflection vs. Transmission
The probability for reflection and transmission depends on the polarization and is determined with the Fresnel equations [21]:
Rs= nicos θi− ntcos θt nicos θi+ ntcos θt 2 (4.3) Rp= ntcos θi− nicos θt ntcos θi+ nicos θt 2 (4.4) where Rsis the reflection coefficient for light with a polarization perpendicular
to the plane of incidence and Rp is the reflection coefficient for light with a
polarization parallel to the plane of incidence.
The probability for reflection is calculated by weighing together these two coefficients with respect to the wave packet’s polarization:
Pr[reflection] =qR2
p(ep· p)2+ Rs(es· p)2 (4.5)
where ep and es are the polarization unit vectors parallel and perpendicular
to the plane of incidence respectively and p is the wave packet’s polarization vector.
4.4.3
Specular vs. Diffuse Deflection
To determine whether the wave packet should be specularly or diffusely de-flected, one of the following two formulas is used [1]:
Pr[diffuse reflection] = 1 − e− 4πRq ni cos θi λ 2 (4.6) Pr[diffuse transmission] = 1 − e−
2πRq ni cos θi−nt cos θt
λ
2
CHAPTER 4. GRACE
where λ is the light wavelength and Rq is the root mean square roughness,
indicating how rough the surface is. See for instance [2] for further details about the surface microroughness.
4.4.4
Diffuse Deflection
If the wave packet is to be diffusely deflected, the new direction is generated using the Lambertian distribution. The angle of reflection θror transmission θt
is determined using the following expression (see appendix E): F−1(R) = arccos(1 − 2R)
2 (4.8) where R is a random uniform variable between 0 and 1. The azimuthal angle ϕ is random and uniform between 0-360 degrees.
In case of a diffuse reflection, the new polarization state is determined with the Q-matrix (see [8]). Due to lack of a better theory, the polarization state is kept the same relative to the traveling direction for diffuse transmission.
4.5
Scattering in Homogeneous Layers
When the wave packet propagates in a radiative transfer material, such as in a radiative transfer bulk or in the fiber wall, in-material scattering and absorption can occur. The wave packet will then be scattered or absorbed before it has reached a surface. This is controlled by three parameters, the scattering and absorption coefficients, and the asymmetry factor (σscat, σabs and g).
The propagating distance before these two processes occur are supposed to be exponentially distributed. This means that the probability of the future realization does not depend on the distance already traveled so far. The distance traveled before a wave packet is scattered or absorbed is known as the mean free path. Every time the wave packet propagates in a scattering medium, the mean free path is generated with the exponential distribution:
f (l) = σexte−σextl (4.9)
where l is the path length and σext = σscat+ σabs is the extinction coefficient
(see appendix B).
The cumulative distribution function is:
F (l) = 1 − e−σextl= 1 − e−(σscat+σabs)l (4.10)
Random mean free path lengths can be generated by taking the inverse cumulative function (see appendix A):
F−1(u) = − ln(1 − u) σscat+ σabs
(4.11) where u is a uniformly distributed random number, u ∈ U (0, 1).
If the distance for the wave packet to travel is greater than the generated distance, an absorption or scattering will occur, and it is easily demonstrated that the probabilities for either of these processes will be:
CHAPTER 4. GRACE Pr[scattering] = σscat σscat+ σabs (4.12) Pr[absorption] = σabs σscat+ σabs (4.13) If the wave packet is absorbed, it is removed from the simulation. If on the other hand a scattering occurs, a new direction is calculated. For this, the Henyey-Greenstein phase function (equation 4.14) [13] is used. Figure 4.6 shows a visualization of this function for g = 0.3.
Θ(θ) = 1 4π
1 − g2
(1 + g2− 2g cos θ)3/2 (4.14)
where θ is the polar angle with respect to the original propagation direction and g is the asymmetry factor. The asymmetry factor has values between -1 and 1. A value of -1 means that all wave packets are backscattered in the same direction they came from and a value of 1 that all wave packets continue in the same direction. A value of 0 gives a uniform distribution.
The Henyey-Greenstein phase function
0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340
Figure 4.6: Polar plot of the Henyey-Greenstein phase function, g = 0.3.
4.6
Light Sources
The initial distribution of the wave packets is determined by the light source. Three standard illuminations are available in Grace:
• Beam: A light source which, if placed distant from the paper, will produce an approximately parallel beam of wave packets.
• Lambertian: A light source that produces wave packets distributed as diffuse light.
• Elrepho: Emulates diffuse incoming light with a gloss trap like in the Elrepho spectrophotometer.
CHAPTER 4. GRACE
The light can either be uniformly distributed over the surface or concentrated in one point.
4.7
Detectors
Detectors can be of the following types: • Power detector
• ARS detector • ARS globe detector • Image detector • Elrepho detector
The simplest detector is the power detector. It measures the total per-centage of reflected, transmitted and absorbed wave packets. This detector is automatically implemented without having to be explicitly set up.
The angle resolved scattering (ARS) detector measures the reflectance or transmission at different angles. Several different apertures are positioned on a hemisphere above or below the sample, depending on if we want to measure light in reflection or transmission, see figure 4.7. The wave packets passing through every aperture are counted and the result can be viewed as a BSDF curve. The ARS globe detector [23] measures the angular distribution of light over the whole hemisphere and produces results like figure 7.4 on page 37.
The image detector emulates a simple camera with an aperture whose posi-tion and diameter can be chosen.
The Elrepho detector simulates the Elrepho instrument and gives reflectance factors comparable to what is experimentally obtained with that instrument.
0 0 0 x y z r Paper Light ray Apertures
Figure 4.7: Schematic view of the ARS detector. The light passing through each aperture is measured.
CHAPTER 4. GRACE
4.8
The Paper Structure
As seen in figure 4.1, the paper model is like a sandwich made up of different layers. A layer can represent a coating, a cellulose fiber network, a plastic film, air, an ink layer, etc. A layer must always be contained between two surfaces.
Surfaces are rectangular grids or height maps. In the simplest case the height map can be all zeros which gives a flat surface. The height map can also be measured from a real paper sample or generated using an algorithm.
There are several different ways to represent a surface in Grace, for example as an interpolated grid, a mesh of triangular polygons or a flat surface with mapped surface normals. Each method has its pros and cons, see [23] for more details.
The layers between the surfaces can be of two different kinds, radiative transfer bulks or basesheets. The radiative transfer bulks are homogenous layers that scatter and absorb light (see section 4.5).
The basesheets are layers containing fibers, pores and fillers. As the wave packet propagates in the basesheet, components are generated randomly accord-ing to probability distributions.
4.9
The Basesheet Layer
The basesheet is a layer containing fibers and fillers. The fiber network is not static with fibers in predetermined positions. A component distribution specifies the probabilities of finding a certain component at different depths in the paper. When a wave packet, propagating through a paper component, reaches the component’s outer surface, the next component is determined by this distribution. Figure 4.8 is an example of a component distribution for a basesheet consisting of only fibers and pores. This example was aquired through analysis of cross-sectional images of paper. [2]
In addition to the component distribution, a contact reduction factor is as-signed to every component type. This factor modifies the probability for contact between components of the same kind. For example, a contact reduction factor of 0.6 for fibers means that the probability for fiber-to-fiber contact is the fiber probability from the component distribution multiplied by 0.6. The probability for fiber-to-fiber contact could for instance increase when the paper is com-pressed and the contact reduction factor allows this to be taken into account.
4.9.1
Fibers
Fibers are implemented as a hollow cylinders, stretched out into an elliptic shape. In Grace the microroughness value is used to achieve the imperfections of the fibers giving unpredictable and random behavior of the wave packets.
Fibers are generated using several parameters (see figures 4.9 and 4.10): the polar angle of the fiber direction (θ), the azimuthal angle (ϕ), the length of the fiber (l), the fiber radius (r), the wall thickness (wt) i.e. the relative thickness
of the fiber wall compared to the minor axis of the elliptic cross-section, the cross-section ellipticity (e), the scattering coefficient in the fiber wall (σscat), the
CHAPTER 4. GRACE -40 -20 0 20 40 60 80 100 120 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Component probability Paper depth Component distribution Pore Fiber
Figure 4.8: Example of a component distribution.
factor for light scattering (g), the refractive index for the fiber wall and the refractive index for the lumen.
4.9.2
Pores and Fillers
Pores and fillers are modeled as outstretched spheres.
The different parameters describing a filler are: the polar angle of the filler direction (θ), the azimuthal angle (ϕ), the filler axis in the ϕ direction (a), the filler axis perpendicular to the ϕ direction (b), the cross-section ellipticity (e), the scattering coefficient in the filler (σscat), the absorption coefficient in the
filler (σabs), the Henyey-Greenstein asymmetry factor for light scattering (g)
and the refractive index.
The only difference between a pore and a filler is that pores have a refractive index equal to one, no scattering and no absorption.
4.10
Boundary Conditions
Since the area of the simulated paper is limited we have to decide what happens when a wave packet reaches the side boundary of the paper structure. The simulation only considers a limited area of the paper sample, while in reality it would often be bigger than this. It would, for instance, not be realistic to just remove wave packets from the simulation if they leave the paper structure. This would produce a dark band close to the edges since energy would leak out. In a real situation, the amount of light entering the simulation area would be approximately the same as the amount of light leaking out.
CHAPTER 4. GRACE
x
y
z
Le
ng
th
Figure 4.9: The geometric layout of a fiber.
e
e
r
r'
+
=
1
1
-r
w
tr'
Figure 4.10: Cross section of a fiber in the simulation model. The fiber is modeled as a hollow cylinder with an elliptic cross-section.
CHAPTER 4. GRACE
In Grace this is solved by bouncing the wave packets on the side boundaries. If we simulate an anisotropic material, i.e. a layer in which the light has a preferred traveling direction, a simple reflection according to the law of reflection would alter the wave packet’s traveling direction compared to the preferred direction. Figure 4.11 shows an anisotropic paper with a preferred direction. To preserve the direction according to the anisotropy of the material, the boundary conditions according to the dashed arrow have to be used. In this case, instead of only changing the sign of the direction cosines for x like a normal reflection, both the x and y directions are reversed.
x
y
z
Top view
Figure 4.11: Boundary conditions. An anisotropic paper sample viewed from above with a preferred traveling direction for light along the diagonal lines. In order to preserve the preferred direction, boundary conditions represented by the dashed arrow are used.
Chapter 5
Implementation
When porting the application from Matlab to C++, it is a good idea to make use of the structures that object oriented programming provides. The program structure is inspired by some basic ideas widely used in programming, known as design patterns [7]. These are basic ways of structuring object oriented code which have proven to give extendable and reusable code.
This chapter discusses the program structure from a programmer’s point of view. It assumes that the reader has some basic knowledge of object oriented programming. Users only interested in the paper optics perspective of the pro-gram can safely skip this chapter. The figures used are according to the UML standard [19].
5.1
Statistical Distributions
The basesheet is a layer with a dynamic network of fibers and fillers. At every position reached by the wave packet it is determined according to a component distribution which paper component it will encounter. When the decision is made, the parameters determining the component are generated according to statistical distributions. Exactly how the random parameters are generated from statistical distributions is explained in appendix A.
It is a good idea to have an abstract class representing statistical tions, to provide an easy access to random numbers following a given distribu-tion. Thus, there is no need to worry about which kind of random variable we are dealing with, whether it is Gaussian, elliptic distributed or has any other thinkable distribution. If we later want to extend the program by adding an-other kind of distribution, e.g. using exponentially distributed pore sizes, this can easily be done by adding a subclass representing the distribution without changing the code anywhere else. The bottom of figure 5.1 shows this structure with the abstract class Distribution and different subclasses representing differ-ent distributions. Note that a constant value can very well be represdiffer-ented in this way, so it is a good idea to use the Distribution class even for parameters that are constant in the current version of Grace but which might be changed to random variables in later versions. This adds no complexity to the program and the performance loss is hardly noticeable. As an example, in the previous version of Grace the fiber wall thickness was constant, but in the C++
imple-CHAPTER 5. IMPLEMENTATION
mentation it is a distribution. In the program interface it is still a constant and in the simulation a DistributionConstant class is used to represent it, but nothing prevents another kind of random variable being used. It is just a matter of adding this to the interface.
The DistributionTable class takes a custom probability density function as a parameter which allows just about any imaginable distribution to be used and this can replace other existing probability distributions. However, this class uses table lookups that are rather slow and it is better to use specialized classes that use faster analytical functions to generate the random numbers.
+GenerateRandomNumber() : double Distribution DistributionHenyeyGreenstein DistributionConstant DistributionElliptic DistributionTable +GenerateComponent() : Component -radius : Distribution -... ComponentGenerator 1 * +process(inout wp : WavePacket) -radius : double -... Component generates
Figure 5.1: Class diagram showing the relation Generator-Component. The gen-erator classes contain statistical distributions and use these to generate random components with fixed numbers.
5.2
Paper Components
The top of figure 5.1 and figure 5.2 show the general idea of the program struc-ture, which is an extended form of the class factory design pattern [7]. Every paper component is a subclass of the Component abstract superclass. The abstract Component class has a method that allows the wave packet to be processed and it is overloaded by every subclass. Thus accessing this method is always the same regardless of the type of component and the basesheet class is independant of the kind of components it contains, it just deals with Com-ponent superclasses. Thus, the program can easily be extended with new kinds of paper components without the need to change the basesheet class or other parts of the program.
Attentive readers might wonder why the pore is a second generation subclass under filler according to figure 5.2. This is because the pore is implemented as a kind of filler where the filler material is air. This was changed in a later improvement described in chapter 6.
Every paper component class has a class factory associated with it. In this 24
CHAPTER 5. IMPLEMENTATION Basesheet +process(inout wp : WavePacket) Component +GenerateComponent() : Component ComponentGenerator FiberGenerator FillerGenerator PoreGenerator Fiber Filler Pore generates generates generates 1 *
Figure 5.2: Class diagram showing the relation between the basesheet class and the component classes. The basesheet contains one generator class for every component type.
context, the word generator is used in all class factory names to distinguish them. All generator classes inherit from a superclass, the ComponentGenerator class. The generator classes deal with generating the paper component objects. They contain probability density functions while the Component classes contain fixed numbers. A basesheet that contains one single type of fiber would normally have one FiberGenerator class associated with it. Figure 5.3 illustrates how the basesheet layer processes a wave packet that encounters a fiber. When a wave packet encounters a fiber, the basesheet calls the GenerateComponent method of the generator class. This method generates a random fiber object according to the statistical distributions and returns it. The basesheet then only has to call the process method of the component which will propagate the wave packet through the fiber. When the wave packet leaves the fiber somewhere, the basesheet generates another component in the same manner and this continues until the wave packet leaves the basesheet layer.
If we want to add another kind of paper component, this can be done by cre-ating a class that inherits from Component and implements the process method. This class should represent a component and be able to process a wave packet propagating through it. A ComponentGenerator class must also be provided that generates all different variations of the paper component found in the base-sheet. Since the basesheet always works through the interface provided in the superclasses, no changes are needed in the Basesheet class.
CHAPTER 5. IMPLEMENTATION
Figure 5.3: Collaboration diagram showing an example of a fiber generation. The basesheet calls the function GenerateComponent and get a Fiber object in return. The function process propagates the wave packet in the fiber. The function returns when the wave packet leaves the fiber and the fiber is then deleted.
Chapter 6
Improvements
This chapters covers the improvements that were made to the model and the new features that were implemented during this master’s project.
6.1
Grammage
The grammage tells us how much each square meter of paper weighs. This will depend on parameters such as component size, shape and density. Since the components are randomly generated, determination of the grammage is not intuitive.
6.1.1
Evalutation of Grammage by Simulation
One way to evaluate the grammage of paper in a simulation is to set all re-fractive indices to one, permit no scattering and set a small absorption for the components. Since all the refractive indices are 1 and there is no scattering, the wave packets will propagate right through the paper without being deflected. Some wave packets will be absorbed in the paper and all the others will propa-gate in a straight line right through it. If the light source is set with an incoming angle of 0◦, the total distance traveled by each wave packet is equal to the paper thickness. The probability of absorption of a wave packet will increase with the distance traveled inside absorbing components. The probability of absorption can then be calculated with the following expression:
A(l) = 1 − e−αl (6.1) where α is the absorption coefficient and l is the distance traveled inside the absorbing component. For readers acquainted with the Kubelka-Munk theory it needs to be pointed out that the absorption coefficient is the parameter used as input to Grace and not the Kubelka-Munk parameter. Equation 6.1 is a direct result of equation 4.9 discussed in chapter 4. α is set as a parameter for the simulation, and the simulation will give A, so it is then possible to calculate l by inverting the equation above:
l = −ln(1 − A)
CHAPTER 6. IMPROVEMENTS
The grammage ω can now be calculated:
ω = lρ (6.3) where ρ is the density of the component and l is the value calculated above. If there are several paper components with different densities the grammage can be calculated by running one simulation for every component. The contribu-tion from every component to the grammage can be calculated by setting all absorption coefficients to zero except for the component being treated. The total grammage is then calculated by adding together the contributions from all components.
6.1.2
Analytical Determination of Grammage
It is often desirable to keep control over the grammage in the simulations. If the paper models used in the simulations have a grammage different from the grammage of the real paper, the amount of light scattering material in the simulations will differ from the real paper and the model will be less realistic. If we then find parameters producing results that agree well with measurements, it is likely that these parameters will not correspond to real physical processes. If the same parameters are used in another context they might produce totally different and erroneous results.
The method described above provides a way to determine grammage, but it is not very practical. The grammage is affected every time most input parameters to Grace are changed and every time we want to know the grammage we have to set up a simulation with all refractive indices to 1, put the light at an incident angle of 0◦, etc. This can of course be automated in the interface, but the simulations are time consuming.
To address this problem, a way to analytically determine grammage in real time was developed during this project. The mathematics are based on Markov theory [16].
Below is the obtained formula to calculate the grammage of a basesheet composed only of fibers and pores. Appendix C covers in detail how to calculate grammage in the general case with several components.
ω = ppfπrf 1 −1 − 1−ef 1+efwt (1 − wt) 1−e f 1+ef ppfrfπ 1−ef 1+ef + 8 3pf prp 1−ep 1+ep bρf (6.4)
where b is the average paper thickness, ρf is the density of the fiber wall, rf
is the average fiber radius, rp is the average pore radius, wt is the fiber wall
thickness, ef is the fiber ellipticity, ep is the pore ellipticity and pf p and ppf
are the average probabilities for fiber-to-pore and pore-to-fiber transition. If a contact reduction factor is used, these two probabilities can be calculated with the following formulas:
ppf = 1 − cppp 1 − pp pf (6.5) pf p= 1 − cfpf 1 − pf pp (6.6) 28
CHAPTER 6. IMPROVEMENTS where pf is the average fiber probability and ppis the average pore probability.
cp and cf are the contact reduction factors for pores and fibers respectively.
6.2
Fluorescence
As mentioned earlier, fluorescent filler materials are often added to paper to increase the whiteness and brightness. The yellow color that can sometimes be perceived on paper is due to absorption of blue light. Fluorescent materials can absorb invisible light in the ultra-violet spectrum and re-emit visible blue light, thereby compensating for the absorption.
During this project, fluorescence was only implemented in the basesheet model, but it can easily be extended to the radiative transfer model as well. An implementation of fluorescence was previously done for the Matlab version of Grace [10]. This version had the drawback of needing several runs, one for each wavelength and the transfer between wavelengths was done by defining internal lightsources for all absorbed wave packets. The current C++ version treats only one wave packet at a time and its wavelength can be updated during the simulation.
Figure 6.1 is an example of a measured fluorescence in a paper coating. The curve to the left is the absorption for different wavelengths and the curve to the right is the re-emitted wavelengths. The curves overlap in the middle. Some light with wavelengths up to 430 nm is absorbed and some light with wavelengths down to 390 nm is re-emitted. However, we need to make sure that no wave packets are transferred from longer to shorter wavelengths within the overlap region because photons of shorter wavelengths have higher energy.
In order to simulate fluorescence, an emission spectrum matrix must be provided, specifying the coupling between wavelengths, see figure 6.2 for an example. This matrix was obtained through measurements on a real fluorescent material [10]. The line in the bottom right marks where the absorbed and emitted light have the same wavelength, and it can be noted that there is no transfer below this and thus no increase in energy. The two curves in figure 6.1 are actually generated from the matrix in figure 6.2 through integration over the rows and columns respectively.
In Grace, parameters can be specified for different wavelengths and the para-meter corresponding to a certain wavelength is calculated through interpolation during simulation. One wave packet only represents one single wavelength and white light must therefore be represented by a multitude of wave packets of dif-ferent wavelengths. To increase efficiency, the difdif-ferent wavelengths that should be simulated are discrete and in the beginning of a simulation all interpolations of wavelength dependant parameters are precalculated and saved in vectors.
Fluorescence was implemented by letting a fraction of the absorbed wave packets be re-emitted at a different wavelength. The wavelength is chosen through the emission spectrum matrix (figure 6.2) where the column corre-sponding to the wavelength of the absorbed light acts as a probability density function.
CHAPTER 6. IMPROVEMENTS 300 350 400 450 500 550 600 0 1 2 3 4 5 6 7 8 9 10
Absorption and Emission Spectrum
Wavelength (nm) E ffi c ie n c y ( a .u .) Overlap Absorption Emission
Figure 6.1: Spectral absorption and emission of a fluorescent material. [10]
Figure 6.2: Example of an emission spectrum matrix of a fluorescent material. This kind of matrix can be produced by illuminating a material with a light with varying wavelength and measuring the re-emitted light. [10]
CHAPTER 6. IMPROVEMENTS
6.3
A New Pore Representation
In the previous version, pores were implemented as fillers with no scattering nor absorption and a refractive index of 1. This means that the concave pore surface was used for light scattering within the pore instead of the convex surface of the adjacent fiber or filler.
The implementation of the basesheet that is the result of this master’s project makes it possible to use an improved representation of pores. Instead of generat-ing random ellipsoids representgenerat-ing pores, the wave packet is instead propagated a random distance before encountering another component. The probability density function for this distance can be chosen to produce approximately the same distribution of light propagation lengths through pores as in the previous ellipsoid model.
If the size of a spherical pore is constant, the new distribution will be (see appendix D):
fh(h) =
h
2r2 (6.7)
where h is the traveled distance inside the pore and r is the pore’s radius. In order to keep a notion of pore ellipticity, this is then modified by the following (see appendix D):
h0= h v u u t (1−e 1+e) 2
1 + ((1−e1+e)2− 1) sin2(θ) (6.8)
where θ is the angle of incidence and e the ellipticity.
This novel pore representation will not give the same results as in the pre-vious one. Even though the pore sizes remain the same, the probability for scattering at the pore boundary is affected. A wave packet reaching the pore boundary will encounter a fiber at a random position. The distribution of in-cident angles and normals will be different and thereby also the probability for reflection.
6.4
A Different Approach to Surface Scattering
Until now, all scattering on the basesheet surface was done with the user-defined topography using general refractive indices defined for every layer. In some cases it can be of interest to let light be reflected directly on the paper components instead. This can for instance be useful to study how the fiber geometry affects the angular distribution of light. This has been studied in [9] with another Monte Carlo simulation model.
In order to perform similar simulations in Grace a feature was added in order to use the components’ surfaces and refractive indices when light scatters on the surface. When the wave packet deflects on the surface, a paper component is generated and the point of interception is calculated. The component’s normal at this point is then used for the scattering instead of the normal from the user-defined surface, see figure 6.3.
This new approach makes it possible to study a fiber network surface and how the geometry of the fibers affect the angular distribution of light. This
CHAPTER 6. IMPROVEMENTS
Surface Component Surface normal
Old version New version
Wave packet
Figure 6.3: New approach to surface scattering. The components’ geometries are used for light scattering at the surface instead of the user-defined topography. is particularly useful if the paper is coated because it is then very difficult to measure the basesheet’s surrounding surfaces.
Chapter 7
Results
This chapter presents the results of different tests. All tests were done on a model of eucalyptus paper sheets [3] and the results were compared to measure-ments and to previous simulation results.
First, a wave packet path visualization is demonstrated - a test especially useful for debugging purposes. Different tests aiming to verify the different improvements are then presented. The chapter concludes with a performance comparison between the Matlab version and the new C++ implementation.
7.1
Basesheet Simulation Example
A test was done where the wave packets’ paths through the fiber structure were logged and visualized. Figure 7.1 shows an example of the wave packet’s path through fibers and pores. This example is one of the very rare cases when a wave packet propagates a very long distance in the basesheet. Usually, the wave packet leaves the layer after only a couple of reflections/transmissions. Parameters are set to unrealistic values to only produce specular deflections. Multiple internal reflection in the fiber wall then produces the long spirals when the wave packet propagates around the lumen in the fiber wall. The fibers have a very small contact reduction factor because the light would otherwise propagate from fiber to fiber after only a few reflections.
The wave packet reaches the simulation boundaries twice, once at the bound-ary at y = 0 and once at x = 1000. The boundbound-ary conditions are then applied and the wave packet bounces back. A short fiber is clearly perceivable where the wave packet propagates back and forth between the fiber ends before it finally escapes.
7.2
Simulation of Real Basesheets Made of
Eu-calyptus Fibers
Grace was originally not conceived to take grammage into account. The fast analytical method for grammage determination described in the previous chap-ter is now implemented in the inchap-terface in Grace and this is very useful when
CHAPTER 7. RESULTS 500 600 700 800 900 50 100 150 200 250 300 350 400 -50 0 X Y Z
Figure 7.1: An example of a simulation showing one wave packet’s path through a fiber network. The red lines shows paths inside fibers and blue dashed lines are paths through pores. The parameters are unrealistic and set to promote light being trapped within fibers with no random scattering.
adjusting the grammage. Often many parameters have to be fine-tuned and the displayed grammage is now instantly updated for every change.
Using the analytical method to determine grammage, it is possible to verify simulations performed previously and try to adjust them to correct grammages. Knowing both the paper thickness and grammage also allows to better optimize unknown parameters.
Simulation models of eucalyptus sheets were set up previously [3] and the results were compared to measurements. The total reflectance of eucalyptus sheets of different grammages was measured before and after calendering.
The surface topographies were measured using a confocal laser scanning microsope. The pore sizes were estimated using mercury porosimetry and fitted against a log-normal distribution. Component distributions were acquired by studying cross-sectional images of the paper before and after the calendering with a scanning electron microscope.
The results obtained in [3] are shown in figure 7.2. The grammages of the simulated papers were too low. Therefore the total amount fiber material in the basesheet model must increase in order to better represent the real paper structure.
Models for a paper sheet of 120 g/m2 were also set up in the new
simula-tions. It was possible to increase the grammages to a correct mean level for the simulation models by modifying the fiber and pore ellipticities in the same way for all paper sheets. The grammages still did not agree perfectly with the measured values so it was necessary to fine-tune the pore sizes for individual paper sheets. With these new parameters the simulations were redone and the results are presented in figure 7.3.
The grammages of the simulated models are now in exact agreement with expected values. The simulated total reflectances are close to measurements but still do not agree perfectly.
The analytical determination of grammage makes it easy to make adjust-ments to the grammage but it is still difficult to optimize the parameters with 34
CHAPTER 7. RESULTS
E60 E60c E90 E90c 30 40 50 60 70 80 90 100
Grammages in previous simulations
Gr ammage (g/m 2) Paper Measured In simulation
E60 E60c E90 E90c 65 70 75 80 85 90 Paper Reflectance (%)
Reflectances in previous simulations Measured
In simulation
Figure 7.2: Grammages and reflectances for measured and simulated eucalyptus sheets before grammage adjustments. Two papers were studied, eucalyptus 60 g/m2 (E60) and eucalyptus 90 g/m2 (E90) before and after calendering
(calendered sheets are suffixed with ’c’).
E60 E60c E90 E90c E120 E120c 50 60 70 80 90 100 110 120 130 Gr ammage (g/m 2) Grammages Paper Measured In simulation
E60 E60c E90 E90c E120 E120c 65 70 75 80 85 90 95 Reflectance Reflectance (%) Paper Measured Simulated
CHAPTER 7. RESULTS respect to total reflectance.
7.3
Using Fiber Geometry for Surface
Reflec-tion
Jensen presented simulation results from paper sheets with aligned fibers using a free-standing simulation software ([14] and [9]). Similar simulations were carried out with Grace using the new feature that allows fibers’ geometry to be used for surface scattering. One paper with aligned fibers is lit from three different directions, and one paper with unaligned fiber is lit from one direction. The resulting Bi-directional Scattering Distribution Functions (BSDF) are shown in figure 7.4.
Jensen had a polygonal and static approach to fiber representation whereas Grace has a stochastic and elliptic representation, but the results show the same pattern.
The simulated BSDF of the paper with unaligned fibers is compared to BSDF measurements of a real paper sheet in figure 7.5, this time visualized in 3D. Even though the simulation model is not optimized to represent the real paper sheet, the two BSDFs have the same general appearance.
7.4
New Pore Representation
In order to verify the analytical distribution for pore sizes derived in appendix D and obtain distributions to the new pore representations for arbitrary pa-rameters used in the old version, simulations were carried out and the actual distances traveled through pores were logged.
Figure 7.6 shows the analytical probability density function (equation 6.7) and a histogram of logged distances actually traveled through the pores when the pores have a constant radius of 5 µm. The diagram shows that there is a good agreement with the analytical distribution.
In other cases when the probability density function for the pore radius is arbitrary it is difficult to analytically derive the new pdf.
However, pore sizes can be approximated by a log-normal distribution ac-cording to measurements [3]. Simulations show that, in this case, the distrib-ution of the actual propagated distances through pores and the distribdistrib-ution of the pore radii are equal. The results from a simulation carried out with spheri-cal pores are presented in figure 7.7. The distributions of pore radii and actual distances traveled through pores are almost identical except for noise, and the difference is small enough that it is within the limits of measurement errors.
This means that the same log-normal curves can be used with the new pore model. Care has to be taken when the allowed radius size is limited. In this case no radius over 20 µm was allowed. With the new representation, distances up to twice the maximum radius must be allowed.
To verify that the actual propagated distances through pores have the same distribution with the new pore representation, the actual propagation distances through pores were logged in the old and new versions and a comparison is shown in in figure 7.8. The two curves overlap almost perfectly.
CHAPTER 7. RESULTS
Light source direction
Figure 7.4: Simulated BSDF in reflection for different fiber orientations. Three images were simulated using aligned fibers and one with fibers with random directions. The light source is placed at an incident angle of 45◦with a direction according to the arrow in the figure.
Figure 7.5: On the left is the same simulated BSDF as in the lower right of the previous figure. This can be compared to measurements of a real paper, shown on the right (courtesy of Hjalmar Granberg [9]).
CHAPTER 7. RESULTS
0 2 4 6 8 10
Pore size (µm)
Probability density
Pore size histogram for constant pore radius
Simulated distances through pores Analytical PDF
Figure 7.6: Histogram of actual traveled distances through pores in simulations when the pores are spheres with constant radius 5.
0 5 10 15 20 25 30 35 40 45 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Pore size (µm) Probability density
Pore size histogram
Traveled distances through pores PDF pore radius
In other cases, when the probability density function for the pore radius is ***godtycklig*** it is difficult to analytically derive the new pdf.
Figure 7.7: The probability density function of the pore radius and wave packets’ actual distances traveled through pores agree well with each other.
CHAPTER 7. RESULTS 0 5 10 15 20 25 30 35 40 45 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Pore size (µm) Probability density
Comparison between old and new model
Traveled distances through pores with the new pore representation Traveled distances through pores with the old pore representation
Figure 7.8: The same probability density function for pore size is obtained with the new pore representation.
We compared the results of simulations with the old pore representation (section 7.2) to the results with the new representation. The results are shown in figure 7.9. As mentioned in section 6.3, we can not expect to obtain exactly the same results with the new pore representation. It seems like the calendered sheets differ more than the uncalendered and this is probably because the pores of the calendered sheets are highly elliptic, which produces higher incident angles at the pore boundaries with the new pore representation. As a result, the wave packets will tend to reflect more often and this will have an influence on the reflectance.
7.5
Fluorescence
For the testing of fluorescence simulations, the model of the 60 g/m2eucalyptus
sheet used in the simulations above was modified and a fluorescent filler mate-rial was added. A simulation with 100,000 wave packets per wavelength with wavelengths ranging from 300 nm to 700 nm was run. The fluorescence matrix in figure 6.2 was used. The results are shown in figure 7.10. We see that some of the light absorbed in the range 350-400 nm is transferred to longer wavelengths with a peak at around 425 nm, just as one would expect with this fluorescence matrix. In the case with no fluorescence, the light was simply absorbed.
A previous study of fluorescent coatings produced similar results with the same kind of curve for the spectral radiance factor [10]. This gives an indication that the fluorescence works well in Grace. However, no measurement data from a real fluorescent paper sample was available so it is difficult to draw any conclusions as to whether the results are realistic. Input parameters for the fluorescence process most certainly need to be adjusted.
