How could radiation combined with hea... Log Out | Topics | Search
Moderators | Register | Edit Profile

FlexPDE User's Forum » User Postings » How could radiation combined with heat conduction be solved by FlexPDE? « Previous Next »

Author Message
Top of pagePrevious messageNext messageBottom of page Link to this message

Peter Andersson (speter)
New member
Username: speter

Post Number: 1
Registered: 10-2004
Posted on Monday, October 11, 2004 - 12:51 pm:   

The combined conduction and radiation problems (in one dimension) can be described by follow equations.

Heat conduction
dx(-k*dx(Temp))= -dqr/dx

(in 2D dy(-k*dy(temp)) is added)

where temp is the temperature, k is the thermal conductivity and qr is the radiation part.


Radiation
dqr/dx= -2*sigma*temp1^4*E2(x)-2*sigma*Temp2^4*E2(kd-x)-2*integral(sigma*[temp(u)]^4*E1( abs(x-u)))

where sigma is the Stefan-Boltzmann constant (5.67e-8 W/(m^2 K^4), E1 and E2 is exponential integral, temp1 and temp2 is temperature at the boundaries, kd is the length of domain.

[The boundary temperature temp1 and temp2 are values at x=0, x=kd.
For black walls the boundary fluxes are
q0,1=sigma*temp1^4/pi
q0,2=sigma*temp2^4/pi
For zero scattering is a gray medium is I(x)=sigma*T^4/pi (which is used in radiation equation above.)]



The integral should be carried out for the whole computaional region between 0 and kd for each x value.

What can I do to solve this integral which also includes the variable temp?
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 235
Registered: 06-2003
Posted on Monday, October 11, 2004 - 02:46 pm:   

Our example "Samples | Steady_State | Heatflow | Radflow.pde" shows the application of FlexPDE to problems in radiative transfer.

This model uses a first-order moments method (in angle) to model radiative energy density and couple it to the local temperature. Other higher-order angular approximations could be used to model more transparent media.

The attached plots show the results of this example, where energy is radiated from a hot slab to a cold slab across an intervening low-density gap.

Radiation Energy Density
Temperature
Top of pagePrevious messageNext messageBottom of page Link to this message

Peter Andersson (speter)
New member
Username: speter

Post Number: 2
Registered: 10-2004
Posted on Tuesday, October 12, 2004 - 11:52 am:   

I was looking for a solution using FlexPDE for combined radiation diffusion with heat conduction for intermediate optical thickness, where the optical thin and thick approximation provides limited accuracy.

Do you suggest that the combined radiation diffusion with heat conduction for intermediate optical thickness could be solved by a set of momentum equations by using for instance general spherical harmonics method?



Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 236
Registered: 06-2003
Posted on Tuesday, October 12, 2004 - 03:33 pm:   

I have not followed the subject for many years, but in the 60's and 70's there was an abundant literature on methods of approximation of radiative transfer.

The core problems in simulating radiative transfer are 1) resolution of the angular dependence of the radiation field, and 2) resolution of the spectral variations in absorption and emission.

1) One of the common techniques for angular resolution is the moments method: Multiply the radiation intensity equation by powers of the direction cosines and integrate over solid angle. (This is equivalent except for some algebra to expansion in spherical harmonics). One is free to keep as many terms as he feels necessary; the problem comes in the termination condition applied to truncate the series. The first three moments are energy density (a scalar), energy flux (a vector) and radiation pressure (a tensor). In thin media, it is sufficient to approximate the radiation pressure as equal to the energy density. In thick media, it is sufficient to approximate it as 1/3 the energy density. Approximations have been used to scale the proportionality as a function of opacity.

Another method which has been applied is to compute the propagation of radiation intensity at selected angles, namely the evaluation points for a Gauss quadrature in angle. Details differ from the moments method, but the spirit is similar.

2) One approach to the treatment of spectral dependence is to use an absorption average cross-section for emission and absorption, and a transmission averaged mean-free-path in the transport term. This process can be broken into frequency groups to give better resolution.

The method described in our example Radflow.pde uses the lowest-order implementation of the moments method, best described as "non-equilibrium diffusion". The example is easily extended to first-order by inclusion of a flux variable. The equation then become
E: dt(E) + dx(Fx) + dy(Fy) + c*sigmaa(T)*E = 4*pi*sigmap(T)*B(T)
Fx: dt(Fx) + c*beta*dx(E) + (1/lambdat)*Fx = 0
Fy: dt(Fy) + c*beta*dy(E) + (1/lambdat)*Fy = 0

Here B is the Planck radiation function integrated over the frequency group represented by the equation. In one-group ("grey") models, this is 4*pi*sigmap(T)*B(T) = a*c*sigmap(T)*T^4

I realize this is sketchy information, but I'm working from memory, and can't point you to any references. I hope it is a clue as to how to proceed.

By the way, all the things I have discussed are amenable to treatment in FlexPDE. The creation of averaged cross-sections is outside the scope of FlexPDE, and must be created in an exterior computation.

Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
New member
Username: vinodk1950

Post Number: 1
Registered: 06-2007
Posted on Friday, June 01, 2007 - 07:47 am:   

If radflow.pde is rewritten in new variables E, Fx and Fy, then, what Boundary conditions will be imposed on E, Fx and Fy to get solution identical to radflow.pde
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 865
Registered: 06-2003
Posted on Friday, June 01, 2007 - 02:51 pm:   

The equations I listed in my previous post are all first order, so there is no definition for NATURAL in any of the equations (see NATURAL in the Help Index).

You must therefore use VALUE conditions on the moments you wish to set. VALUE(E)=<expr> will set the radiation energy density. VALUE(Fx) will set the x-flux, etc.

In the RADFLOW example, NATURAL(Erad) is the outward normal component of beta*lambda*grad(Erad), or the inward surface-normal flux. This must be partitioned between Fx and Fy depending on the boundary direction.

It is unlikely that this formulation will give a solution identical to RADFLOW. Including the Fx, Fy equations creates a model with more angular information than that used in RADFLOW. If the solutions agree, then the Fx and Fy equations were unnecessary.
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
New member
Username: vinodk1950

Post Number: 2
Registered: 06-2007
Posted on Wednesday, June 06, 2007 - 05:59 am:   

I have written simplest discrete ordinates equations in x-y geometry to find the first eigen mode. I used the student key in window version of Flexpde4.15 (eigenvalue=1.234, correct), flexpde5.15 (eigenvalue=-0.318, wrong) and flexpde5.16(eigenvalue=-0.318, wrong). File is attached and also given below.
TITLE 'Discrete Ordinates-S2-Mode '

SELECT
gridlimit=1
ngrid=5
modes=1

VARIABLES
e1, e2, e3, e4 { flux in four directions}

DEFINITIONS
sigt=1
sigs=1
w=0.25
mu=0.57735
nu=0.57735
a=2.0 { half of square side}
phi=w*(e1+e2+e3+e4) { angle integrated flux }

EQUATIONS { equations for angular fluxes }
e1: div(mu*e1,nu*e1)+sigt*e1-sigs*phi*lambda=0
e2: div(-mu*e2,nu*e2)+sigt*e2-sigs*phi*lambda=0
e3: div(-mu*e3,-nu*e3)+sigt*e3-sigs*phi*lambda=0
e4: div(mu*e4,-nu*e4)+sigt*e4-sigs*phi*lambda=0

BOUNDARIES

REGION 1 'box'
START(-a,-a)
VALUE(e1)=0 VALUE(e2)=0 LINE TO(a,-a) {incoming flux zero along the boundary}
VALUE(e2)=0 VALUE(e3)=0 LINE TO(a,a) {incoming flux zero along the boundary}
VALUE(e3)=0 VALUE(e4)=0 LINE TO(-a,a) {incoming flux zero along the boundary}
VALUE(e4)=0 VALUE(e1)=0 LINE TO(-a,-a) {incoming flux zero along the boundary}


PLOTS
CONTOUR(Phi)
ELEVATION(Phi) FROM (-a,-a) to (a,-a)
ELEVATION(Phi) FROM (a,-a) to (a,a)
ELEVATION(Phi) FROM (a,a) to (-a,a)
ELEVATION(Phi) FROM (-a,a) to (-a,-a)
ELEVATION(Phi) FROM (-a,0) to (a,0)
ELEVATION(Phi) FROM (0,a) to (0,-a)
ELEVATION(Phi) FROM (-a,-a) to (a,a)
ELEVATION(Phi) FROM (a,-a) to (-a,a)
END
application/octet-streamrad2d-s2-mode.pde file
rad2d-s2-mode.pde (1.3 k)
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 873
Registered: 06-2003
Posted on Wednesday, June 06, 2007 - 06:24 pm:   

The rules for boundary conditions state that a given condition continues in force until changed, or until the path is complete.

This means, for example, that you have applied VALUE(E1)=0 around the entire perimeter.

I don't believe this is what you want, and it distorts the eigenvalues of the problem.

If you declare natural(e1)=0 on the second leg, natural(e2)=0 on the third leg, etc, you will get the result 1.2265 from version 5.0.16.

You can look at the plots and see whether the solutions are obeying the boundary conditions you want.
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 874
Registered: 06-2003
Posted on Wednesday, June 06, 2007 - 06:27 pm:   

Correction:

IF YOU HAD PLOTTED THE VARIABLES, you could have seen that the desired BC's were not met.

Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Junior Member
Username: vinodk1950

Post Number: 3
Registered: 06-2007
Posted on Wednesday, June 06, 2007 - 08:56 pm:   

Many thanks for your attention:
We have professional version of Flexpde2. and 3. Some papers are appearing on applying FEM to Discrete ordinates method in radiation transport. Why not Flexpde?
First, I understand that you mentioned that in first order equations, there are no Natural boundary conditions in Flexpde. However in part integration of first order Div operator the vector is available for Natural boundary conditions along the outer boundaries and that is what you mean now. Angular flux e1 points in north-east direction, e2 in north-west direction, e3 in south-west direction and e4 in south-east direction. Mathematically, the problem requires that there is no incoming radiation on the outer boundary of the convex body. This is applied as follows(exclude corners) (1) no incoming flux at lower boundary ( e1 & e2=0)(2) no incoming at the right boundary (e2 & e3=0) (3) no incoming at top boundary ( e3 &e4=0) and (4) no incoming at left boundary ( e4 & e1=0). Further, in FEM method special mention is needed at the four boundary corners and I had ignored them. In the bottom left corner no incoming means (e1, e2 and e4=0), in bottom right corner (e2,e3 and e1 =0), in top right corner (e3,e4 and e2=0) and in top left corner (e4, e1 and e3=0). I thought, the boundary condition is approximately visible from the inner contours of phi=0.25(e1+e2+e3+e4) in flexpde4. , flexpde3 & flexpde2. I shall investigate and discuss later
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 876
Registered: 06-2003
Posted on Wednesday, June 06, 2007 - 09:25 pm:   

I don't see any reason that FlexPDE cannot be used for Discrete Ordinate equations. As I recall, there are second-order forms of DO equations as well. (My books have not been unpacked since our move, so I can't check right now.)

It is true that for first-order equations, FlexPDE does not by default define the NATURAL BC. (This can be changed by using SELECT FIRSTPARTS, q.v.) But my suggestion of defining NATURAL BCs in your problem was simply to turn off the VALUE BCs. If you force E1=0 on all sides, then nothing can come in and nothing can go out. This is a problem.

You could use NOBC(E1) if you prefer, but the effect is the same.

I would have to check to be sure, but I think a VALUE condition will take precedence on any joint of edges on which it is applied. This should cover you special corner conditions automatically.

Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 4
Registered: 06-2007
Posted on Wednesday, June 06, 2007 - 11:55 pm:   

Thanks
Sorry to say, but I have not understood your statement that of e1 being set to zero on all sides. As, I have understood from my input, I have first set e1 & e2 zero along the bottom boundary, e2 & e3 zero along the right boundary, e3 & e4 zero along the top boundary and e4 & e1 zero along the left boundary.
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 5
Registered: 06-2007
Posted on Thursday, June 07, 2007 - 12:21 am:   

Many thanks
You meant to say, that unspecified variables in the boundary conditions are automatically set by NATURAL(unspecified variable)=0 and therefore NOBC(unspecified variable) is to be always used.
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 6
Registered: 06-2007
Posted on Thursday, June 07, 2007 - 01:18 am:   

Following is the corrected input

TITLE 'Discrete Ordinates-S2-Mode '

SELECT
gridlimit=1
ngrid=5
modes=1

VARIABLES
e1, e2, e3, e4 { flux in four directions}

DEFINITIONS
sigt=1
sigs=1
w=0.25
mu=0.57735
nu=0.57735
a=2.0 { half of square side}
phi=w*(e1+e2+e3+e4) { angle integrated flux }

EQUATIONS { equations for angular fluxes }
e1: div(mu*e1,nu*e1)+sigt*e1-sigs*phi*lambda=0 { e1- angular flux in north-east direction}
e2: div(-mu*e2,nu*e2)+sigt*e2-sigs*phi*lambda=0 { e2- angular flux in north-west direction}
e3: div(-mu*e3,-nu*e3)+sigt*e3-sigs*phi*lambda=0 { e3- angular flux in south-west direction}
e4: div(mu*e4,-nu*e4)+sigt*e4-sigs*phi*lambda=0 { e4- angular flux in south-east direction}

BOUNDARIES

REGION 1 'box'
START(-a,-a)
VALUE(e1)=0 VALUE(e2)=0 NOBC(e3) NOBC(e4) LINE TO(a,-a) {incoming flux zero along the bottom boundary}
NOBC(e1) VALUE(e2)=0 VALUE(e3)=0 NOBC(e4) LINE TO(a,a) {incoming flux zero along the right boundary}
NOBC(e1) NOBC(e2) VALUE(e3)=0 VALUE(e4)=0 LINE TO(-a,a) {incoming flux zero along the top boundary}
NOBC(e2) NOBC(e3) VALUE(e4)=0 VALUE(e1)=0 LINE TO(-a,-a) {incoming flux zero along the left boundary}


PLOTS
CONTOUR(Phi)
ELEVATION(Phi) FROM (-a,-a) to (a,-a)
ELEVATION(Phi) FROM (a,-a) to (a,a)
ELEVATION(Phi) FROM (a,a) to (-a,a)
ELEVATION(Phi) FROM (-a,a) to (-a,-a)
ELEVATION(Phi) FROM (-a,0) to (a,0)
ELEVATION(Phi) FROM (0,a) to (0,-a)
ELEVATION(Phi) FROM (-a,-a) to (a,a)
ELEVATION(Phi) FROM (a,-a) to (-a,a)
END
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 7
Registered: 06-2007
Posted on Thursday, June 07, 2007 - 01:22 am:   

Following is the corrected input

TITLE 'Discrete Ordinates-S2-Mode '

SELECT
gridlimit=1
ngrid=5
modes=1

VARIABLES
e1, e2, e3, e4 { flux in four directions}

DEFINITIONS
sigt=1
sigs=1
w=0.25
mu=0.57735
nu=0.57735
a=2.0 { half of square side}
phi=w*(e1+e2+e3+e4) { angle integrated flux }

EQUATIONS { equations for angular fluxes }
e1: div(mu*e1,nu*e1)+sigt*e1-sigs*phi*lambda=0 { e1- angular flux in north-east direction}
e2: div(-mu*e2,nu*e2)+sigt*e2-sigs*phi*lambda=0 { e2- angular flux in north-west direction}
e3: div(-mu*e3,-nu*e3)+sigt*e3-sigs*phi*lambda=0 { e3- angular flux in south-west direction}
e4: div(mu*e4,-nu*e4)+sigt*e4-sigs*phi*lambda=0 { e4- angular flux in south-east direction}

BOUNDARIES

REGION 1 'box'
START(-a,-a)
VALUE(e1)=0 VALUE(e2)=0 NOBC(e3) NOBC(e4) LINE TO(a,-a) {incoming flux zero along the bottom boundary}
NOBC(e1) VALUE(e2)=0 VALUE(e3)=0 NOBC(e4) LINE TO(a,a) {incoming flux zero along the right boundary}
NOBC(e1) NOBC(e2) VALUE(e3)=0 VALUE(e4)=0 LINE TO(-a,a) {incoming flux zero along the top boundary}
NOBC(e2) NOBC(e3) VALUE(e4)=0 VALUE(e1)=0 LINE TO(-a,-a) {incoming flux zero along the left boundary}


PLOTS
CONTOUR(Phi)
ELEVATION(Phi) FROM (-a,-a) to (a,-a)
ELEVATION(Phi) FROM (a,-a) to (a,a)
ELEVATION(Phi) FROM (a,a) to (-a,a)
ELEVATION(Phi) FROM (-a,a) to (-a,-a)
ELEVATION(Phi) FROM (-a,0) to (a,0)
ELEVATION(Phi) FROM (0,a) to (0,-a)
ELEVATION(Phi) FROM (-a,-a) to (a,a)
ELEVATION(Phi) FROM (a,-a) to (-a,a)
END
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 877
Registered: 06-2003
Posted on Thursday, June 07, 2007 - 02:59 pm:   

1. The default boundary condition when none is stated is NATURAL()=0.
2. Once stated, a boundary condition continues in force until explicitly changed, or until the path ends.
3. NOBC() may be used to terminate an applied condition. It is usually equivalent to NATURAL()=0.

See "Boundary Conditions -> Syntax" in the help index.
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 8
Registered: 06-2007
Posted on Friday, June 08, 2007 - 05:13 am:   

Thanks
I now understand and following are the boundary cards for the above problem

START(-a,-a)
VALUE(e1)=0 VALUE(e2)=0 LINE TO(a,-a)
NOBC(e1) VALUE(e3)=0 LINE TO(a,a)
NOBC(e2) VALUE(e4)=0 LINE TO(-a,a)
NOBC(e3) VALUE(e1)=0 LINE TO(-a,-a)

I now have a circular ( X-Y ) domain of radius R described below. How do I implement no incoming flux ( i.e. in terms of e1, e2, e3 and e4) on this circular outer boundary

REGION 1 'circle'
start(R,0) arc(center=0,0) angle =360
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 9
Registered: 06-2007
Posted on Friday, June 08, 2007 - 06:37 am:   

Following is the extension of boundary cards given above for square

REGION 1'circle'
start "Outer_Boundary" (-R/sqrt(2.),-R/sqrt(2.))
VALUE(e1)=0 VALUE(e2)=0
arc(center=0,0) to (R/sqrt(2.),-R/sqrt(2.))
NOBC(e1) VALUE(e3)=0
arc(center=0,0) to (R/sqrt(2.),R/sqrt(2.))
NOBC(e2) VALUE(e4)=0
arc(center=0,0) to (-R/sqrt(2.),R/sqrt(2.))
NOBC(e3) VALUE(e1)=0
arc(center=0,0) to (-R/sqrt(2.),-R/sqrt(2.))
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 879
Registered: 06-2003
Posted on Friday, June 08, 2007 - 03:15 pm:   

I assume that the "e" variables have the units of radiation intensity, in which case they are already fluxes, and a VALUE() BC would dictate flux.

Presumably these are directed fluxes, so you would need to weight the BC manually by boundary direction.

Alternatively, you can use the control SELECT FIRSTPARTS, which causes first-order terms to be integrated by parts. With this option, NATURAL(e1) will specify the outward boundary-normal component of (mu*e1,nu*e1), with similar rules for the other equations. This follows from the Divergence Theorem: Vol_Integral(DIV(Vector)) = Surf_Integral(normal<dot>Vector).

Notice that this is a global selector, and all first-order terms must be provided appropriate NATURAL() values (or defaulted to zero).
Top of pagePrevious messageNext messageBottom of page Link to this message

Johan H. Zietsman (johanzietsman)
New member
Username: johanzietsman

Post Number: 1
Registered: 08-2007
Posted on Tuesday, August 14, 2007 - 09:55 am:   

I am modelling a metallurgical smelting furnace and had a look at the Radflow example. I am not familiar with the methods used to model radiation heat transfer in this example. I am used to using a view-factor based approach to do it. However, view factors make it difficult to combine radiation and conduction modelling.

I would like to familiarise myself with the technique used in Radflow. Can you please direct to some literature references that thoroughly covers this topic?

Thank you for the help.
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 933
Registered: 06-2003
Posted on Tuesday, August 14, 2007 - 01:03 pm:   

The book I liked best years ago is "The Equations of Radiation Hydrodynamics" by Gerald C. Pomraning.
I haven't read it recently, but his chapter "Approximate Description of Radiative Transfer" seems to cover the subject in the way I remember.

It's available from Amazon for $16.

Top of pagePrevious messageNext messageBottom of page Link to this message

Johan H. Zietsman (johanzietsman)
New member
Username: johanzietsman

Post Number: 2
Registered: 08-2007
Posted on Thursday, August 16, 2007 - 02:56 am:   

Dear Robert. What is the method called that you use in the Radflow example? I will not get the Pomraning book soon, but I have access to other resources. If I know what the method/model is called, I can scan literature for it.

I have ran into the following models that may be possible to implement in FlexPDE:
- Discrete ordinates radiation model
- Rosseland radiation model
- P-1 radiation model
- Discrete transfer radiation model

I also know of the following models that may not be applicable to FlexPDE implementations:
- Zone radiation model (Hottel and Sarofim)
- Monte Carlo radiation model
- Radiosity model (I think it is similar to the zone model)
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 938
Registered: 06-2003
Posted on Thursday, August 16, 2007 - 02:29 pm:   

The equation I used in the Radflow example is a low-order Angular Moments method.

It is derived by multiplying the transfer equation by increasing powers of the direction cosines and integrating over angle. The first three moments are then energy density (E), energy flux (F) and radiation pressure (P).
Making simplifying assumptions about the flux and pressure allows the system to collapse to a single equation in energy density. Analysis of radiation pressure in the extreme cases of isotropic vs unidirectional radiation allows reduction of the pressure to either E/3 or E, respectively. Assuming that time derivative of flux is small allows reduction to Fick's Law
F=-beta*grad(E). The energy density equation can be combined with the material energy equation to determine material temperature. (In my Radflow example, I merely assumed the radiation was in equilibrium with the material temperature and ignored the material energy Cp*Temp. This is legal in steady-state, but not for time-dependence. I also ignored conduction.)

This is closely related to Spherical Harmonic (P-n) expansion, since the first two spherical harmonics are the same as the first two moments.

The issue about Rosseland and Plank mean cross-sections is fundamental to Gray-body approximation, and concerns methods of averaging frequency-dependent absorption and emission cross-sections.


Top of pagePrevious messageNext messageBottom of page Link to this message

Jacques Muller (jmuller)
New member
Username: jmuller

Post Number: 1
Registered: 08-2007
Posted on Friday, August 17, 2007 - 03:29 am:   

Mr Nelson,

relating to the topic of these posts, how do you suggest we can specify pde's to use the discrete ordinates method above for radiation and also include heat conduction?

Would temperature need to be specified as another variable and the standard conduction heat transfer equation be added? Or should the conduction term be added to the d.o. equations? How would the radiation density and temperature then be linked in either case?

We are modeling a furnace that has a radiative source in the center, radiating energy through a gaseous medium to the sides, where the energy is then conducted to outside.
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 940
Registered: 06-2003
Posted on Friday, August 17, 2007 - 06:02 pm:   

There are a lot of ways to configure radiation transport systems, and the
choice depends on what really needs to be modeled in your system.
Much of the exotic work in transport involves stellar atmospheres, where a
lot of nasty things go on that won't appear in a furnace.

1.
It sounds to me like perhaps the radiation source is rather distributed,
which implies no sharply peaked angular distributions. Also, I assume that
the system maintains a stable configuration for long periods of time, so you
don't have to worry about the propagation speed of radiation.

Under these conditions, you can probably use a variation of the moments
methods that I described before.
If the changes in temperature are slow, then it is reasonable to assume
that the radiation field will be in thermal equilibrium with the material
temperature (ie, absorption approximately balances emission at each point),
and the gas will be at the same temperature as the radiation passing through
it.

Under these simplifying assumptions, you can collapse it all to a single
equation in T by asserting
E = a*T^4 { the radiation energy density, a=4*sigma/c=7.5634e-15
erg/cm3/degreeK^4 }
F = -c*beta*LambdaR*grad(E) { vector flux, beta=1/3 in opaque
materials }
dt(Em) + dt(E) + div(F) = 0 { Em = material energy, rho*Cv*T }

Since the basic variable is T, you can include a conduction term in the Em
equation as well. If I were doing this, I would certainly start here, as it
is far the least work and probably is sufficient.

2.
The moments model can be generalized from this equilibrium diffusion model
if you want, by separating the E and T variables:
dt(E) + grad(F) + c*SigmaP*(E-a*T^4) = 0
dt(Em) = c*SigmaP*(E-a*T^4) {add conduction term, if you want }
I seriously doubt that E will deviate from a*T^4 for very long.

3.
If you want to use Discrete Ordinates, you will be able to model more kinds
of shadowing effects. This will be similar to (2), in that you will need to
carry a separate temperature variable. You will have one variable for each
ordinate, rather than a single E variable, and E will be replaced by a
summation over ordinates in the Em equation.
This model entials the most work.

Top of pagePrevious messageNext messageBottom of page Link to this message

Johan H. Zietsman (johanzietsman)
Junior Member
Username: johanzietsman

Post Number: 3
Registered: 08-2007
Posted on Friday, August 24, 2007 - 05:44 am:   

Dear Robert. Thank you for your assistance this far.

We have been able to set up a problem descriptor that solves combined radiation and conduction in a gas using the discrete-ordinates method. It is stable and solves quite quickly.

I now want to surround this gas domain with a wall through which conduction takes place. I would like to do this in a single problem descriptor.

I have 4 pdes for the four directions of radiation used in the S2 mode of the discrete-ordinates method. I also have one pde to handle conduction.

All the pdes should be active in the gas part of the domain, but only the conduction pde should be active in the wall part. I have tried this by eliminating (setting to zero) all terms except one in the radiation pdes. This results in very slow conversion rates.

How can I deactivate the radiation pdes in the wall domain without making the solution so unstable? Please help. I include the problem descriptor below.

Thank you again.

TITLE 'Combined Radiation-Conduction Model 2'

COORDINATES
CARTESIAN2

SELECT
NGRID = 5
ERRLIM = 1E-2

VARIABLES
I1, I2, I3, I4 { flux in four directions [W/m^2/Sr] }
T { temperature [K] }

DEFINITIONS
sigma = 5.669E-8 { Boltzmann's constant [W/m^2/K^4] }
sigmam = 0.1 { scattering coefficient of medium [1/m] }
kappam = 0.2 { absorption coefficient of medium [1/m] }
betam = sigmam + kappam { extinction coefficient of medium [1/m] }
epsilonm = 0.2 { emission coefficient of medium [1/m] }
phi = 1 { phase function }

{ Properties of Air }
kAir = 0.02227 ! W/m/K

{ Direction Ordinates }
mu1=-SQRT(0.5)
mu2=SQRT(0.5)
mu3=-SQRT(0.5)
mu4=SQRT(0.5)

xi1=-SQRT(0.5)
xi2=-SQRT(0.5)
xi3=SQRT(0.5)
xi4=SQRT(0.5)

eta=0.57735026

{ Direction Weigths [m^2] }
w=pi

{ In-scatter Summations [W/m^3/Sr] }
Inscatter = (w*I1 + w*I2 + w*I3 + w*I4) * phi * sigmam / (4*pi)

{ Surface Emissivity }
epsilon = 0.5

{ Surface Reflectivity }
rho = 0.5

{ Surface Temperatures [K] }
TS = 100 + 273
TN = 25 + 273
TW = (TS - (TS - TN)*Y)
TE = (TS - (TS - TN)*Y)

{ Heat Flux Vectors [W/m2] }
qr1 = w*VECTOR(mu1, xi1)*I1
qr2 = w*VECTOR(mu2, xi2)*I2
qr3 = w*VECTOR(mu3, xi3)*I3
qr4 = w*VECTOR(mu4, xi4)*I4
qr = qr1 + qr2 + qr3 + qr4

{ Absorption Heat Flux Vector [W/m3] }
qabsorb = kappam*w*(I1 + I2 + I3 + I4)

{ Emission Heat Flux Vector [W/m3] }
qemit = epsilonm*sigma*T^4

{ Radiation Heat Flow Rates on Walls [W] }
qrS = SINTEGRAL(NORMAL(qr), 'S Wall')
qrE = SINTEGRAL(NORMAL(qr), 'E Wall')
qrN = SINTEGRAL(NORMAL(qr), 'N Wall')
qrW = SINTEGRAL(NORMAL(qr), 'W Wall')
qrBalance = qrS + qrE + qrN + qrW

qemitfromgas = VOL_INTEGRAL(qemit, 'Cavity')
qabsorbfromradiation = VOL_INTEGRAL(qabsorb, 'Cavity')

qconduct = -kAir*grad(T)
qcS = SINTEGRAL(NORMAL(qconduct), 'S Wall')
qcE = SINTEGRAL(NORMAL(qconduct), 'E Wall')
qcN = SINTEGRAL(NORMAL(qconduct), 'N Wall')
qcW = SINTEGRAL(NORMAL(qconduct), 'W Wall')
qcBalance = qcS + qcE + qcN + qcW

qBalance = qrBalance + qcBalance

{ Switches }
absorbactive = 1
emitactive = 1
radiationactive = 1

{ Boundary Condition Functions [W/m2/Sr] }
I3S = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(xi1)*I1 + w*abs(xi2)*I2)
I4S = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(xi1)*I1 + w*abs(xi2)*I2)

I1N = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(xi3)*I3 + w*abs(xi4)*I4)
I2N = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(xi3)*I3 + w*abs(xi4)*I4)

I2W = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(mu1)*I1 + w*abs(mu3)*I3)
I4W = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(mu1)*I1 + w*abs(mu3)*I3)

I1E = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(mu2)*I2 + w*abs(mu4)*I4)
I3E = (epsilon*sigma*T^4)/(2*pi) + (rho/pi)*(w*abs(mu2)*I2 + w*abs(mu4)*I4)

PS = -epsilon*sigma*T^4 + (1-rho)*(w*I1 + w*I2) + NORMAL(qconduct)
PN = -epsilon*sigma*T^4 + (1-rho)*(w*I3 + w*I4) + NORMAL(qconduct)
PW = -epsilon*sigma*T^4 + (1-rho)*(w*I1 + w*I3) + NORMAL(qconduct)
PE = -epsilon*sigma*T^4 + (1-rho)*(w*I2 + w*I4) + NORMAL(qconduct)

xs = 0

EQUATIONS { equations for angular fluxes }
I1: div(mu1*I1, xi1*I1) + sigmam*I1 + kappam*I1 - epsilonm*(sigma*T^4)/(4*PI) - Inscatter + I1*xs = 0
I2: div(mu2*I2, xi2*I2) + sigmam*I2 + kappam*I2 - epsilonm*(sigma*T^4)/(4*PI) - Inscatter + I2*xs = 0
I3: div(mu3*I3, xi3*I3) + sigmam*I3 + kappam*I3 - epsilonm*(sigma*T^4)/(4*PI) - Inscatter + I3*xs = 0
I4: div(mu4*I4, xi4*I4) + sigmam*I4 + kappam*I4 - epsilonm*(sigma*T^4)/(4*PI) - Inscatter + I4*xs = 0
T: div(kAir*grad(T)) + absorbactive*qabsorb - emitactive*qemit = 0

BOUNDARIES

REGION 1 'Walls'
kAir = 1
absorbactive = 0
emitactive = 0
radiationactive = 0
epsilonm = 0
kappam = 0
sigmam = 0
mu1=0
mu2=0
mu3=0
mu4=0
xi1=0
xi2=0
xi3=0
xi4=0
xs = 1
START (0, 0)
VALUE(I1)=0 VALUE(I2)=0 VALUE(I3)=0 VALUE(I4)=0 VALUE(T) = TS LINE TO (1, 0)
VALUE(I1)=0 VALUE(I2)=0 VALUE(I3)=0 VALUE(I4)=0 VALUE(T) = TE LINE TO (1, 1)
VALUE(I1)=0 VALUE(I2)=0 VALUE(I3)=0 VALUE(I4)=0 VALUE(T) = TN LINE TO (0, 1)
VALUE(I1)=0 VALUE(I2)=0 VALUE(I3)=0 VALUE(I4)=0 VALUE(T) = TW LINE TO (0, 0)

REGION 2 'Cavity'
START(0.1, 0.1)
NATURAL(I1)=0 NATURAL(I2)=0 VALUE(I3)=I3S VALUE(I4)=I4S NATURAL(T) = PS LINE TO (0.9, 0.1)
VALUE(I1)=I1E NATURAL(I2)=0 VALUE(I3)=I3E NATURAL(I4)=0 NATURAL(T) = PE LINE TO (0.9, 0.9)
VALUE(I1)=I1N VALUE(I2)=I2N NATURAL(I3)=0 NATURAL(I4)=0 NATURAL(T) = PN LINE TO (0.1, 0.9)
NATURAL(I1)=0 VALUE(I2)=I2W NATURAL(I3)=0 VALUE(I4)=I4W NATURAL(T) = PW LINE TO (0.1, 0.1)

FEATURE 1 'S Wall' START (0, 0) LINE TO (1, 0)
FEATURE 1 'E Wall' START (1, 0) LINE TO (1, 1)
FEATURE 1 'N Wall' START (1, 1) LINE TO (0, 1)
FEATURE 1 'W Wall' START (0, 1) LINE TO (0, 0)

MONITORS
CONTOUR(error)

PLOTS
CONTOUR(T-273)
CONTOUR(qemit)
CONTOUR(qabsorb)
VECTOR(qconduct)
VECTOR(qr)
VECTOR(qr1)
VECTOR(qr2)
VECTOR(qr3)
VECTOR(qr4)
CONTOUR(I1)
CONTOUR(I2)
CONTOUR(I3)
CONTOUR(I4)
ELEVATION(I1, I2, I3, I4) ON 'S Wall'
ELEVATION(I1, I2, I3, I4) ON 'E Wall'
ELEVATION(I1, I2, I3, I4) ON 'N Wall'
ELEVATION(I1, I2, I3, I4) ON 'W Wall'
ELEVATION(NORMAL(qr)) FROM (1, 0) TO (1, 1)

SUMMARY
REPORT(qrS) AS "S Wall Radiation Heat Flow [W]"
REPORT(qrE) AS "E Wall Radiation Heat Flow [W]"
REPORT(qrN) AS "N Wall Radiation Heat Flow [W]"
REPORT(qrW) AS "W Wall Radiation Heat Flow [W]"
REPORT(qrBalance) AS "Balance of Radiation Heat Flow [W]"
REPORT(qcS) AS "S Wall Conduction Heat Flow [W]"
REPORT(qcE) AS "E Wall Conduction Heat Flow [W]"
REPORT(qcN) AS "N Wall Conduction Heat Flow [W]"
REPORT(qcW) AS "W Wall Conduction Heat Flow [W]"
REPORT(qcBalance) AS "Balance of Conduction Heat Flow [W]"
REPORT(qBalance) AS "Balance of Heat Flow [W]"
REPORT(VOL_INTEGRAL(epsilonm*sigma*T^4, 'Cavity')) AS "Emission loss from gas [W]"
REPORT(qabsorbfromradiation) AS "Absorption loss from radiation [W]"
END

Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 941
Registered: 06-2003
Posted on Friday, August 24, 2007 - 12:37 pm:   

FlexPDE solves all equations in all regions.

By setting the I coefficients to zero in the wall, you have made the equations 0=0, which is rather hard to solve for I.

But there is no reason to zero the I equations. There is, after all, a radiation field in the wall, it's just isotropic and equilibrated with the wall material temperature. If you give the wall large and equal absorption and emission coefficients, the I equations should result in qabsorb=qemit in the wall, leaving you with the conduction part as the active temperature transfer mechanism. Or, you can go ahead and zero qabsorb-qemit in the Temperature equation in the wall, but leave the I equations non-zero, so they can be solved.

If you really want to separate the wall material, you should use a CONTACT boundary with high resistance. Leave some coefficients active in the I equations, so they can be solved, even if the solution is 0.


Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 942
Registered: 06-2003
Posted on Friday, August 24, 2007 - 12:40 pm:   

Now that I think about it, you should always leave I active in the wall, because there will be a back-radiation from the wall due to the angular redistribution and re-radiation.
Top of pagePrevious messageNext messageBottom of page Link to this message

Johan H. Zietsman (johanzietsman)
Member
Username: johanzietsman

Post Number: 4
Registered: 08-2007
Posted on Thursday, August 30, 2007 - 08:19 am:   

I had huge instability problems when trying to solve combined radiation and conduction, so I am back doing radiation alone to try and verify the radiation part of the model. I have set the model up like it is described in a paper by Fiveland using the S2 discrete ordinates method. See his results below.

Fiveland results

If you run my problem descriptor, you will see that the GG elevation plots at x=0.1, 0.3 and 0.5 have some resemblance with the attached graphs from the Fiveland paper. However, my results contain what seems to be "discontinuities" compared with the published ones. When you look at the other graphs in my results it is quite clear that those feature will appear in the GG plots.

Are the "discontinuities" on the diagonals of my solutions caused by the triangular mesh that FlexPDE uses? The reason why I ask is because the S2 mode of the discrete ordinates method uses directions that are in line with the two diagonals of the square.

If this is the case, I will have to try and use a rectangular grid. Is FlexPDE capable of generating meshes other than triangular in 2D?

I attach my problem descriptor below:
application/octet-streamDOM Radiation Problem Descriptor
Discrete Ordinates Radiation Model 001.pde (4.9 k)


Thank your for your help.
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 944
Registered: 06-2003
Posted on Thursday, August 30, 2007 - 02:30 pm:   

Looks to me like it's telling you the truth.

Take I3, for example. This beam travels diagonally up-left. The bottom value is Eb/pi and the right wall value is 0. So you have a sharply shadowed beam, with the upper right half dark.

Similarly I4, travelling up-right sharply shadowed.

I1 and I2 start out zero, but are increased by scatter from I3 and I4.

The mesh tries to refine along the shadow lines to resolve the (discontinuous) dropoffs.

All perfectly legitimate.

I don't know where the reference plots come in, they don't seem to be addressing the same problem. The P3 errors could be caused by improper termination of the Pn sequence.

Top of pagePrevious messageNext messageBottom of page Link to this message

Johan H. Zietsman (johanzietsman)
Member
Username: johanzietsman

Post Number: 5
Registered: 08-2007
Posted on Thursday, August 30, 2007 - 03:56 pm:   

I tend to agree with you, in that the results are not far from what I expected. However, I have to compare the solution with a published one (e.g. the one from Fiveland.

I will explain the graphs. There are four series on each. The first series is from the zone method, and the second from the P3 method. The last two were done with the discrete ordinates method that I am using. I used the S2 configuration of the method. In other words, the results from my descriptor must be compared with the S2 results on each of the graphs.

What is concerning me is the fact that curves from FlexPDE are a lot more irregular than the published ones.

When Fiveland describes how he executed the numerical solution, he explains that he used a specific set of initial values, and then used a very specific sequence with which he recalculated the nodal values to work towards convergence. This makes sense when one considers the interdependencies between the various I terms.

I don't know how FlexPDE calculates the numerical solution, but because it is a general-purpose solver it is quite certain that it is not done in the way that Fiveland described. Could this have an influence the results that it arrives at?

I am assuming that Fiveland's published results are a good reference, since they, in turn, correspond well with that of other workers. It is the discrepancy between my results and his that concerns me.

Thank you again for the help.

P.S. Is FlexPDE able to create meshes other than triangular?
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 945
Registered: 06-2003
Posted on Friday, August 31, 2007 - 01:26 pm:   

Since the FlexPDE solution conforms to logic about what the answer should be, it is pointless to conjecture what details of the solution process might be tweaked to get a different answer.

The solution you are getting is the one dictated by your equations and boundary conditions. If the solution does not agree with a published result, then you are most likely not solving the same problem.

Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 946
Registered: 06-2003
Posted on Friday, August 31, 2007 - 04:06 pm:   

You could make the sidewalls periodic. Then it would simulate a 1D propagation in Y.
Exiting I3 on the left would feed in as entering I3 on the right, etc.
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 10
Registered: 06-2007
Posted on Monday, September 03, 2007 - 12:46 pm:   

Dear Johan
In Discrete_Ordinates_Radiation_Model_001.pde, set Eb as zero, kappam as positive and rhow as non-negative to get normal results. In Flexpde, non-zero external source (Eb) in the incoming boundary condition in the South face does not seem to work. This problem can be circumvented by adding a thin medium at the prescribed temperature. In fact earlier, I also encountered this problem and struggled to replace it as above. In short, Flexpde does not incorporate inhomogeneous boundary conditions in discrete ordinates method. This problem should be corrected in the software.
Vinod Kumar
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 947
Registered: 06-2003
Posted on Tuesday, September 04, 2007 - 02:02 pm:   

Vinod -

You'll have to clarify what you mean by this post.
I see no failure in applying Eb at the south boundary.
Furthermore, FlexPDE does not implement the Discrete ordinates method, the user does in his script.

Clarify what you mean by "FlexPDE does not incorporate inhomogeneous boundary conditions...".
I am not aware of any boundary conditions you cannot apply in FlexPDE.


Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 11
Registered: 06-2007
Posted on Wednesday, September 05, 2007 - 06:38 am:   

Dear Nelson
No adverse comments meant for Flexpde. I have the following observations:
(1) The treatmeant of the boundary conditions at the corner nodes, when incoming vector is non zero because of diffused reflection and/or external sources ( as is the case in Johan's problem).
(2) The oscillations seen in plots and particularly of I3 & I4 at the North boundary ( are these oscillations result of lagrangian interpolation of nodal values). If this is true, then all the plots should be carefully interpreted.
(3) All standard discrete ordinates scheme are based on sweeping the domain along various directions and these include a set to zero fix-up of any negative directional component till convergence over all directions is achieved. On the other hand, FEM solvers ( including Flexpde) assemble the full matrix and obtain the solution vectors by using an appropriate solver and you might end up with some results different from standard discrete ordinates method.
(4) By replacing the following boundary cards ( reflecting boundary conditions along the east and west walls) in Johan's problem, one gets 1-D results as correct straight lines
START(0, 0)
NATURAL(I1)=0 NATURAL(I2)=0 VALUE(I3)=I3S VALUE(I4)=I4S
LINE TO(a, 0)
VALUE(I1)=I2 NATURAL(I2)=0 VALUE(I3)=I4 NATURAL(I4)=0
LINE TO(a, b)
VALUE(I1)=0 VALUE(I2)=0 NATURAL(I3)=0 NATURAL(I4)=0
LINE TO(0, b)
NATURAL(I1)=0 VALUE(I2)=I1 NATURAL(I3)=0 VALUE(I4)=I3
LINE TO(0, 0)
Thanks for your comments.
Vinod
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 948
Registered: 06-2003
Posted on Wednesday, September 05, 2007 - 02:13 pm:   

1) This item has a subject but no predicate, so I can only guess. It is true that you can only define one VALUE BC at any node, and there is an ambiguity at a corner where you intend to apply a discontinuity. Which of two conflicting BC's do you apply there? By default, FlexPDE applies the one that occurs last in walking the boundary.

But we sternly warn against discontinuous BC's in the first place, as discontinuities are not representable by polynomial interpolation.

Unless the software has rotted, VALUE BC's should always take precedence over NATURAL's at such corner. A VALUE meeting a NATURAL is not really a discontinuity, since the NATURAL applies to the boundary leg, not to a point.

2) Using your modified BCs (item 4) I don't see any oscillations in I3 and I4 at the North boundary (see attached graphic). Oscillations in the original formulation were because the BCs demanded a discontinuous field, and discontinuities cannot be interpolated by the low-order polynomial interpolation fundamental to the finite element method. You could, of course, use a small ERRLIM or a FRONT statement to force dense meshing at the discontinuities and thereby limit the range (but not the amplitude) of the oscillations. A tiny diffusion term would also smear the discontinuity and eliminate the ringing.

3) True.

4) This BC definition correctly represents reflective walls, and turns the 2D solution into a 1D problem in Y. Periodic sidewall boundaries should have the same result.

I3I3
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 12
Registered: 06-2007
Posted on Thursday, September 06, 2007 - 01:14 am:   

Dear Nelson
Thanks for your graphics. This example is important in following ways.
(1) The planar geometry homogeneous linear transport equation with pure scattering admits two simple solutions (a) I(x, mu)=c1, and (b) I(x, mu)=c2(x-mu/sigma), where mu*dI(x, mu)/dx + sigma*I(x, mu)=(sigma/2)*Integral(I(x,mup)d(mup)), the integral is from -1<mup<1, I(x, mu) is angular radiation density, x is spatial variable and mu & mup are angular variables; c1 & c2 are two constants to be determined by taking the linear combination of (a) and (b) and satisfying the boundary values. ( verify (a) and (b) by substitution in Equation)
(2) The above equation is being programmed in Johan’s script with the reflecting boundary condition given in my previous post. Only difference is that we are using S2 approximations ( I1 &I2 represent mu=-0.57735, while I3 & I4 represent mu=0.57735 directions). The conclusions in (1) given above continue to hold good for S2-Discrete Ordinates Method. A linear combinations of above two solutions (a) and (b) is determined by satisfying the boundary conditions on south and north walls. This is an exact solution to S2 equations.
(3) Flexpde is able to reproduce this result after quite a few number crunching and thus Johan’s script with reflecting boundary conditions on east and west walls matches Flexpde results with exact solution.
Vinod Kumar
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 13
Registered: 06-2007
Posted on Thursday, September 06, 2007 - 01:17 am:   

Dear Nelson
Thanks for your graphics. This example is important in following ways.
(1) The planar geometry homogeneous linear transport equation with pure scattering admits two simple solutions (a) I(x, mu)=c1, and (b) I(x, mu)=c2(x-mu/sigma), where mu*dI(x, mu)/dx + sigma*I(x, mu)=(sigma/2)*Integral(I(x,mup)d(mup)), the integral is from -1<mup<1, I(x, mu) is angular radiation density, x is spatial variable and mu & mup are angular variables; c1 & c2 are two constants to be determined by taking the linear combination of (a) and (b) and satisfying the boundary values. ( verify (a) and (b) by substitution in Equation)
(2) The above equation is being programmed in Johan’s script with the reflecting boundary condition given in my previous post. Only difference is that we are using S2 approximations ( I1 &I2 represent mu=-0.57735, while I3 & I4 represent mu=0.57735 directions). The conclusions in (1) given above continue to hold good for S2-Discrete Ordinates Method. A linear combinations of above two solutions (a) and (b) is determined by satisfying the boundary conditions on south and north walls. This is an exact solution to S2 equations.
(3) Flexpde is able to reproduce this result after quite a few number crunching and thus Johan’s script with reflecting boundary conditions on east and west walls matches Flexpde results with exact solution.
Vinod Kumar
Top of pagePrevious messageNext messageBottom of page Link to this message

Vinod Kumar (vinodk1950)
Member
Username: vinodk1950

Post Number: 14
Registered: 06-2007
Posted on Tuesday, November 13, 2007 - 12:30 am:   

Dear Nelson
As mentioned in help, Flexpde uses the method of subspace iteration to solve for selected number of eigenvalues ( bathe & Wilson, 1976). In discrete ordinates method, the matrix formed will not be symmetric. What type of modified scheme are you using for unsymmetric matrix.
Top of pagePrevious messageNext messageBottom of page Link to this message

Robert G. Nelson (rgnelson)
Moderator
Username: rgnelson

Post Number: 991
Registered: 06-2003
Posted on Tuesday, November 13, 2007 - 01:14 pm:   

I don't remember if there any assumptions of symmetry in the basic reduction operations
A* = YtAY and B* = YtBY
But the solution AY = BX and the reduced eigenvalue problem A*Q = B*Q*lambda are both done for general matrices.

I will guess then that the solution will work and produce the right-eigenvectors of the system.

FlexPDE assumes that the eigenvalues are real, so that might be a problem.

Add Your Message Here
Post:
Username: Posting Information:
This is a private posting area. Only registered users and moderators may post messages here.
Password:
Options: Enable HTML code in message
Automatically activate URLs in message
Action:

Topics | Last Day | Last Week | Tree View | Search | Help/Instructions | Program Credits Administration