Development of an Open Source Platform for Monte Carlo Simulation of Paper Structures

(1)

IT 09 013

Examensarbete 30 hp April 2009

Development of an Open Source

Platform for Monte Carlo Simulation of Paper Structures

Jie Shang

Institutionen för informationsteknologi

(2)

(3)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Development of an Open Source Platform for Monte Carlo Simulation of Paper Structures

Jie Shang

Paper is a very complex optical material, and for optimization of the optical properties of paper, different advanced models of the interaction of light with complex structures have been studied. Among these models, the Monte Carlo model has been chosen and its routines have been developed for different purposes in previous projects. In order to make this type of simulations more accessible, and future contributions from the scientific communities easier to apply, the main goal of this project is to develop an open source simulation platform named OpenPaperOpt for light scattering simulation from fiber networks. As part of this project, this thesis aims to integrate the P3D model which was previously implemented and tested on the OpenPaperOpt. After the integration, the OpenPaperOpt platform has a general lights scattering process on the static fiber network. The thesis then gives guidelines for further development of different network models.

Tryckt av: Reprocentralen ITC IT 09 013

Examinator: Anders Jansson Ämnesgranskare: Olle Eriksson

Handledare: Ludovic Coppel / Per Edström

(4)

(5)

Acknowledgements

I would like to thank my thesis advisors, Ludovic Coppel and Per Edström, for their encouragement, knowledge, and guidance. I am grateful for their endless patience and understanding taking the time to work through any problems or questions that I en- countered throughout this project. I would also like to thank Olle Eriksson for being my thesis reviewer and for his support in my efforts. I thank him for his patience and time.

I am very grateful to the Kempere foundation which provide the scholarship to us, I also would like to thank STFI-Packforsk AB provided computer resource I could access.

Finally, I would like to thank my parents who supported me though my education.

Without them, this would not have been possible.

(6)

(7)

Introduction

1.1 Background

In order to optimize the optical properties of paper, advanced model for the interaction of light with complex paper structures has been studied. Monte Carlo simulation routines have been developed for different purposes in previous projects. These make it possible to continue development of the simulation models in an Open Source platform for faster development and broader dissemination. It will make this type of simulations more accessible, and future contributions from the scientific community will hopefully make the model even more applicable. This thesis is a part of a large five-year national project, PaperOpt, financed by Vinnova and eight industrial partners. It was performed in collaboration with two other Master’s thesis projects, one closely related, about statistical fiber networks, and one focusing on data format for exchange of inputs and outputs between different models and programs.

1.2 Project Objectives

This thesis is part of the project which is aiming to develop an open source simulation platform for light scattering simulation from different models of the fiber network. One previously implemented model (P3D) was fully implemented and tested. The focus of the thesis is to embed the P3D into the open source simulation platform and give the guidelines for further development of different fiber network models.

(10)

Chapter 1. Introduction

1.3 Previous Work

Two light scatterings in paper models have been developed, one is called GRACE [1], and the other is P3D [2]. In principle, both of them perform the similar light scattering procedure on the paper models, but the paper models they used are quite different.

1.3.1 GRACE

In GRACE, paper is described as a three dimensional structure including rough surfaces, coating layers, transparent layers, ink layers, and basesheet layers containing fibers, fillers, pores, and fines. The computations of light scattering from a paper are performed in three dimensions using a Monte-Carlo model. They follow fundamental laws of optics as well as some approximations; the input to the light scattering model includes a detailed description of the paper structure and knowledge of the scattering from real papers.

1.3.2 P3D

As a variant of GRACE, P3D was developed under the dogma that the structure of P3D should be usable in the light scattering model GRACE. So P3D have the same light scattering model as GRACE, but the generated paper in the P3D is quite different from the GRACE. The paper in P3D should initially consist only of fibers. Another difference is that the simulated paper in P3D was completely generated before the light scattering process. However, In GRACE, the paper structure is randomly generated when the light hit the paper structure.

1.4 Thesis Overview

Chapter 2 introduces some basic properties of the real paper. In chapter 3 some different models used in the paper optics are presented and the model applied in the OpenPa- perOpt is discussed. Chapter 4 presents two models of the fiber network, one is the fiber network currently implemented and integrated into the OpenPaperOpt, and another is the model which could be implemented in the future. Chapter 5 addresses the general light scattering procedure in the OpenPaperOpt and gives the specification on how to reuse the light scattering procedure. Chapter 6 gives implementation related to embedding works. Chapter 7 presents the results obtained from the OpenPaperOpt and

(11)

Chapter 1. Introduction

suggests some future improvements. The left appendices show the input specifications to generate the results referred in the chapter 7.

(12)

Chapter 2

About Paper

Although paper is a thin material, it is a quite complex object. In order to understand the paper model used in this thesis, it’s better to have some knowledge about the structure and components of the real paper. This chapter will briefly illustrate the components of the paper.

2.1 Paper

Paper is mainly made of fibers. The fibers can be acquired from the wood, rags, and grasses. The process which extracted the fibers from wood is called pulping. Pulping can be done in several methods, which range from purely mechanical, in which the wood is ground into fibers by disk refiners or grindstones, to chemical, in which the fibers are separated by chemically degrading and dissolving the lignin that binds them together in the tree. After the initial step, the pulp is formed. In modern paper industries, pulps contain not only fibers but also some other fillers, such as chalk or china clay. These fillers are mainly used to improve the characteristics of the paper for printing. Then, the pulp will be shaped to a sheet in the paper mill and follow the several operations:

pressing, drying, calendaring, sizing and coating. These operations are mainly to remove the water from the sheet and make the paper be smooth. The purpose of sizing is to curb fibers’ tendency to absorb liquids. Coating different surfaces on paper could impart certain qualities to the paper, including weight and surface gloss, smoothness.

(13)

Chapter 2. About paper

2.2 Natural Fiber

The most common ingredient in the paper is fiber, so we should know much about the fiber. The wood fibers are basically long rectangular like pipes with tapered edges. The fibers are glued together in the wood by the middle lamella. Each fiber is divided into several layers as shown in the figure 2.1. The cavity in the middle of the fiber is called lumen. Due to the several mechanical processes used in the paper making, the thin fiber walls tend to collapse and form different shapes of the lumen.

Figure 2.1: The most upfront cell has been sectioned in order to show the layered structure characteristic of each cell, comprising the primary wall, the three secondary wall layers (S1, S2, S3) and the lumen. The grey lines in the secondary walls layers represent idealized cellulose microfibrils. The cellulose microfibers in the middle secondary layer appear twisted with respect to the cell’s vertical axis. The angle they form, called the microfibril angle (MFA), plays a crucial role in determining the stiffness of

the wood.

The thicknesses of the fiber walls vary a lot with seasons. In spring, fibers with thin walls and large cavities are formed to manage the large need for water transportation, in the late summer; fibers with thicker walls are created. This can be viewed as annual rings. The thick walled fibers grown in the late autumn behave quite differently from the thin fibers grown in the spring. The angle of the building blocks in the S2 layer is an important factor for the fiber stiffness. A small angle will generate a stiff fiber and a large angle will generate a more elastic fiber. The angle is mostly a random variation but is constant within one fiber [2].

(14)

Chapter 2. About paper

The length of the wood fibers varies from about 0.5mm to 5mm. The variation in fiber length comes both from natural variation and from the types of wood.

The coefficient of refraction for cellulose is, according to [2], about 1.55.

(15)

Chapter 3

Light Scattering Model

This chapter briefly presents the light scattering model used in the OpenPaperOpt.

3.1 Overview

Although the paper is a very thin object, when we saw it under the microscope, we will find paper actually is a very complex material. In order to study the optical properties of the paper, we should make some simplifications for paper. As the time goes by, different models are proposed for this purpose. In the 1930s, the Kubelka-Munk (KM) theory [3] was developed. 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 [4]. In the 1971, the Discrete Ordinate Radiative Transfer Theory (DORT) [5] was developed and this model could be considered as the extension of the KM-model. Both of these two models lack physical reality and have some defects in the light scattering procedure. Besides them, there is another model named as a city Monte-Carlo. The Monte-Carlo model used in the light scattering model is a more physical model than KM or DORT. It uses experimental data describing the actual paper sample structure and the properties of its components, to simulate the light propagation inside the paper [6]. The Monte-Carlo model is applied to this project and will be further discussed in the following sections.

3.2 Monte-Carlo Model

The city of Monte-Carlo is well known for its casinos and also has been a name of computation method which imitates the random behavior of the games of chances in

(16)

Chapter 3. Light Scattering Model

a casino. The variables in the Monte-Carlo model are randomly generated but have a known range of values with a probability distribution. Using the Monte-Carlo model in light scattering, we can use the model in the same way to approximate the properties of a multiple-scattering medium. This could be achieved by making assumptions on the random behavior of light. The light is simulated as randomly generated wave packets.

It is simulated how they propagate through and interact with the paper. One of the advantages of using the Monte-Carlo model is that there is no limit to the complexity of the model of digital paper. For example, if we want to know how much of the incoming light is reflected back from a paper sheet, the KM-model would be to describe the paper sheet by a mathematical model and analytically calculated the ratio reflection/trans- mission. In order to achieve this, the model of the digital paper needs to be simple enough to solve the equations exactly. When using a Monte-Carlo model, there is no such constraint on the complexity of the model of the digital paper, it is a matter of having enough computer power to simulate enough wave packets to sufficiently minimize the statistical errors [4].

3.3 Paper Structure

In order to simplify the problem of the light scattering on the paper, we could consider paper as a layered material. Thus, we could focus on the light scattering of different layers individually. The layers have different types, for example, coating layer, ink layer and basesheet. Between different layers there is a surface to separate them.

Figure 3.1: Paper structure in light scattering model, the paper contains different layers and surface separate different layers [6].

(17)

Chapter 3. Light Scattering Model 3.3.1 Statistical Structure Layer

Statistical structure layer can be used to study how different components affect the optical properties of the paper. When using the statistical structure layer to model the basesheet, it could contain three types of components: fiber, filler and pore. Fibers are modeled to be hollow cylinders with ellipsoidal cross-sections; pores are modeled as ellipsoids; and fillers are modeled in the same way as pores but made of a material that can scatter and absorb light. The position of each type of components is not static but determined by the distribution functions. As the light scattering goes in depth in the layer, the probability of the occurrence of different components are calculated by the distribution functions.

3.3.2 Static Description

Compared to the layer mentioned above, the positions and geometrical properties of the components in this kind of layer are first generated and then light scattering is simulated in this structure. Section 4 will present the details about how to construct the static fiber network in the application.

3.3.3 Surfaces

A surface is an interface separating two layers in the paper structure. In the current light scattering model, Surface itself doesn’t have the thickness but it is used to define the thickness variation of layers. The surface topography is defined by a height map;

these are matrices of height values. The sampling points are often evenly spaced but one could also consider using uneven spacing to achieve higher spatial resolution in parts of the surface while having a lower resolution in other parts to save space and computation time [7].

3.4 Light Source

The light scattering simulation can’t be processed without light source. The light source emit light ray to illuminate the simulation volume. We can think of a light ray as the straight path followed by a light particle (called as a photon) as it travels through space.

In the physical world, a photon, carries energy, and then enters the different layers of the digital paper, energy will decrease gradually. In our simulation model, we model a photon as a wave packet. Light source could initialize the wave packets and send

(18)

them to the paper volume. currently, Three types of light source are implemented in the program.

• Beam Source acts mostly as an incident laser beam. With different initial distribution parameters, it can generate three different types of light distributions on the paper volume.

1. pencil all wave packets are incident at one point.

2. gaussian the wave packets are distributed over a circle on the paper sample according to the Gaussian distribution.

3. uniform the wave packets are uniformly distributed over the paper sample.

• Lambert Source is used to generate the diffuse the light ray to the paper sample.

It has the same initial parameters as the Beam Source.

• Elrepho Source is the paper industry’s own spectrophotometer. It measures colour, brightness, opacity and whiteness of paper, paperboard, pulp, coating inks and fillers. It will generate wave packets according to the specification of an Elrepho light source.

3.5 Detector

Detector is used to collect the wave packets which have left the paper volume and generate the result according to the status of wave packets. Detector implemented in the simulation could gather the wave packet whether it is transmitted or reflected from the paper volume. Moreover, the program supports several detectors working simultaneously. In order to fully analysis the properties of the paper, several types of detector are implemented. Some of the typical detectors are listed below.

• ARS (Angle Resolved Scatterometer) detector measures the reflectance or trans- mission at different aperture. The power reflected or transmitted is transformed into a Bidirectional Scatter Distribution Function (BSDF ) [8].

• Image detector works like a simple camera. The wave packets reflected from or transmitted through the paper structure have to pass through an aperture and are then counted in an image plane forming a two-dimensional image of the light coming from the paper sample.

(19)

• Elrepho detector returns the reflectance factor which is compatible with an El- repho instrument. The power reflected within the Elrepho detector aperture is normalized by the power that could be measured from a diffuse sample.

Figure 3.2: Image generated by different type of detectors with different types of the light sources. Top left is the combination of lambert light source and the ARS detector, Top right is the combination of elrepho light source and the ARS detector, Bottom right is the combination of the beam light source and image detector, and the bottom

left is the comparison of different light sources.

3.6 Light Scattering

Light scattering is one of the physical process to analysis the optical properties of the materials. In this project, the scattering is conducted for the paper in a Monte-Carlo ap- proach. The light generated from the light sources hits the paper object and is collected by the detectors. In order to model the path that light travels to the paper object, the scattering procedure in the project is to simulate the travel of a large number of wave packets propagated by light source. During the simulation of one wave packet, the status of wave packet is changed by what it hit in the paper. The wave packet may hit the different surface, for example, the coating layer, the ink layer and the basesheet.

In the different layer, there are different probabilities for the scattering and absorption occurring. The detector will collect the wave packets reflected from the paper object. It could generate different result image by different types of the detectors.

(20)

Chapter 3. Light Scattering Model 3.6.1 Wave Packet Description

The wave packets in the simulation act alone and don’t interact with others. Each wave packet contains information about its current state. This information is updated as the wave packet propagates through the paper structure. The status of a wave packet includes:

• Position vector

• Direction vector

• Polarization vector

• 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 large number of wave packets with random polarization directions the effect of unpolarized can still be achieved.

3.6.2 Scattering at Surfaces

Light scattering from surfaces is a large area of study in optics, radar, and antenna research. If we suppose that the surface is perfectly flat, the angle of the reflected and the transmitted light can be determined with the well-known law of reflection 3.1 and Snell’s law 3.2, Both are shown in the figure 3.3, where n_i and n_t are the indices of refraction for the two media at each side of the interface, and the θ_i and θ_t are the angles from the surface normal at which the incident and transmitted waves propagate.

According to the equation 3.1, the ray reflected by the surface will propagate at the same angle from the surface normal as the incident wave.

θi = θr (3.1)

nisinθi = ntsinθt (3.2)

The amount of reflected and refracted light are determined by the Fresnel relations 3.3

(21)

Figure 3.3: An illustration of Law of reflection and refraction (Snell’s law)

R_s =

n_icosθ_i− n_tcosθ_t n_icosθ_i+ n_tcosθ_t

2

(3.3)

Rp =

ntcosθi− n_icosθt

n_tcosθ_i+ n_icosθ_t

2

(3.4)

Where R_s is the reflection coefficient for light with a polarization perpendicular to the plane of incidence and R_p is the reflection coefficient for light with a polarization parallel to the plane of incidence [4].

The light scattering on a surface will be conducted in the following steps. First the intersection of every wave packet with the surface waviness will be calculated. At each intersection, the normal of the surface will be calculated and the Fresnel equations are applied together with Snell’s law to determine the direction, magnitude, and polarization status of the refracted and the reflected light ray. The scattering by the micro-roughness is introduced by letting a fraction of the reflected and refracted light ray distributed in a Lambertian manner. The ratio between the scattered and totally reflected or refracted light on a surface with a Gaussian surface structure could be calculated by the following equations [9].

diffusely reflected power

totally reflected power = 1 − e⁻⁽4πRq nicosθi

λ )²

(3.5)

(22)

diffusely refracted power

totally reflected power = 1 − e⁻⁽2πRq (nicosθi−nicosθt)

λ )²

(3.6)

In the equation3.5and 3.6, the n_iand n_tare the indices of refraction for the two media delimited by the surface. The θ_i and θ_tare the angles from the surface normal where the incident and transmitted light ray generated, The wave length of incident light ray are given as λ in the above equation and the Rq is the roughness of the surface where the reflection and refraction occurred. The polarization of the reflected or refracted light is determined by the Fresnel formula3.3,3.4.

For the scattered transmitted light, the polarization factor cannot be described by an- alytical functions, for the lack of better theory; the polarization of the scattered transmitted light could be determined by the Fresnel equations then rotated along with the direction vector of the light after a scattering event occurs.

3.6.3 Scattering at Static Fiber Network

Scattering at the static fiber network likes surface described above. It also needs to find the intersections, Due to the large amount of fibers exist in the fiber network, It’s important to have some methods to accelerate finding the intersections. So the light scattering on the static fiber network will be conducted in fast ray tracing fashion. Two main strategies are applied to create a suitable data structure for this purpose.

• Bounding volume is the most fundamental and ubiquitous tool for fast ray trac- ing acceleration. This volume contains a given object and permits a simpler ray intersection check than the object itself [10].

• Spatial Subdivision is processed top-down, partitioning a volume bounding the environment into smaller pieces. The smaller volumes thus formed are assigned collections of objects which are totally or partially contained with them [10].

By using the above strategies, the first intersection point could be found quickly. During the simulation the program will access the fiber network for the closest intersection of a ray and any surface in the paper volume. A surface could be a surface between air and fiber or between fiber and lumen. At the very beginning of the simulation, the ray will begin at a known intersection, which could also accelerate the simulation.

(23)

Chapter 4

Fiber Network Generation

In this thesis, my work is to implement the fiber network model from P3D and fit it into the OpenPaperOpt program with newly developed XML input module [7]. So this chapter is dealing with the static fiber network. It firstly introduces the fiber network model used in the P3D, and then discusses another fiber network model which could be used in the later in the OpenPaperOpt. The fiber network generated here will not contain fillers or fines. It will only contain the fibers. The position of the fibers in the static fiber network is predefined. In the process of making a geometrical mode of fiber material, we first have to consider the geometry of a single fiber. Then with help of certain stochastic algorithm, we distribute the fibers into the paper volume to form our desire fiber network.

4.1 Staticsheet

This model is used and implemented in the current project. The model Staticsheet consists of three procedures to generate the final static fiber network. The first step is to model the basic unit in the fiber network, the basic unit in this model is the fiber segment, and the second step is to use the fiber segments to generate a desire fiber by using certain algorithms, the last step is to determine the position of each fiber in the paper volume to form the 3D paper structure. The algorithms described in this section referred to [2], if no other resources are specified.

4.1.1 Fiber Segment Model

In the real world, the wood fiber will appear in various shape and curvature. In order to model the both the sharp kinks of fibers originating from the wood and the smooth

(24)

Chapter 4. Fiber Network Generation

Table 4.1: The properties of the fiber segment

P arameter Description remarks

length length of the fiber segment relative to the length of the fiber radius radius of the fiber segment relative to the radius of the fiber

angle orientation relative to the previous segment

position the coordinates of the fiber segment relative to the previous segment

flowing of an industrially manufactured polymer fiber, we should model the fiber as general as possible to make it useful in a wide variety of cases. Based on the above discussion, from three dimension view, we could model every single fiber as a hollow cylinder. If we view the fiber in two dimensions, it is just a single straight line. However, a single fiber in real life should have a varying curvature and also out of plane curvature.

So we could link several cylinders with different orientation together to approximate the shape of a natural fiber.

Thus, the basic unit of a fiber should be called as fiber segment, which is defined as a single straight cylinder. See figure4.1

Figure 4.1: An illustration of fiber segment shows the properties of the fiber segment, fiber length and radius

To define a fiber segment, we need several parameters to describe the necessary properties of the fiber segment. See table 4.1.

In the real world, the above parameters are only a part of the properties of a fiber segment. For the simplicity of the mathematics and programming needed, when we implement the fiber segment, actually we only consider two parameters listed in the table 4.1. They are the angle and the position of the fiber segment. The length of fiber segment could be considered as a constant value. Another reason is that the length of a whole fiber is also a constant value. By using these two constant values, we could calculate how many segments a fiber contains. And sometime, the length of the last fiber segment will be shorter than other segments contained in the same fiber. The radius of a fiber segment should be the same as the radius of the fiber. So this parameter could derive from the fiber which contains the current fiber segment.

(25)

Table 4.2: The properties of the fiber

fiber length length of the fiber independent

radius radius of the fiber independent

form factor to calculate the curvature independent

angle orientation angle generated by distribution function segment length length of the fiber segment independent

number of segments the number of segment relative length of the fiber and fiber segment angles for segments orientation angle for each segment relative to the form factor and previous segment

segments fiber segments relative to each individual fiber segment

4.1.2 Fiber Generation

The complete fiber consists of one or more fiber segments. They will be joined together so that they are continuous in tangent in the connecting joint or kinked, i.e. discontinuous in tangent. The table 4.2shows the parameters to generate a fiber.

The fiber length and segment length determine the number of fiber segments contained in one fiber. The length, radius and form factor of each fiber could be generated by different distribution functions. Data from a machine called Fibermaster [11] is information about several individually measured fibers. For each measured fiber there are length, width and form factor measurements. This makes it possible to use these data to construct a distribution function, i.e. we could use data from Fibermaster [11] to define the length, radius and form factor for a fiber.

The most important characteristic of an individual fiber are the curl and angles for numerous fiber segments. The curl of the fiber is dominated by the angle of the fiber and the angles of the segments. Here, Curl is defined by the relationship between the length of the fiber and its maximum extent which is the diameter of the smallest sphere that the fiber could fit into, as shown in the figure4.2

The following equations are used to calculate the curl and other factors. The variables used in the equations are shown in the figure 4.2

curl = l

D − l (4.1)

curl factor = l

D (4.2)

(26)

Figure 4.2: Definition of measures used for different definition of fiber curl. The marked lengths are used in the equations4.1to4.3. l is the fiber length, d is the length of the diagonal in a bounding box and D is the largest extent measurable with a vernier

caliper.

form factor = d

l (4.3)

The angle θ for the fiber will be generated by the von Mises distribution. This distri- bution has the distribution function

f (θ) = e^{k cos (2θ)}

πI0(k) (4.4)

where I_n(k)is the modified Bessel function of order n. The parameter k could be estimated from the elastic modulus of the paper by expanding the equation in the following form.

f (θ) = 1 π

∞

X

i=0

aicos(2iθ) (4.5)

If the equation4.4 is expanded by using the following term

a_i = 2I_i(k)

I (k) (4.6)

(27)

Neglecting terms (a) of order greater than two greater with the relation

E_{M D}

E_CD = 6 + 4a₁+ a₂

6 − 4a₁+ a₂ (4.7)

The parameter k could be calculated out. In the equation 4.7, E_{M D} and E_CD define the elastic moduli in the machine and cross-machine directions.

4.1.3 Fiber Network Generation

When the fiber is complete it must be prepared for insertion into the network. The model described here will use the P AKKA [12] approximation to estimate the bending of a fiber, however the P AKKA [12] only provides the depth of the each fiber segment.

The placing of the fibers in the plane (two dimensions) has to be generated by other methods. Thus, the whole process of generating the fiber network will contain two stages. First, the fiber is placed in the paper plane (two dimensions). Secondly, the depth of the fiber is estimated by using the P AKKA [12] estimation.

4.1.3.1 Step One

Several factors, length, radius and angle, could involve deciding the placement of the fiber in the plane of the paper. The fibers will be placed in several random positions and a measure of goodness is calculated for each placing. The placement with the best measure will be selected.

The measure of the goodness, G, used the following the equation

G(Af, ω(¯p), ωmean) = Z Z

p¯

∈ A_fPg(ω(¯p) − ωmean)d¯p (4.8)

The A_f in the equation 4.8 is the area covered by the fibers, ω(¯p) is weight per square meter at the point ¯p and the ωmean is the mean value of the weight. P_g is a polynome of degree 3 satisfying the conditions:

(28)

Pg(0) = 0 (4.9)

P_g(ω_mean) = 1 (4.10)

P_g⁰(ω_mean) = 0 (4.11)

P_g⁰⁰(ω_mean) = 0 (4.12)

A small value of G is considered to be good.

4.1.3.2 Step Two

The step one deals with placing fibers in two-dimension. If we put the fibers in Cartesian coordinate system, after the step one, the fiber will get the position in XY-plane. Step two is to give the Z-value to each fiber. That is to determine the depth of each fiber in the paper volume. The principle of this step is to place each segment as low as possible.

It’s only constraint by the previously placed material, the P AKKA [12] approximation and the rest of the fibers. The fiber are placed one at a time and one segment at a time starting in one end of the fiber being placed.

Figure 4.3: Schematic image of the restrictions on the placing of a fiber segment,zi.

The figure4.3illustrates the restrictions when placing one segment, z_i. The first pair of constraints, a and b are calculated using the P AKKA [12] approximation which limits the elevation of the fiber segment z_i to:

(29)

|z_i− z_i−1| ≤ t_fTf (4.13)

In the equation4.13, the t_f is the fiber thickness and T_f is calculated by

Tf = (C · ω_f · W F F )^1/4 (4.14)

Where C is a constant value given in Nm, ω_f is the segment length in m, and W F F is the wet fiber flexibility in 1/N m². c in the figure4.3is the next constraint which comes from the local geometry and exists between a and b. The last constraints are the d1 and d2 which come from the segments that have already been placed in the paper volume.

4.1.4 Boundary Conditions

During the procedure of generating the digital paper, the boundary condition used in static fiber network is periodic in the plane of the paper. This means that the structure can be tiled with itself to form an infinitely large sheet of paper. When a fiber passes a boundary it is translated to the other side of the structure. Long fibers may however cross themselves after the application of the boundary conditions. Fiber segments normally also reside on both sides of the boundary being crossed. See figure4.4

Figure 4.4: Illustration of the boundary conditions used in generating the fiber network. Fiber 1, 2 and 3 show example of the different situations in the network, Fiber 3 is an example of self-intersection mentioned in the text. The square represents the

boundary of the network.

(30)

4.2 Other Network Generation Model

This section introduces another model which could be used later in the OpenPaperOpt program. This model has been fully studied during the thesis work, but it is not used and implemented in the current program. It’s not difficult to implement this model based on the model Staticsheet. Because both of the models have the same procedure to generate the fiber network. Compared with the model Staticsheet, the fibers in this model have more curvature . This is because the basic unit in this model is different from the model Staticsheet. The basic unit in this model has more features which make it easy to generate the fibers with realistic curvature. The model presented in this section referred from [13].

4.2.1 Fiber Segment Model

In this model, the fiber segment is defined as a single circle arc which is shown in the figure4.5

Figure 4.5: An illustration of fiber segment showing the parameters defining fiber segment

In the figure4.5, c is the central point of circle arc, r is the radius of the circle, u and v are the base vectors and the angles between the u vector and the start and the end point of the arc, α₁ and α₂, respectively. The points in the fiber segment could be calculated by the following equation 4.15, and the variable α is defined in the scope [α1, α2] we could use these points to visualize the fiber segment.

x = c + r · cosα · u + r · sinα · v (4.15)

(31)

Table 4.3: The properties of the fiber segment

l length of the fiber segment independent r radius of the fiber segment independent

θ orientation angle relative to the previous segment θ_k1 kink angle 1 relative to the previous segment θk2 kink angle 2 relative to the previous segment

t_ec,prev unit vector belong to the previous segment

n_prev plane normal belong to the previous segment de,prev ending tangent vector belong to the previous segment p_e,prev segment end point belong to the previous segment

To define a fiber segment, we need several parameters to describe the necessary properties of the fiber segment. See table4.3

In the table4.3, θ_k1 is defined as the angle between the project of d_s(see figure4.6) onto the previous segment plane and the vector d_e,prev of the previous segment while θ_k2 is the angle between d_s and the previous segment plane. θ, θ_k1 and θ_k2 determined the position of the fiber segment relative to the previous segment.

In the figure 4.6, It illustrates the parameters t_ec,prev, d_e,prev and p_e,prev. n_prev defines the segment plane normal.

Figure 4.6: Illustration of parameters for generating fiber segment, tsc is the unit vector pointing from the start point to the center, dsis starting tangent vector and ps

is the segment starting point.

The first five parameters in the table 4.3 are for every fiber and the rest parameters are derived from the previous fiber segment. If the segment is the first of the fiber, no previous exists and the values of the last four parameters in the table 4.3 are given arbitrarily.

(32)

Table 4.4: The properties of the fiber

l length of the fiber independent nseg the number of segments independent

4.2.2 Fiber Generation

The fiber generation process in this model is slightly different from the model Staticsheet.

In model Staticsheet, we firstly have the fibers and then generate the fiber segments on the specific fiber with the desire parameters. In this model, the procedure is reversed;

we firstly have the fiber segment and then generating an integral fiber by rotating fiber segment around its neighbors. By this method, the curl of a fiber depends on segment- to-segment orientation θ, θ_k1 and θ_k2, segment opening angle α, number of the segments and ratio between the segment radii. Due to the complexity of curl dependencies when more than one fiber segment, curl factor cannot be used as an input parameter. Beyond the curl factor, another important parameter of the fiber is the fiber length. The fiber length in this model is easy to control, with segment length as an input data for the fiber segment, the total fiber length is simply split between its segments. Based on the above discussion, if you want to generate the fiber with the designated curl factor, you should choose model one. So the parameters for generating the fiber are shown in the table 4.4.

To generate a fiber, the process involves heavy work on the generating the fiber segment through rotating vector in space around a specific axis. To rotate a vector (x, y, z) φ radians around an axis defined by a unit vector (n_x, n_y, n_z) the following equation is applied.

v⁰ = R(a, φ) · v (4.16)

where R and A are given by

R(a, φ) = I − sinφ · A + (1 − cosφ) · A² (4.17)

(33)

A =

0 −n_z −n_y

−n_z 0 n_x

−n_y −n_x 0

(4.18)

And I is identity matrix given by

I =

1 0 0 0 1 0 0 0 1

(4.19)

By applying the above equations, it’s easy to calculate the parameters defined in the figure 4.5 from the given parameters listed in the table 4.3, The further details about the calculation are fully presented in the [13].

4.2.3 Fiber Network Generation

In order to control the placement of the fiber, a center point between the end points is defined. When placing the fiber in the paper volume, the center point is positioned randomly into the paper volume without consideration that the fibers have already been placed in the paper volume. This implies that there may be several fibers that could take up the same space. So this model is suitable to construct low-density three-dimensional fiber network.

4.2.4 Boundary Conditions

The network is periodic. If the fibers are out of the boundary, it will be trimmed to fit into the network and the trimmings are also replaced. The boundary conditions guarantee that the original fiber will be complete and unbroken when placing in the paper volume.

The trimming procedure is done by segment-wise and works by redefining the opening angle α1, α2 and the center point of the circle arc. The figure 4.7 show an example of trimming process on a fiber segment.

(34)

Figure 4.7: Illustration of procedure of trimming a fiber segment, the original fiber will be replaced by five new segments identical to the original fiber segment. The

number indicates the identical segments.

(35)

Chapter 5

Light Scattering Process

OpenPaperOpt as an open source simulation platform for light scattering simulation on the fiber network, it provides a general light scattering process for different fiber networks. This makes the researchers keep focusing on generating different structures of the fiber network and study the optical properties of the different papers. This chapter shall introduce the data model OpenPaperOpt used to perform the light scattering, present the data specification and methods in order to use this data model. This chapter assumes that readers understand the basic rationale of the ray tracing.

5.1 Overview

When the static fiber network has been generated, it shall be prepared for light scattering simulation. OpenPaperOpt platform doesn’t perform light scattering directly on the raw data which represented the generated fiber network. First, it will convert the raw data to a more suitable format. This kind of data format can guarantee that light scattering can be performed on different types of fiber network. The second step is to create the data model for the fast ray tracing. The last step is to perform the ray tracing. After the first step, the following steps completely performed on the formatted data.

5.2 Data Model

The basic element in fiber network generated by the generator is segment. The data model which OpenPaperOpt used is built up based on the segment. Each segment in the final structure is represented by an instance of base element. When all the segments have been converted to the instances of base element, the fibers constructed from these

(36)

Chapter 5. Light Scattering Process

segments have also been converted and then the fiber network has been prepared for the light scattering process.

Before the light scattering process, the OpenPaperOpt will convert the raw data according to the data model. It started from converting the fiber segment. Fiber segments themselves perform the most of calculation.

Figure 5.1: A shows the data model used for converting the fiber segment. B illustrates the position of converted fiber segment related to the original fiber segment.

Numbers represent the faces; a base element has 16 triangle faces.

The base element used is shown in the figure 5.1. As shown in the image A, the cross section of fiber segment is abstracted to an octagon, and the inner octagon represents the fiber lumen. The size of the lumen is decided by the shrink factor. The default shrink factor used in the OpenPaperOpt is 0.4 if the user doesn’t use a specific value. The scope of the shrink factor is between 0.0 and 1.0. In the image B, Firstly it shows that the positions of the base element are calculated inside of the original fiber segment. The benefit of using this kind of positioning is that it makes that base element can be used widely no matter what the shapes used to model the fiber segment, i.e. In the current fiber network, the fiber segment is model as hallow cylinder, even it is modeled as cube, we could also use the base element to convert the raw data. Secondly, the surfaces of the lumen and segment are flat triangles. Flat surface are used to make detecting the intersection fast. Triangles used instead of the polygons are due to that if the polygon with more than three points is used, all the corners of the polygon can’t be in the same plane in various circumstances. This kind of base element converts the raw data as abstract as possible and makes the light scattering process in the OpenPaperOpt used in different fiber network as much as possible.

In order to use the data model defined in the OpenPaperOpt, the fiber network generator

(37)

side of the fiber segment. l_p and r_p shown in the figure 5.2are needed. They are easily calculated out by applying the some linear algebra formula.

Figure 5.2: Calculation of the two points used for the data model defined in the OpenPaperOpt

In the figure5.2, two circles denote the two cross sections of the fiber segment. Position vector c0 and c1 are the center of the circles. Suppose we know the radius of the fiber segment and represent it as variable r.

The unit normal vector n can be calculated by rotating the vectorc0^→c1 90 degree around c0.

c→0c1 = c1− c₀ (5.1)

n = Rotate(c0, 90)·c0^→c1 (5.2)

After the n has been calculated, it should be scale the by the radius r and the following equations applied to calculate the l_p and r_p.

lp = c0+ n · r (5.3)

rp = c0− n · r (5.4)

We now have all the desired parameters that the OpenPaperOpt required.

(38)

5.3 Data Preparation for Light Scattering

The light scattering in OpenPaperOpt is performed in a fast ray tracing fashion. When a wave packet is traveling in the paper structure outside any fibers, OpenPaperOpt could find all the fiber elements close to the ray. This procedure since should take long time to detect a ray against all the fibers in the fiber network. A two-level access structure is created for this purpose to accelerate the detecting.

First a number of segments in a fiber are grouped together and stored in a bounding box. The capacity of bounding box is controlled by the box length. The default box length is set to 5, i.e. each bounding box could contain five fiber segments if the users don’t give other specific value to the box length. After that, each bounding box will be added to a lattice box which is formed by dividing the whole paper volume into several sub volumes. The figure5.3 shows an example.

Figure 5.3: An example showing a paper volume divided into eight sub volumes, the size of each sub volume is decided by side length.

The OpenPaperOpt decided on how many sub volumes should be divided according to the original size of the fiber network and the side length which is shown in the figure 5.3. The default side length used is 30mm if user don’t give a specific value. And the generator of the fiber network should define an interface which could provide the size of the fiber network to the OpenPaperOpt. The size should be in 3-dimension, the size in x-axis, the size in the y-axis and the size in z-axis.

Until now, the light scattering process in the OpenPaperOpt can be started with the

(39)

5.4 Optimization

The data model used in OpenPaperOpt can be adopted for different types of fiber segments. No matter what types of the fiber segment are, the OpenPaperOpt converted it to a straight segment. This would affect the realistic of the result generated by the light scattering process if the model of fiber segment is a circle arc.

Figure 5.4: Illustration of using sub straight segment to displace the arc fiber segment for light scattering, numbers represent the straight sub segments.

In the figure5.4, the fiber segment is a circle arc. In order to increase the realistic of the result generated by the light scattering, we can first convert each arc to several straight sub segments before we forward the raw data to the OpenPaperOpt. Another reason we could do the conversion is that OpenPaperOpt only require the cross section of the fiber segment to format the raw data. The final shape which is consisted by the sub segments should be similar to the original fiber segment. Thus, after the processing, the OpenPaperOpt could give a more realistic result. When we face to different types of the fiber segments, the abstract data model used in OpenPaperOpt makes various optimization methods come out reasonably and naturally.

(40)

Chapter 6

Implementation

This chapter is concerned with integrating P3D which implemented the Staticsheet model stated in the chapter 4 into the OpenPaperOpt platform. Most of the classes from the P3D are reusable when integrating to the OpenPaperOpt platform. In order to keep the consistency of the program structure, the structure of some classes from P3D has been rearranged. Some new methods have been developed for providing the new features and cooperating the IO [7] of the OpenPaperOpt. As an open source platform, lots of developer would join in this community to add or modify the program; a well documented source code could provide the better readability to facilitate developers familiar with the program quickly. So new format documentation comment and documentation generator tool have been applied to program.

6.1 Documentation Code

Java documentation comment has been chosen as the format of the comments used in the OpenPaperOpt. Even the Java document comment format is designed for the Java programming language and the OpenPaperOpt is implemented in the C++ programming language, the well defined and structured Java documentation comment could be applied to the OpenPaperOpt perfectly. There are several advantages of the Java document comment.

1. It has define different tags for class, functions and parameter and return value of the functions which increase the readability of the source code.

2. It is widely used and parsed by different document generator tools.

(41)

Chapter 6. Implementation

3. It is easy to learn, new comers could learn it quickly and write the nice formatted the comments for the code fragments.

OpenPaperOpt is an open source platform, that means lots of developers would con- tribute to it, it should have a mechanism which could provide the development document quickly after the changes have been made to the source code. The code documentation generator tool, Doxygen which is a documentation system for C++, C, Java, and so on, have also been used. Doxygen can generate the on-line documentation (in HTML) and/or off-line reference manual in several file formats from a set of documented source files. It could also generate dependency graphs, inheritance diagrams, and collaboration diagrams, which could give more help to the developers. Another reason of choosing the Doxygen is that both OpenPaperOpt and Doxygen are developed under the GNU General Public License [14].

6.2 Naming Space

All the classes from the P3D are now defined in the namespace p3d, which is inspired from the definition of various class distributions that are defined under the namespace Distributions. The namespace could provide a more controlled mechanism for preventing name collisions. By using the namespace mechanism, for example, the different fiber model could be implemented under the same class name. This greatly increases the readability of the code and classes could have more meaningful name. In the future, if we want to add a network simulator, the developer should also consider having namespace with a proper name.

6.3 Network Simulator

A network simulator is to simulate the structure of the paper. Class Staticsheet is designed for this purpose. In a typical simulation an instance of Class Staticsheet is created and several member functions are called in turn to perform the necessary steps to create a sheet. These steps include:

1. Adding the fibers into the simulation volume, this operation could be repeated for several times.

2. Shuffle the fibers that have already added to the simulation volume.

3. Calculate the curl for each fiber.

(42)

4. Approximate the local angle for segments of each fiber.

5. Place the fibers in the plane of the simulation volume, i.e. estimate the X, Y coordinates for each fiber.

6. Estimate the depth of the fiber in the simulation with given surface, i.e. define the Z coordinate for each fiber. If without the surface, that means the generated the fiber network will place on a flat surface, otherwise it will be placed the on an uneven surface.

Each step could be achieved by a unique member function, this kind of design give the most flexibility to the user to define structure of the paper. In order to define a more realistic simulated paper, a new member function which is used to add the fibers to the simulation volume is implemented, this function make use of the distribution functions to generate the fibers with different radius, length and form factor. These properties of the fiber could be acquired from the distribution functions.

For cooperating the IO, a new member function which could save the generated fiber network into binary file and an XML description file are implemented. This is due to that a whole light scattering process on the static fiber network consists two steps. First OpenPaperOpt generate the fiber network according to the input specification files and save the fiber network into the files. Second OpenPaperOpt reads the fiber network from the files and performs light scattering on it with desired parameters. Another main reason to have this new function is that user could save the intermediated work at any time and resume the generating work later.

The binary file contains the detail information about each fiber in the simulation volume.

And the XML description file contains the information about the simulation volume, it includes the size of the generated the fiber network, the top surface and bottom surface for the fiber network.

6.3.1 Network Generated on an Uneven Surface

No matter the fiber network is generated on a flat or uneven surface, the essence of this operation is the same; it is to estimate the depth of each fiber in the simulation volume. Generating the fiber network on a flat surface has already implemented in the P3D, the rationale is to create the height buffer with the desirable resolution for the current simulation volume, the height buffer is an instance of class SurfAcc. Through it, we could find the lowest possible of placing a fiber and search the highest point in

(43)

use a height map with the associated resolution to initialize the height buffer. By using the member functions of the class SurfAcc we could find the depth for each fiber in the simulation volume. The higher resolution of the height map could generate more accurate and realistic result but it will also cost more running time. The height map could also be generated from the OpenPaperOpt, when generating the fiber network on flat or uneven surface, the surface on the top of the fiber network will be automatically generated and this surface could be used as height map and for generating the fiber network which on top of it.

6.3.2 Surface Generation

Each time when a fiber network is generated, there are two surfaces coming out together.

One is the surface on the top of the fiber network; another is the surface at the bottom of the network. Both of them are saved in the binary files with XML description files.

The binary files contain the value of the height points and the XML files contain the information like the binary file name, the resolution of the surface and so on. The resolution of the surface is calculated out automatically, it is determined by the size of the fiber network and the value of the probe size which is shown in the figure 6.1, currently the default probe size is 5.0mm. The XML description could be directly copy/paste to the input specification which is used to run the light scattering simulation.

Figure 6.1: Illustration of the resolution of the surface, this is an example of surface with resolution 4 by 4, the probe size is the maximum value of the length of each block,

the class SurfAcc could find out the highest and lowest point in each block.

(44)

6.4 General Light Scattering Procedure

A general light scattering procedure has been integrated into the OpenPaperOpt, this procedure is focusing on performing the light scattering on the static fiber network.

Figure 6.2: Class diagram showing the general program structure. The class Basesheet represents the layer in statistic description and the class Staticsheet is the layer in the

static description. The class Surface stands for the surfaces in the paper volume.

The figure6.2shows the general idea of program structure, different surfaces and layers in the paper volume are the subclasses of the abstract superclass StructureObject . The abstract StructureObject 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 surfaces or layers. But when dealing with the static fiber network, things are slight different. Class Staticsheet is the implementation of the network generator, and as we have discussed in the chapter 6, the general light scattering procedure in the OpenPaperOpt doesn’t directly work on the raw data. So the member method Process in the class Staticsheet is only a declaration without specific implementation. And the class StaticVolume is created for preparing the data which is generated by the class Staticsheet for the light scattering and the member method Process in it could invoke the general light scattering procedure. This kind of implementation not only reflects the idea of the object-oriented programming but also separating the network generator with the light scattering process well. Different classes have their own focusing.

Light scattering process starts with preparing the data for the ray tracing. The figure 6.3 illustrates the procedure of data preparation for the light scattering. This part of work has already been implemented and tested in the previous project [2]. The main data model used in data preparation procedure has already been presented in the chapter 6. Here is to give introduction of different classes to facilitate the future

(45)

Figure 6.3: Class diagram showing the relation between the Staticsheet class and the RndVolume class. The RndVolume contains the fiber and segments which are converted

from the counterparts in the Staticsheet.

the data preparation. So first an instance of class RndVolume is created according to the fiber network which is generated by the class Staticsheet. Then the instances of class Fiber and class FiberSegment are encapsulated as the instances of class RndFiber and RndFiberSeg respectively. To achieve this step, the class Fiber has the method which could retrieve the coordinates of each fiber segment. In order to perform the fast ray tracing the fiber segments in the class RndVolume will be put into bounding boxes which are implemented by the class OrientClipBox, this class will accelerate the ray tracing process to find the intersection points. The instance of the class RndVolume actually represents the fiber network which is suitable for the ray tracing procedure in the OpenPaperOpt.

The ray tracing is implemented in the class RndRay which is totally relied on the instance of the class RndVolume, that means if the fiber network could be successfully converted to the instance of the class RndVolume, the ray tracing procedure in the OpenPaperOpt could be performed.

Development of an Open Source Platform for Monte Carlo Simulation of Paper Structures

Examensarbete 30 hp April 2009

Development of an Open Source

Platform for Monte Carlo Simulation of Paper Structures

Jie Shang

Institutionen för informationsteknologi

Abstract

Development of an Open Source Platform for Monte Carlo Simulation of Paper Structures

Acknowledgements

Contents

Chapter 1

Introduction

1.1 Background

1.2 Project Objectives

1.3 Previous Work

1.4 Thesis Overview

Chapter 2

About Paper

2.1 Paper

2.2 Natural Fiber

Chapter 3

Light Scattering Model

3.1 Overview

3.2 Monte-Carlo Model

3.3 Paper Structure

3.4 Light Source

3.5 Detector

3.6 Light Scattering

Chapter 4

Fiber Network Generation

4.1 Staticsheet

4.2 Other Network Generation Model

Chapter 5

Light Scattering Process

5.1 Overview

5.2 Data Model

5.3 Data Preparation for Light Scattering

5.4 Optimization

Chapter 6

Implementation

6.1 Documentation Code

6.2 Naming Space

6.3 Network Simulator

6.4 General Light Scattering Procedure