import theme
theme.load_style()
This lecture by Tim Fuller is licensed under the Creative Commons Attribution 4.0 International License. All code examples are also licensed under the MIT license.
Using simple axial bar elements and transformation equations we were able to solve for the deformation of a 3D space truss under prescribed load/support conditions.
Using the same connectivity table we were able to solve for the temperature distribution assuming a prescribed temperature at the support and loaded nodes, and neglecting the effect of convection to the surrounds.
In reality, the heat loss to the surrounding air would be significant and should be considered in the design process.
Thermal expansion and self weight may significantly affect the structural response.
For the bar shown, two heat transfer mechanisms occur at each point along the bar
Surface convection: The area specific heat to the surroundings is equal to the product of the heat transfer coefficient and the difference between the surface temperature and free stream temperature.
$$ q_{\text{conv}} = h\left(T_x - T_{\infty}\right) $$For free convection to air, typical values for are 2 - 25 w/m$^2$ K and can be estimated based on model geometry and gas properties.
Average heat transfer coefficient
$$ \overline{h} = \frac{\overline{Nu_D}k}{D} $$Nusselt number
$$ \overline{Nu_D} = \left(.6+ \frac{.387Ra_D^{1/6}}{\left(1+\left(.599/Pr\right)^{9/16}\right)^{8/27}}\right)^2 $$Rayleigh number:
$$ Ra_D = \frac{g\beta\left(T_s-T_{\infty}\right)D^3}{\nu\alpha} $$Ideal gas compressibility:
$$ \beta = \frac{1}{\left(T_s-T_{\infty}\right)/2} $$Evaluating this with properties for air at 300K, and allowing that the surface temperature may vary with position along the bar, we obtain the following:
$$ h(x) = \frac{.026}{d}\left(19.8\left(\lvert T(x)-T_{\infty}\rvert \right)^{1/6} \left(\frac{d^3}{T(x)+T_{\infty}}\right)^{1/6}+.6\right)^2 $$At each point, the axial conduction is equal to the conduction out plus the convection to the surroundings.
For some differential length:
$$ q(x)A(x) - \beta h\left(T(x) - T_\infty(x)\right)\Delta x - q(x+\Delta x)A\left(x+\Delta x\right)=0 $$Taking the limit as $\Delta x\rightarrow 0$, we obtain
$$ \beta h\left(T(x) - T_{\infty}\right)+A(x)\frac{dq}{dx}=0 $$For heat conduction, the flux-potential relationship is modeled by Fourier's law:
$$ T(x) - T(x+\Delta x) = q(x)\frac{\Delta x}{k_{th}A} $$or, in differential form
$$ q(x) = -k_{th}\frac{dT}{dx} $$Substituting the flux potential relationship into the balance law we obtain the strong form for axial conduction in a bar with surface convection
$$ k_{th}A(x)\frac{d^2T}{dx^2}-\beta h T\left(T(x)-T_{\infty}\right)=0, \quad 0 \leq x \leq L $$with Dirichlet boundary conditions
$$ T(0) = T_h, \quad T(L) = T_c $$Assuming the problem data (what are all the problem data?) are constant, the strong form is
$$ kAT''-\beta h\left(T-T_{\infty}\right)=0, \quad 0 \leq x \leq L $$$$ T(0) = T_h, \quad T(L) = T_c $$Following the 3 step procedure
So, the weak form for axial conduction in a bar with surface convection and Dirichlet boundaries is
$$ \int_0^L w'kAT'dx + \int_0^L\beta hTdx - \int_0^L\beta hT_{\infty}dx=0, \quad \forall w, \ w(0) = w(L) = 0 $$Express the integral over the domain as a sum of integrals over elements
$$ \begin{multline} \sum_e \left( \int_{x_0^{(e)}}^{x_L^{(e)}} {w^{(e)}}' k^{(e)}A^{(e)} {T^{(e)}}' dx + \int_{x_0^{(e)}}^{x_L^{(e)}} w^{(e)} \beta^{(e)} h^{(e)} T^{(e)} dx - \\ \int_{x_0^{(e)}}^{x_L^{(e)}} w^{(e)} \beta^{(e)} h^{(e)} T_{\infty} dx \right)=0, \quad \forall w, \ w(0) = w(L) = 0 \end{multline} $$Introduce an approximation for the "trial" solution and weight "test" function
$$ T^e = N_i^e T_i^e, \quad \frac{dT^e}{dx} = \frac{dN_i^e}{dx}T_i^e = B_i^eT_i^e, \quad w^e = w_i^eN_i^e, \quad \frac{dw^e}{dx} = w_i^e\frac{dN_i^e}{dx} $$Substituting in to the weak form gives
$$ \begin{multline} \sum_e\left(\int_{x_0^{(e)}}^{x_L^{(e)}} w_i^{(e)}\frac{dN_i^{(e)}}{dx}k^{(e)}A^{(e)}\frac{dN_j^{(e)}}{dx} T_j^{(e)} + \int_{x_0^{(e)}}^{x_L^{(e)}} w_i^{(e)}N_i^{(e)}\beta^{(e)} h^{(e)} N_j^{(e)}T_j^{(e)}dx \right. \\ \left. - \int_{x_0^{(e)}}^{x_L^{(e)}} w_i^{(e)}N_i^{(e)} \beta^{(e)} h^{(e)} T_{\infty}dx\right)=0, \quad \forall w, \ w_1 = w_n = 0 \end{multline} $$We factor out the $w_i^{(e)}$ of each term and rearrange
$$ \begin{multline} \sum_ew_i^{(e)}\left[ \left( \int_{x_0^{(e)}}^{x_L^{(e)}} k^{(e)}A^{(e)}\frac{dN_i^{(e)}}{dx}\frac{dN_j^{(e)}}{dx} + \int_{x_0^{(e)}}^{x_L^{(e)}} \beta^{(e)} h^{(e)} N_i^{(e)} N_j^{(e)}dx\right)T_j^{(e)} \right.\\ \left.- \int_{x_0^{(e)}}^{x_L^{(e)}} \beta^{(e)} h^{(e)} N_i^{(e)} T_{\infty}dx\right] = 0, \quad \forall w, \ w_1 = w_n = 0 \end{multline} $$Writing this in a more compact form
$$ \sum_e w_i^{(e)} \left(k_{ij}^{(e)}T_j^{(e)} - f_i^{(e)}\right) = 0 $$where the stiffness $k_{ij}^{(e)}$ is
$$ k_{ij}^{(e)}=\int_{x_0^{(e)}}^{x_L^{(e)}} k^{(e)}A^{(e)}\frac{dN_i^{(e)}}{dx}\frac{dN_j^{(e)}}{dx} + \int_{x_0^{(e)}}^{x_L^{(e)}} \beta^{(e)} h^{(e)} N_i^{(e)} N_j^{(e)}dx $$and the element external flux array, (equivalent to the force array for our structural problem), which is a sum of boundary and body forces
$$ f_i^{(e)} = f_{\Gamma i}^{(e)} + f_{\Omega i}^{(e)} = \int_{x_0^{(e)}}^{x_L^{(e)}} \beta^{(e)} h^{(e)} N_i^{(e)} T_{\infty}dx $$For this problem the boundary source array is zero, as we have specified Dirichlet boundaries at each end of the domain. This is often done in discussion of FE when the focus is the development of an element stiffness matrix. We have not yet specified the function $T_{\infty}$ that determines how the free stream temperature varies along the bar.
From the temperature distribution in the truss members, we can compute the thermal expansion strains to see how these affect the stresses and displacements in the structure.
Strong Form
$$ \frac{d}{dx}\left(EA(x)\frac{du}{dx}-E\alpha A(x) T(x)\right) = 0 $$Weak Form
$$ \int_0^Lw'(x) EA u'(x) dx + \int_0^L w(x) EA\alpha T'(x) dx = 0, \quad \forall w(x), w(x)_{x=0,L}=0 $$Finite Element Equations
$$ k_{ij}^{(e)}u_j^{(e)} - f_i^{(e)} = 0 $$where
$$ k_{ij}^{(e)} = \int_{x_0^{(e)}}^{x_L^{(e)}} A^{(e)}E^{(e)}\frac{dN_i}{dx}\frac{dN_j}{dx} dx $$$$ f_i^{(e)} = -\int_{x_0^{(e)}}^{x_L^{(e)}}N_i^{(e)}(x)EA\alpha\frac{dT^{(e)}}{dx}dx $$The thermal stresses appear in the form of a body force, adding a force that counteracts that normally required to achieve some amount of thermal expansion.
Modified from Fish & Belytschko, Problem 5.1 For this problem do NOT use the direct stiffness method.
Consider a heat conduction problem in the domain [0, 20]m. The bar has a unit cross section, constant thermal conductivity $k=5$W/$^{\text{o}}C$$\cdot$m and a uniform heat source $s=100$ W/m. The boundary conditions are $T(0)=0^{\text{o}}C$ and $\overline{q}(20)=0$W/m$^2$. Solve the problem with two equal linear elements. Plot the finite element solution $T_N(x)$ and $dT_N(x)/dx$ and compare to the exact solution which is given by $T(x)=-10x^2+400x$.
For the bar shown,
a) State the strong form representing heat flow and solve analytically. Give the symbolic solution $T(x)$ and total heat flux $Q(x)$ along the length of the bar.
b) Construct the element body source arrays, the boundary force arrays, and assemble the global external force array.
c) Construct the element stiffness matrices and assemble the global stiffness matrix.
d) Solve for the temperature distribution $T(x)$. Generate a plot of the global approximate solution $T_N(x)$ along the length of the bar overlaid with the analytic solution.