Steady state thermal problem formulation for ISO 10211
pb...@o2.pl
I try to implement a solution for steady state thermal conditions as required by the ISO 10211 norm. There has been already 3 threads on this subject started in 2013 by Helmut Müller. The last one ends with a very promising post (from Wolfgang) suggesting a final answer, but... nothing follows, unfortunately.
Anyway, the problem is basically to calculate steady state temperatures inside a multi-layer wall, where each layer has a different conduction coefficient, and internal/external surfaces have their own convection coefficients. We have outside and inside temperatures given (and the mentioned coefficients of course).
According to Wolfgang’s lectures I started with a basic physical equation, which I suppose for this kind of problem is:
-c∆T = 0
where c is a coefficient and T is temperature.
So, the weak form would be something like:
c*∫grad(fi)*grad(T) – c*∫fi*δT = 0
Ω δΩ
Now, I need to split the boundary to the inner and outer one, and I just do nothing to the 2 other boundaries which form the shorter edges of the rectangle which my wall is (Helmut called them adiabatic). After replacing the c with the given conduction/convection coefficients (K and h) and replacing δT with (T – Tin) and (T – Tout), which I suppose could be allowed; I get the final weak form, which after moving the known part to the rhs is:
K*∫grad(fi)*grad(T) – h1*∫fi*T – h2*∫fi*T = -h1*Tin*∫fi – h2*Tout*∫fi
Ω δΓ1 δΓ2 δΓ1 δΓ2
And this is exactly what appears in Helmut’s code, except that he does not have any minuses before the boundary integrals.
Nevertheless, with or without minuses the results are wrong. And here is the code I use to build my matrix and rhs: DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end();
for( ; cell != endc; cell++ )
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit( cell );
for( unsigned int q_index = 0; q_index < n_q_points; ++q_index )
for( unsigned int i = 0; i < dofs_per_cell; ++i )
for( unsigned int j = 0; j < dofs_per_cell; ++j )
cell_matrix( i, j ) += ( (Karray[cell->material_id( )]) * fe_values.shape_grad( i, q_index ) *
fe_values.shape_grad( j, q_index ) * fe_values.JxW( q_index ) );
for( unsigned int f=0; f < GeometryInfo<2>::faces_per_cell; ++f )
if( cell->face( f )->at_boundary( ) && cell->face( f )->boundary_id( ) > 0 )//== 2 )
{
double h = 2 == cell->face( f )->boundary_id( ) ? hout : hin;
double T = 2 == cell->face( f )->boundary_id( ) ? Tout : Tin;
for( unsigned int q_index=0; q_index < n_face_q_points; ++q_index )
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += ( h * fe_face_values.shape_value( i, q_index ) * fe_face_values.shape_value( j, q_index ) *
fe_face_values.JxW (q_index) );
for( unsigned int q_index = 0; q_index < n_face_q_points; ++q_index )
for( unsigned int i = 0; i < dofs_per_cell; ++i )
cell_rhs( i ) += ( h * T * fe_face_values.shape_value( i, q_index ) * fe_face_values.JxW( q_index ) );
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global( cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs );
I provided it in case I made a mistake translating the weak form to the deal II code. Anyway, as I mentioned, the results are wrong. That is for example:
Number of active cells: 7254
Total number of cells: 10211
582 CG iterations needed to obtain
convergence.
Temperature at point A (7.1) : 4.00151
Temperature at point A (7.1) : 4.00151
Temperature at point B (0.8) : 0.236359
Temperature at point C (7.9) : 5.43241
Temperature at point D (6.3) : 3.84004
Temperature at point E (0.8) : 0.310349
Temperature at point F (16.4) : 18.5076
In brackets there are given the correct results, while mine afterwards. I tried many things. I triple-checked every part of the code: triangulation, finer triangulation, boundary ids, material ids, solvers, etc. In fact I’ve spent the last 3 weeks trying to fix it myself. Without success, as one might deduce from my post…
Any help would be greatly appreciated !
Piotr