Coupled thermo elasticity

[pcarruscag]

Continuation of #2398 and #2399 for @Vaish-W to work on.

The idea is to instantiate the heat solver alongside the FEA solver and then access the temperature at the nodes when we compute stresses, stiffness matrix, etc.

The instantiation of solvers happens in CSolverFactory::CreateSolverContainer, so look for places with solver[FEA_SOL] and add the creation of a heat solver conditional on some option, we can re-use WEAKLY_COUPLED_HEAT_EQUATION.

This should automatically create a SingleGridIntegration class for the heat problem together with the StructuralIntegration for the normal FEA problem.

Now we need to modify the CFEAIteration class to also integrate and advance the heat problem. This is in CFEAIteration::Iterate, before the /*--- FEA equations ---*/ section we can add a call to the SingleGrid_Iteration of the heat integration class, see CFluidIteration.cpp around line 117 for what you need to do.

At this point we should have temperatures being computed by the heat solver and we can start using them in the FEA solver.
We need to do this in all functions that compute stresses Compute_StiffMatrix, Compute_StiffMatrix_NodalStressRes, Compute_NodalStressRes, and Compute_NodalStress.
Where these functions do element->SetRef_Coord we need to also SetTemperature (if it is available). These functions do not have access to heat solver at the moment, but CFEASolver::Preprocessing sees the entire solver container, so here we can store a pointer to the heat_solver or it's solution.

The last thing to do is to change CElasticityOutput to also output residual information and solution for the heat solver. You can look at CHeatOutput.cpp to see what is available in the heat solver.

Let me know if this makes sense to you and open a draft PR once you have some code so I can comment early on.
We also have weekly developer meetings we can use to talk about the project, the schedule and instructions are in our Slack channel.

[Vaish-W]

Thank you @pcarruscag for your detailed feedback on the pull request. I will surely work on this soon!!

[Vaish-W]

Hello,
I did first two steps as you suggested in PR #2404 but now stuck at preprocessing update.

I have:

1. Introduced a CHeatSolver* heat_solver member to the CFEASolver class.
2. initialized it to nullptr in the constructor and deallocated it in the destructor.

Now moving to Preprocessor:
What would be approprate logic to preprocessing function which includes initializing the heat solver when it is enabled from config?

[pcarruscag]

Hi,
Make it just a CSolver* to make it simpler, then in preprocessing of the FEASolver you set this pointer to solver to the heat solver that was already created
heat_solver = solver_container[...][HEAT_SOL];
You don't need to deallocate because this is just a pointer to a solver that is managed (created and destroyed) outside of the FEA solver.

[Vaish-W]

Hey, how can I test the code we completed other than regression and unit tests? CI pipeline passed but should I prepare a problem for thermal stress simulation of a simple beam or plate?

[pcarruscag]

Yes, I would say the first verification step should be with a setup for the heat solver that gives an isothermal solution and comparing this with just setting the free stream temperature.
Then I would check an unloaded beam with linear temperature distribution, which has an easy analytical solution for the deformation due to thermal expansion.

[Vaish-W]

Can I use cantilever beam config_linear.cfg tutorial setup for this? And do these modifications:
Enable WEAKLY_COUPLED_HEAT_EQUATION.
Set all boundaries to isothermal (e.g., MARKER_ISOTHERMAL = (clamped, 300.0)).
Use FREESTREAM_TEMPERATURE = 300.0 for comparison

is it correct to do so?

[pcarruscag]

Yes, but in the case with weakly coupled heat set the free stream temperature to a different value from the isothermal value to make sure it's not used by mistake.

[Vaish-W]

thermoelasticity.txt
I am getting following output where it gives 0 von mises stress for this cfg.

-------------------- Solver Preprocessing ( Zone 0 ) --------------------
Filename: default_element_properties.dat.
There is no element-based properties file.
The structural domain has uniform properties.
Initialize Jacobian structure (Non-Linear Elasticity).
Initialize Jacobian structure (heat equation) MG level: 0.

------------------- Numerics Preprocessing ( Zone 0 ) -------------------

----------------- Integration Preprocessing ( Zone 0 ) ------------------

------------------- Iteration Preprocessing ( Zone 0 ) ------------------
FEM iteration.

------------------------------ Begin Solver -----------------------------

Simulation Run using the Single-zone Driver
CSysSolve::ConjugateGradient(): system solved by initial guess.
CSysSolve::ConjugateGradient(): system solved by initial guess.
+----------------------------------------------------------------+
|  Inner_Iter|  rms[DispX]|  rms[DispY]|      rms[T]|    VonMises|
+----------------------------------------------------------------+
|           0|  -32.000000|  -32.000000|  -32.000000|  0.0000e+00|

----------------------------- Solver Exit -------------------------------
All convergence criteria satisfied.
+-----------------------------------------------------------------------+
|      Convergence Field     |     Value    |   Criterion  |  Converged |
+-----------------------------------------------------------------------+
|                  rms[DispX]|           -32|          < -8|         Yes|
+-----------------------------------------------------------------------+
-------------------------------------------------------------------------
+-----------------------------------------------------------------------+
|        File Writing Summary       |              Filename             |
+-----------------------------------------------------------------------+
|SU2 binary restart                 |restart_linear_elastic.dat         |
|Paraview                           |linear_elastic.vtu                 |
+-----------------------------------------------------------------------+

------------------------- Exit Success (SU2_CFD) ------------------------

I have applied boundary condition for thermal as well as fea. What mistake did I make in cfg or could that be related to solver problem?

[Vaish-W]

If I remove left and clamped boundary condition from MARKER_ISOTHERMAL then I am getting stress as 1.1541e+06. Is that because fea and thermal boundary conditions are conflicting with each other and hence not accepting both condition for one marker?

+----------------------------------------------------------------+
| Inner_Iter| rms[DispX]| rms[DispY]| rms[T]| VonMises|
+----------------------------------------------------------------+
| 0| -9.644128| -9.726810| -32.000000| 1.1541e+06|

[pcarruscag]

Please reopen the pull request with the new branch

[pcarruscag]

Don't remove the markers, just make the load 0.

[Vaish-W]

Still getting same output with no convergence

[pcarruscag]

Well, this is what it's like to develop CFD codes :) you just have to debug things...

[Vaish-W]

Yes, I was just trying to understand if my .cfg file contains any errors before debugging the code... But will figure it out surely :)

[pcarruscag]

Some.
MARKER_ISOTHERMAL = (left_thermal, 350, right_thermal, 300, right, 0, left, 0, upper, 0)
You should not put right, left, etc. as 0 temperature, they should be 0 heat flux.