Neumann boundary conditions for thermal problem

[helmut....@ingenieurbuero-mbi.de]

Hello to all members of the list, escpecially to the inventors of deal.ii !

First of all I´d like to remark that this lib seems to be very impressive and usefull. 

At the moment I evaluate it to decide if we could use it for our calculations in building physics or not. So I worked through some examples, and tried to implement some stuff I use in my daily business. 

We have very often to calculate 2 and 3 dimensional problems with thermal conductivity and a solution dependant flux through the boundarys, due to ambiente temperature and a given heat transfer coefficient, eg 20°C and 0.25 m^2K/W inside and lower temperatures and a lower surface resitance on the outer boundaries. 

At the moment I´m not able to assemble the system in a way that this effect ist calculated the right way. The only solution I have found ist to solve the problem iterative by calculating the flux by hand and to adjust neumann value step by step. As you can imagine this is a true performance killer.

So my question at the moment is: is there a way to assemble such solution dependant quantities direct into the system-matrix or into the system-rhs to avoid the itterations ?

Unfortunately I did not find any example dealing with such bc that was understandable to me. With our current solution (getDP) there were no problems to assemble the system in such a way, so I assume there must also be a way in deal.II I just don´t see.

Any help would be very appreciated, thank you in advance

sincerly Helmut

Sorry for my bad wording and spelling, but I finished school a long time ago...

[Guido Kanschat]

Dear Helmut,

what you are describing is usually referred to as mixed or Robin boundary condition. Please have a look at the recent thread on albedo boundary conditions.

Best,
Guido

[helmut....@ingenieurbuero-mbi.de]

Dear Guido,

thank you for your answer and your time,

this was an Approach I have already tried, but unfortunately without any success...

my setup looks like this

for (; cell != endc; ++cell)

      {

               fe_values.reinit(cell);

               cell_matrix = 0;

               cell_rhs = 0;

               for (unsigned int i = 0; i < dofs_per_cell; ++i)

                       for (unsigned int j = 0; j < dofs_per_cell; ++j)

                               for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)

                                       cell_matrix(i, j) += (fe_values.shape_grad(i, q_point)  
                                                        * fe_values.shape_grad(j, q_point)

                                                       * fe_values.JxW(q_point));

to assemble the cell related part of the system-matrix (lambda as heat conductivitie set to 1, so it does not appear )

Then, on the boundary 1 I tried:

for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face)

                       if (cell->face(face)->at_boundary() && (cell->face(face)->boundary_indicator() == 1))

                       {

                               fe_face_values.reinit(cell, face);

                               for (unsigned int i = 0; i < dofs_per_face; ++i)

                                       for (unsigned int j = 0; j < dofs_per_face; ++j)

                                               for (unsigned int q_point = 0; q_point < n_face_q_points;
                                                                                                ++q_point)

                                               {

                                                       // on boundary 1 : h*(theta_s - theta_e)

                                                       cell_matrix(i, j) += 1 / R_SI_1

                                                        * ((fe_face_values.shape_value(i, q_point)

                                                        * fe_face_values.shape_value(j, q_point))

                                                         * fe_face_values.JxW(q_point) - THETA_1);

                                               }

                       }

assuming that the part 

((fe_face_values.shape_value(i, q_point) * fe_face_values.shape_value(j, q_point)) * fe_face_values.JxW(q_point)

would represent the (later) calculated Temperature on the surface. But this did not work for me. 

With lambda as the conductivitie set to 1, a mesh  ( square 2 units by 2 units), robin bc with r_si = 0.25 and theta_1 = 10 degrees, opposite edge bound ( by dirichlet to 1 degree :

std::map<unsigned int, double> boundary_values;

VectorTools::interpolate_boundary_values(dof_handler, 0, ConstantFunction<2>(1), boundary_values);

MatrixTools::apply_boundary_values(boundary_values, system_matrix, solution, system_rhs);

)

 the boundary with robin conditions should be pulled up to 7 degree (analytical solution) , but it drops below 0...

Does anybody has an example program, that simply solves a heat conduction problem with those boundary conditions (convektive heat transfer, described by coefficent and ambient temperature ?

Sincerly,

Helmut

[Andrew McBride]

Hi

why is the the term involving the external temperature included in the system matrix? Should it not be included in the RHS (i.e. the force vector)?

Cheers

Andrew

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[helmut....@ingenieurbuero-mbi.de]

Hi Anrew,

I tried to do it this way, too. Unfortunately without success.

That was the reason why I asked for a simple example someone could have on his (or her) disk.

I´m just looking for a Programm that solves the problem on a simple region, eg a square just to see (and learn) how it is done.

My itterative solution does indeed apply to rhs, but in the thread that Guido has mentioned, this is assembled (as far as I see) into the cell_matrix and then into the system_matrix.

Thanks, I´ll give it further tries (brute force attack)

Cheers, Helmut

[helmut....@ingenieurbuero-mbi.de]

Sorry, obviously I have misunderstood something:

my working itterative solution looks like this:

neumann = 1.0 / R_SI_1 * (THETA_1 - myOldValues[q_point]); // calculated outside loop

...
            for (unsigned int i = 0; i < dofs_per_cell; ++i)

              cell_rhs(i) += neumann * fe_face_values.shape_value(i, q_point) * fe_face_values.JxW(q_point);

and this at least works. So the term fe_face_values.shape_value(i, q_point) * fe_face_values.JxW(q_point) can´t be the value of the solution that I require.

Does anybody have a working example for stupid newbies like me ?

Does anybody solve problems like mine ?

I thought it would be a common problem, as fas as robin bc are the common case for calculations in building physics...

Sincerly,

Helmut

[helmut....@ingenieurbuero-mbi.de]

Ok,

it looks like I would have found the solution :

for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face)

                        if (cell->face(face)->at_boundary() && (cell->face(face)->boundary_indicator() == 1))

                        {

                                fe_face_values.reinit(cell, face);

                                for (unsigned int i = 0; i < dofs_per_cell; ++i) //

                                        for (unsigned int j = 0; j < dofs_per_cell; ++j)

                                                for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)

                                                        cell_matrix(i, j) += 1 / R_SI_1* fe_face_values.shape_value(i, q_point) * fe_face_values.shape_value(j, q_point) * fe_face_values.JxW(q_point);

                                for (unsigned int j = 0; j < dofs_per_cell; ++j)

                                        for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)

                                                cell_rhs(j) += 1/ R_SI_1 * THETA_1 * fe_face_values.shape_value(j, q_point) * fe_face_values.JxW(q_point);

                        }

this does the job. I had to add contribution to rhs and system_matrix ....

Now it does what it should and I´ll proceed with the evaluation.

Never the less, does nobody else use those Robin-BCs ?

Cheers, Helmut

[Guido Kanschat]

Dear Helmut,

even if your boundary term is on the face, you have to do the integration for all dofs per cell, not just dofs per face. Could you first try with zero external temperature? Indeed, as Andrew has pointed out, the external temperature should be on the right hand side. I hope that works as well, if you run the loop over all dofs.

Let me know if that fixes it!

Best,

Guido

[helmut....@ingenieurbuero-mbi.de]

Hi Guido,

as you can see, the working version does integrate over the dofs of the whole cell. So I think my problem is solved by the last code snippet I have posted. It works like a charm, never the less, perhaps it would be helpfull if you or the other gurus could have a look to confirm that this is the right way... The adaptive refinement of the mesh was an easy part after this. I love it, as I assume that it will save me a lot of time during future calculations.

By the way, I have learned programming about 25 years ago, avoided c++ until now, and I´m very impressed by the beauty of the code of the deal library and the published examples. I have to learn a lot !

Cheers, 

Helmut

[Thomas Wick]

Hallo Helmut,

I also did this face integration some time ago. If you think of
future extensions of your code you should be careful at one point.

The FE-matrix (also on faces) assembling reads:

a_{ij} = a(\phi_j, \phi_i)

Thus, the indices must be switched. This no problem for symmetric
equations like Laplace.

However, if you have non-symmetric equations like I do, then it plays a
role. Therefore, I would only
recommend to switch i and j in the matrix entries local_matrix(j,i)
instead of local_matrix(i,j), such that (my example):


for (unsigned int i=0; i<dofs_per_cell; ++i)

for (unsigned int j=0; j<dofs_per_cell; ++j)

{

local_matrix(j,i) -= (phi_v[i] * phi_v[j]
) * fe_face_values.JxW(q);

}


Best Thomas

--
++-------------------------------------------------------++
Dr. Thomas Wick
The University of Texas at Austin
Institute for Computational Engineering and Sciences (ICES)
201 East 24th Street
ACE 5.332 | Campus Mail C0200
Austin, TX 78712, USA

Email: tw...@ices.utexas.edu
Tel.: +1 (512) 232-7763
www: http://numerik.uni-hd.de/~twick/
++-------------------------------------------------------++
--

> --
> The deal.II project is located at http://www.dealii.org/ [1]

> For mailing list/forum options, see

> https://groups.google.com/d/forum/dealii?hl=en [2]

> ---
> You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it,
> send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out

> [3].
>
>
>
> Links:
> ------
> [1] http://www.dealii.org/
> [2] https://groups.google.com/d/forum/dealii?hl=en
> [3] https://groups.google.com/groups/opt_out

--
++-------------------------------------------------------++
Dr. Thomas Wick
The University of Texas at Austin
Institute for Computational Engineering and Sciences (ICES)
201 East 24th Street
ACE 5.332 | Campus Mail C0200
Austin, TX 78712, USA

Email: tw...@ices.utexas.edu
Tel.: +1 (512) 232-7763
www: http://numerik.uni-hd.de/~twick/
++-------------------------------------------------------++
--

[Wolfgang Bangerth]

> The FE-matrix (also on faces) assembling reads:
>
> a_{ij} = a(\phi_j, \phi_i)
>
> Thus, the indices must be switched. This no problem for symmetric equations
> like Laplace.
>
> However, if you have non-symmetric equations like I do, then it plays a role.
> Therefore, I would only
> recommend to switch i and j in the matrix entries local_matrix(j,i)
> instead of local_matrix(i,j), such that (my example):
>
>
> for (unsigned int i=0; i<dofs_per_cell; ++i)
> for (unsigned int j=0; j<dofs_per_cell; ++j)
> {
>
> local_matrix(j,i) -= (phi_v[i] * phi_v[j]
> ) * fe_face_values.JxW(q);
>
> }

I always teach my students to derive the bilinear form by multiplying the
strong form by a test function *from the left*. This avoids the confusion with
having to transpose the matrix.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Thomas Wick]

> I always teach my students to derive the bilinear form by multiplying
> the strong form by a test function *from the left*. This avoids the
> confusion with having to transpose the matrix.
>

This is of coarse the right idea.
But unfortunately, there are many classes and books who do not so.

And then you are asking several days why your results are wrong ... ;)

Best Thomas

> Best
> W.
>
>
> --
>
> ------------------------------------------------------------------------
> Wolfgang Bangerth email:
> bang...@math.tamu.edu
> www:
> http://www.math.tamu.edu/~bangerth/

[helmut....@ingenieurbuero-mbi.de]

Hi Thomas,

thank you for the remark. As I only intend to solve laplace in the context of building physics, I see no problem there. But of course I´ll change the order as you mentioned, just to be correct. 

That seems to be the problem for me : I do not know enough of the background math, I´m just a user who wants to solve well defined problems. So normally it is not necessary for me to look behind the scene. I just have to answer the questions judges have, useing the tools I have. But at the moment I look around to find new tools...

Sincerly 

Helmut

[Wolfgang Bangerth]

> |
> for(unsignedintface =0;face <GeometryInfo<2>::faces_per_cell;++face)

> if(cell->face(face)->at_boundary()&&(cell->face(face)->boundary_indicator()==1))
> {
> fe_face_values.reinit(cell,face);
> for(unsignedinti =0;i <dofs_per_cell;++i)//
> for(unsignedintj =0;j <dofs_per_cell;++j)
> for(unsignedintq_point =0;q_point <n_face_q_points;++q_point)
> cell_matrix(i,j)+=1/R_SI_1*fe_face_values.shape_value(i,q_point)*fe_face_values.shape_value(j,q_point)*fe_face_values.JxW(q_point);
> for(unsignedintj =0;j <dofs_per_cell;++j)
> for(unsignedintq_point =0;q_point <n_face_q_points;++q_point)
> cell_rhs(j)+=1/R_SI_1*THETA_1*fe_face_values.shape_value(j,q_point)*fe_face_values.JxW(q_point);
> }
> |

This looks correct to me.

> this does the job. I had to add contribution to rhs and system_matrix ....
>
> Now it does what it should and I�ll proceed with the evaluation.
>
> Never the less, does nobody else use those Robin-BCs ?

Yes, it's quite common in cases where the flux across the boundary depends on
the value of the solution. Examples are osmosis through a membrane, the flux
of light from a turbid medium into air or vacuum, or (a nonlinear example) the
Planck radiation condition of a hot body.

As a general remark -- you'd have found your solution quicker if you write
down the weak form that goes with your solution, integrate by part, and bring
all parts that are known to the right hand side and leave everything that
contains the solution on the left. This way you see what should go into the
matrix and which should go into the right hand side.

[helmut....@ingenieurbuero-mbi.de]

Hello Wolfgang,

thank you for the time reading and answering to me.

Of course, you are right. It would have been the better way to derive the solution from the problem, as you mentioned. But unfortunately, PDEs are not my native language, my time at the university has happened in the 90th of the last century, so working through weak and strong formulation would have costed me more time than just trial and error.

At the moment I spend a lot of time with deal, and I think I learn a lot. Not only about c++, but also about the FE-method that I frequently use but not fully understand.

So, have a nice weekend, and thank you again.

Cheers, Helmut