I am attempting to solve for a displacement field where I know the thermal stresses on my discretized domain, which consists of hex cells.

My first question is: is easier way to solve for the displacements besides integrating the linear elasticity equations: $$ \nabla \cdot \sigma = 0 $$ where $\sigma$ is stress (contains both the mechanical and thermal stresses).

I currently am using the finite volume method to solve the above equation, and the integral form is $$ \int \int (\sigma \cdot \hat{n}) dA = 0 \ \ \ \ \text{ applied divergence theorem}\\ \int \int (\sigma \cdot \hat{n}) dA \approx \sum_{f=1}^N \sigma_fA_f\hat{n}_f $$ where the integral is approximated using the midpoint rule that reads as the sum of the product of the face stress, face area, and face normal over all faces of a cell. The face stress is computed by taking an average of the neighboring cell's stress.

The stress uses the duhamel-neumann relationship for an anisotropic material (note that I am currently testing the code on a simple isotropic material) $$ \sigma = C:(\epsilon_m - \epsilon_t) $$ Where $\epsilon_m$ is the mechanical strain that is computed using the linear strain displacement relationship and $\epsilon_t$ is the thermal strain, computed using $$ \epsilon_t = \alpha (T-T_0) $$ where $\alpha$ is the CTE.

I can formulate the above system of equations as a linear system $Au=b$, where $u$ are the displacements to be solved for. $b$ is the forcing and essentially all the thermal stress term is contained inside $b$ (thermal stress has no impact on $A$).

The good news is that if I assign the same thermal stress to the entire domain (i.e., same temperature at all the cells), then I get a solution that is "correct." By correct, I mean it matches commercial software (ANSYS). So something is correct.

However, if my temperature distribution across the domain is not the same, then I get solutions that not only do not come close to matching ANSYS, but also do not match my intuition (certain cells with larger temperature changes are expanding less or even contracting(!), as indicated by the computed displacement field, than parts with smaller or no temperature changes.

I have checked my code and it seems I am implementing everything correctly, so I am wondering if my numerics are incorrect. But I am getting the correct solution if the whole field is initialized to the same temperature, so something must be correct. e.g., does adding the thermal stress term somehow affect how I should approach the discretization? I would think not. This term is easy and it's really just a forcing term..

As previously mentioned, the thermal stress is manifested mathematically in the forcing vector $b$ and has no impact on $A$. So changing the temperature just affects the magnitude of my forcing vector, and my $A$ remains the same. So this tells me that my $A$ matrix is formulated correctly, or that it could be formulated incorrectly, but in a way that has not affected my solutions for other problems. I think the latter is very unlikely because I have also tested my code on cantilever and various other bending/compression problems (with NO thermal stress), and I get solutions that match well with analytical solutions or with commercial software. I have also done some manufactured solutions studies to look at convergence and the tests succeeded.

It's just this thermal stress term that's causing me issues, and only when the temperature is not uniform across the discretized domain.

Any suggestions? I'm happy to provide more info if needed.

Edit #2: Here is the particular geometry I am working with right now: Geometry Mesh

The length is in the x direction. I apply a uniform temperature of $1900K$ to all the cells. The following picture is the displacement contours for the y-direction, which is the more refined direction in the above picture. enter image description here

The max y-displacement from the contours for the +y and -y surfaces are about 0.0006m and -0.0006m. This material is steel, using CTE = $1.2\times 10^{-5} K^{-1}$. The dimension of the block is about 0.6098m by 0.0047625m by 0.00508m (sorry multiple decimals. The original grid was created in inches) in the x,y,z directions. The mesh step size in the 3 directions is constant, at 0.03048m (1.2inches) by 0.004318m (0.17 inches) by 0.0254m (1 inch).

So here I get a solution that matches what ANSYS gives me.

Next I proceeded to run a different simulation where I only specify a temperature of 1900K in the domain in the top two layers (see the region in the white boxes below) and the rest of the domain is at 298K. enter image description here

In this simulation where only part of the domain experiences a temperature change, I get a solution that is way off from expected or from what ANSYS gives me.

Edit #3: Adding displacement contours for the aforementioned test case.

Here are the 3 displacement contours for the x-,y-, and z-displacements for the second case. My intuition tells me that the x- and z- displacement contours show trends that are expected (not sure if the magnitudes are correct, but at least the trends are reasonable). But I would not expect the y-displacement to behave the way it is shown in the y-displacement contours.

Z-displacement contour

X-displacement contours

Y-Displacement contours

  • 1
    $\begingroup$ In your equation for $\sigma$, I think that plus sign should be a minus. $\endgroup$ Nov 3 '18 at 17:53
  • $\begingroup$ Yes, sorry about that. I changed it to the minus now. I have it correctly in my code. $\endgroup$
    – roulette01
    Nov 3 '18 at 18:02
  • $\begingroup$ So, just to make sure I understand your situation, if you take a block of material and, say, heat it uniformly, the displacements you calculate are basically $\alpha \Delta T$ times the length dimension? $\endgroup$ Nov 3 '18 at 18:42
  • 1
    $\begingroup$ Have you carried out verification of your code using something like the method of manufactured solution ? It would be an easy place to start $\endgroup$
    – BlaB
    Nov 3 '18 at 18:42
  • 1
    $\begingroup$ Another easy test would be to set Poisson's ratio equal zero, apply a temperature that varies linearly in one direction and see if your code matches the easily-calculated values; 1D solutions are particularly easy to solve by hand. $\endgroup$ Nov 3 '18 at 18:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.