Thermoelastic Problem

[Lisa Collins]

Hi,

I want to start deal.ii and model a thermoelastic problem in 2D & 3D.

Can you please guide me which steps I should study.
Which step should be considered as the base of my work and by modifying it I can model a thermoelastic problem?

Best  regards
Lisa

[deal.II-MailingList]

Hello Lisa,

I suppose that it depends on how you want to solve your thermoelastic problem.

Do you want to assemble one big system which contains one block for the elastic equation,

another block for the heat equation and coupling blocks? Then you could solve that system

in a monolithic way. Therefore you should have a look at step 20, 21 and 22.

If you want to solve the elastic equation and the heat equation alternately a look to step 18

and 26 might be enough.

step 18 shows an example for how to solve an elastic equation as a time depend problem.

And step 26 gives an example for how to deal with the heat equation.

Since I am currently working on a thermoelastic problem let me know if you have further

questions.

Best,

Joerg

[Sudharsan Parthasarathy]

Hi,

   I am working on a similar problem and have a follow-up question. If one were to use the second approach using step-18 and step-26, how should one go about using the temperature field obtained from the heat equation in solving the elastic problem? Can VectorTools::point_value() be used along with the DOFHandler object of the TopLevel class(of step-18) or should one somehow share the DOFHandler object of the HeatEquation class in order to use the point_value function? Is there any other way of doing it?

Thanks,

Sudharsan

[Daniel Arndt]

Hello Sudharsan,

what you want to do is to step through all the cells and reinitialize both DoFHandlers simultaneously, i.e. something like

cell_elastic=dh_elastic.begin_active();

cell_temperature=dh_temperature.begin_active();

for (;cell_elastic != dh_elastic.end(); ++cell_elastic, ++cell_temperature)

{

  fe_values_elastic.reinit(cell_elastic);

  fe_values_temperature.reinit(cell_temperature);

  fe_values_temperature.get_function_values(temperature, local_temperature);

...

}

This of course can only work if both DoFHandlers are based on the same Triangulation.

Best regards,

Daniel

[Sudharsan Parthasarathy]

Hi Daniel,

   Thanks for the reply. 

    I used to have my code in the same way as you have stated in your reply. However, this approach requires the DOFHandler object of both the thermal and elastic problems to be shared and I want to avoid such an approach. Can it be avoided? Is there another approach? For example, will I be able to use VectorTools::point_value() using the temperature Vector as the fe_function given in the documentation and the dh_elastic as the DOFHandler given in the documentation?

Thanks,

Sudharsan

[Daniel Arndt]

Sudharsan,

I am not entirely clear, what you mean by "this approach requires the DoFHandler object of both the thermal and elastic problems to be shared". In what I was suggesting you have two separate DoFHandler objects that only share a common triangulation.

Of course, you can also use VectorTools::point_value for evaluating the temperature Vector. This however means that for each point the corresponding cell is searched for and the a FEValues object is created and initialized on that specific cell. Therefore, this approach is much more expensive.

Evaluating a FEFunction with a non corresponding DoFHandler simply doesn't make sense.

Best,

Daniel

[hamedb...@gmail.com]

Hello Joerg,

I want to deal with a thermoelastic problem in a way that the elastic equation and the heat equation are solved alternately.
I was wondering if I can define a FESystem<dim> and Extractors as follows :

fe (FE_Q<dim>(degree+1), dim,
FE_Q<dim>(degree), 1),

const FEValuesExtractors::Vector displacement (0);

const FEValuesExtractors::Scalar temprature (dim);

for displacement and temperature fields as done in step 20 or 22 for mixed vector-valued problems but don't solve problem in a monolithic way.
I am struggling to combine step 18 and 26 with step 20 to define one mixed FESystem and DOFHandler for elastic and heat equation. For example in step 18 there is a function to compute the strain as follows:

template <int dim>

inline

SymmetricTensor<2,dim>

get_strain (const FEValues<dim> &fe_values,

const unsigned int shape_func,

const unsigned int q_point)

{

SymmetricTensor<2,dim> tmp;

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

tmp[i][i] = fe_values.shape_grad_component (shape_func,q_point,i)[i];

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

for (unsigned int j=i+1; j<dim; ++j)

tmp[i][j]

= (fe_values.shape_grad_component (shape_func,q_point,i)[j] +

fe_values.shape_grad_component (shape_func,q_point,j)[i]) / 2;

return tmp;

}

and I don't know if I can change the last lines to something like this:
...

= (fe_values[displacement].shape_grad_component (shape_func,q_point,i)[j] +

fe_values[displacement].shape_grad_component (shape_func,q_point,j)[i]) / 2;

Thanks for your time,
Hamed

[Daniel Arndt]

Hamed,

There shouldn't be a problem to do what you following step-22. In particular, you can modify step-18::get_strain as you proposed.

Did you face any difficulties so far?

Best,

Daniel

[hamedb...@gmail.com]

Daniel,

Since in the FEValuesViews::Vector< dim, spacedim > Class there isn't any shape_grad_component function, I am not sure that I can update the get_strain as I mentioned.
I was wondering if you have already  written a code for a thermoelastic problem in dealii, solving heat equation and elastic equation separately.

Best,
Hamed

[Daniel Arndt]

You're right, Hamed. There is FEValuesViews::Vector< dim, spacedim >::symmetric_gradient [1]
and FEValuesViews::Vector< dim, spacedim >::gradient [2] instead. The first should be exactly what you are looking for.

Code that I have written normally uses multiple DoFHandlers if I solve for multiple equations separately.

Best,
Daniel

[1] https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#a4e5dfbb49d284a368acac6ef5dd4b2ef
[2] https://dealii.org/8.4.1/doxygen/deal.II/classFEValuesViews_1_1Vector.html#ad7b4df64147989229f566eff541ddebc

[Hamed Babaei]

Hello Daniel,

Since although we can solve two elastic and heat equations separately, we still need to enter stress computed from elastic equation into the heat equation in every time step and keep updating stress field.
I was wondering when we have two different DoFHandlers for each of the equations, how to use, for instance, stress field stored on a quadrature point related to elastic problem DofHandler, in order to solve stress-induced heat equation problem. In other words, I don't know how to exchange information between two equation.

Thanks,
Hamed

[Daniel Arndt]

Hamed,

assuming that both DoFHandlers are based on the same triangulation, you can do something like the following:

FEValues fe_stress_values(...);
stress_cell = dof_handler_stress.begin_active();
temperature_cell = dof_handler_temperature.begin_active();
for (; stress_cell != dof_handler_stress.end(); ++stress_cell, ++temperature_cell)
{
  fe_stress_values.reinit(stress_cell);
  fe_stress_values.get_function_values(solution_stress, local_stress);
  //do all the temperature related stuff and use local_stress
}

Instead of having two iterators separated, you can also use SynchronousIterators [1]. They arec partly used in step-35 for that purpose.

Best,
Daniel

[1] https://dealii.org/8.4.1/doxygen/deal.II/structSynchronousIterators.html

[ha.ba...@gmail.com]

I have solve a thermoelastic problem using block matrices and block vectors and a same dofhandler for the displacement and temperature fields.

Fully coupled thermal-stress analysis is needed when the stress analysis is dependent on the temperature distribution and the temperature distribution depends on the stress solution. For example, metalworking problems may include significant heating due to inelastic deformation of the material which, in turn, changes the material properties. Some problems require a fully coupled analysis in the sense that the mechanical and thermal solutions evolve simultaneously, but with a weak coupling between the two solutions. In other words, the components in the off-diagonal submatrices are small compared to the components in the diagonal submatrices. For the thermoelastic problem the off-diagonal terms are zero and solution is straight and stable. You can solve the block systems simultaneously or using staggered solution technique.  

Setup system as follow:

dof_handler.distribute_dofs (fe);

DoFRenumbering::Cuthill_McKee(dof_handler);

BlockDynamicSparsityPattern dsp(2, 2);

std::vector<unsigned int> block_component(3, 0);

block_component[dim] = 1; // temperature

DoFRenumbering::component_wise(dof_handler, block_component);

DoFTools::count_dofs_per_block(dof_handler, dofs_per_block,block_component);

const types::global_dof_index n_dofs_u = dofs_per_block[0];

// determine number of temperature DOFs if requested.

types::global_dof_index n_dofs_temp;

n_dofs_temp = dofs_per_block[1];

dsp.block(0, 0).reinit(n_dofs_u, n_dofs_u);

dsp.block(0, 1).reinit(n_dofs_u, n_dofs_temp);

dsp.block(1, 0).reinit(n_dofs_temp, n_dofs_u);

dsp.block(1, 1).reinit(n_dofs_temp, n_dofs_temp);

dsp.collect_sizes();

...

use the following approach for system matrix and RHS:

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

{

// extract temperature gradient and value

for (unsigned int k = 0; k < dofs_per_cell; ++k){

temp_grads_N[k] = fe_values[temp_fe].gradient (k, q);

temp_N[k] = fe_values[temp_fe].value (k, q);

}

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

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

{

const unsigned int comp_j = fe.system_to_component_index(j).first;

if (comp_j < dim){

cell_matrix(i,j) += determine block(0,0)

}else if (comp_j == dim){

// temperature

cell_matrix(i,j) += determine block(1,1)

}

}

}

cell->get_dof_indices (local_dof_indices);

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

{

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

system_matrix.add (local_dof_indices[i],local_dof_indices[j],

cell_matrix(i,j));

}

for rhs:

loop over cells:

PointHistory<dim> *local_quadrature_points_data

= reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());

// extract nodal total and incremental displacement

// will be used for strain and delta_strain calculation

fe_values[u_fe].get_function_gradients (solution_delta, displacement_increment_grads);

fe_values[u_fe].get_function_gradients (solution, displacement_grads);

// extract temperature values and gradients

if(parameters.is_temperature){

fe_values[temp_fe].get_function_values (solution, temp_value);

fe_values[temp_fe].get_function_values (solution_n, temp_n_value);

fe_values[temp_fe].get_function_gradients (solution, temp_grads);

}

// Then we can loop over all degrees of freedom on this cell and

// compute local contributions to the right hand side:

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

{

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

const unsigned int comp_i = fe.system_to_component_index(i).first;

if ((comp_i < dim) ){

 cell_rhs(i) +=

}else if (comp_i == dim)

{

// temperature

cell_rhs(i) +=

}

}

}

cell->get_dof_indices (local_dof_indices);

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

system_rhs(local_dof_indices[i]) += cell_rhs(i);

[Hamed Babaei]

Thank you for your helpful information.

I was wondering if there would be any difficulties for parallelization of the code when using block matrices and vectors. would you please let me know if you have paralleled your code.

Thanks,
Hamed

[ha.ba...@gmail.com]

Hello, dont worry about extension to parallel solution. I have used dealii::LinearAlgebraTrilinos. Assembly process is exactly similar. The only change is cell->is_locally_owned() at the start of the loop over cells. 

[Muhammad Mashhood]

Hi Joerg and Lisa!

Thank you for informative reply and posting this concern on the forum. I am also interested in thermoelastic problem and new use of deal.ii.

My question is that other than the tutorial steps 26 & 18 or 20,21 & 22, is there pre-developed application at "https://www.dealii.org/" for this field of study? Thank you!  

[Wolfgang Bangerth]

On 5/24/19 2:46 AM, Muhammad Mashhood wrote:
> Thank you for informative reply and posting this concern on the forum. I am
> also interested in thermoelastic problem and new use of deal.ii.
> My question is that other than the tutorial steps 26 & 18 or 20,21 & 22, is
> there pre-developed application at "https://www.dealii.org/" for this field of
> study?

None that immediately implements the thermoelastic problem. But there are of
course many building blocks you can find in a variety of tutorial programs and
code gallery programs.

Best
W.

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

[Muhammad Mashhood]

Hi Daniel,

               Thank you for the suggestion and I also think it can be one of the way to successfully run both solutions by sharing information on correct corresponding q_points.

I tried the same approach (as I am also using partitioned approach for thermomechanical coupling) but so far my (mesh files being read for both thermal and solid mechanics parts are same but ) triangulation are different.

So I also want to use the same triangulation with different dof_handler for both parts.

The scenario of workflow in my case is following:

The thermal code uses the Triangulation<dim> triangulation_thermal while the solid mechanics code is using parallel::distributed::Triangulation<dim> triangulation_solid. Yes the thermal part is programmed in serial and the solid mechanics part in parallel. The solid mechanics code runs and makes the grid with it's triangulation, As I want to use the same triangulation so I try to copy this triangulation with function triangulation_solid.copy_triangulation(triangulation_thermal) and here the triangulation_thermal is the one which I will pass to the name space and class of thermal analysis part before running it. Moreover during my analysis of solid mechanics part I also refine the mesh (and would also want to use the same refined mesh for the thermal analysis too by copying it). 

After describing the scenario following are my concerns where I need your expert opinion :

1) when I use the copy_triangulation() function, I come across the error (even here the mesh is not refined so far but still) mentioning :

  void dealii::Triangulation<dim, spacedim>::copy_triangulation(const dealii::Triangulation<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
The violated condition was:
    (vertices.size() == 0) && (levels.size() == 0) && (faces == nullptr)
Additional information:
    You are trying to perform an operation on a triangulation that is only allowed if the triangulation is currently empty. However, it currently stores 100 vertices and has cells on 1 levels.

2) As far as the documentation of tria.h is concerned it also informs that the copy_triangulation() cannot copy the refined mesh (where in my case the mesh might also get refined during analysis based on certain criteria and I would like to use same refined triangulation by copying it for thermal analysis so that triangulation is same in both problems).

considering the above scenario as well as the concerns, I would be grateful to receive any suggestion from your side. Hope I am clear in my description.

Waiting for your kind response. Thank you in advance!

Best regards,

Muhammad

[Wolfgang Bangerth]

On 8/2/19 4:01 AM, Muhammad Mashhood wrote:
>
> considering the above scenario as well as the concerns, I would be grateful to
> receive any suggestion from your side. Hope I am clear in my description.
> Waiting for your kind response. Thank you in advance!

Instead of copying triangulations, you always have the option of just creating
the two objects the same, and then refining them in exactly the same way.
Would that solve your problem?

[Muhammad Mashhood]

Dear Prof. Wolfgang,

                                  Thank you for your concern and understanding.

Yes you are right. I tried something similar instead of copying the triangulation.

I used const SmartPointer<const Triangulation<dim> > to transfer triangulation between both thermal and the solid mechanics classes and seemingly the objective of using same triangulation for both parts is full filled also with same refinement and coarsening.

While for the sake of using different dof_handler for temperature and solid mechanics, I currently have to define again the temperature_dof_handler as well as the fe_temperature objects in the solid mechanics (class) part to do something like following:

          dof_handler_temperature.distribute_dofs(fe_temperature);

          cell_solid_mech = dof_handler.begin_active();

          cell_temperature = dof_handler_temperature.begin_active();

(for loop for cell_solid_mech and cell_temperature )

          fe_values_solid_mech.reinit(cell_solid_mech);
          fe_values_temperature.reinit(cell_temperature);
          cell_matrix = 0;
          cell_rhs = 0;
          fe_values_temperature.get_function_values(temperature_solution, temperature_solution_qpoint);
          fe_values_solid_mech[displacement].get_function_symmetric_gradients(displacement_solution,
                                                                   strain_tensor); 

So in this way, I am successful in using the common triangulation for both parts but with different corresponding dof_handlers.

The only thing is that I was wondering if I might be using extra processing and memory by defining the fe_temperature, dof_handler_temperature and fe_values_temperature by defining them again in the solid mechanics program. And if I am, then would there be any alternative method where I can use the const SmartPointer<const DoFHandler<dim> > and const SmartPointer<const FE_Q<dim><dim> > or something similar to return  fe_temperature and dof_handler_temperature from HeatEquation class to SolidMechanics class (where they are being used for evaluating q_point temperature as I mentioned above) ? Thank you!

[Wolfgang Bangerth]

On 8/5/19 3:40 AM, Muhammad Mashhood wrote:
> The only thing is that I was wondering if I might be using extra processing
> and memory by defining the fe_temperature, dof_handler_temperature and
> fe_values_temperature by defining them again in the solid mechanics program.
> And if I am, then would there be any alternative method where I can use the

> *const SmartPointer<const DoFHandler<dim> > *and *const SmartPointer<const
> FE_Q<dim><dim> > *or something similar**to return *fe_temperature* and
> *dof_handler_temperature* from *HeatEquation* class to *SolidMechanics* class

> (where they are being used for evaluating q_point temperature as I mentioned
> above) ? Thank you!

Muhammad,
you are in essence asking general C++ programming questions about the
difference between regular member variables and pointers, and how to return
objects by value, by reference, or by pointer. The answer is that yes, you can
do as you suggest above. But I would suggest you take a look at a good C++
book or online resource to understand *why* the answer is yes :-)

[Muhammad Mashhood]

Prof. Wolfgang,

                         Thank you for the concern and response. I got your point there. I will try it in this way then. :)