diff --git a/docs/febioheat_theory.md b/docs/febioheat_theory.md new file mode 100644 index 0000000..188fafc --- /dev/null +++ b/docs/febioheat_theory.md @@ -0,0 +1,126 @@ +## Introduction +Theory Manual +The heat transfer equation is defined as follows. + +\[ +\begin{equation} +\label{eq:heat-transfer} +\rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q}=Q +\end{equation} +\] + +Here, $u$ is the _temperature_, $\rho$ the _density_, $C$ the _specific heat_, $\mathbf{q}$ the _heat flux_, and $Q$ the _heat source_. + +The following boundary conditions are assumed. + +\[\begin{align} + & u={{u}_{0}}\quad\text{on}\,{{\Gamma }_{D}} \\ + & \mathbf{q}\cdot \mathbf{n}={{q}_{0}}\quad\text{on }\,{{\Gamma }_{N}} \\ +\end{align}\] + + +Furthermore, it is assumed that the heat flux is defined via Fourier’s law, + +\[ +\begin{equation} +\label{eq:heat-flux} +\mathbf{q}=-k\,\nabla u +\end{equation} +\] + +where $k$ is the thermal conductivity. + +## Weak Formulation +To construct the weak form, first multiply the left- and right-hand side with an arbitrary function $w$ (that satisfies the homogenous boundary conditions, i.e. $w=0$ on ${{\Gamma }_{D}}$). + +\[w\left( \rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q} \right)=wQ\] + +Then, integrate over the problem domain. + +\[\int\limits_{\Omega }{w\left( \rho C\frac{\partial u}{\partial t}+\nabla \cdot \mathbf{q} \right)d\Omega }=\int\limits_{\Omega }{wQd\Omega }\] + +Using integration by parts and the boundary conditions, we can rewrite the second term on the left-hand side. + +\[\int\limits_{\Omega }{w\rho C\frac{\partial u}{\partial t}d\Omega }+\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }-\int\limits_{\Omega }{\nabla w\cdot \mathbf{q}d\Omega }=\int\limits_{\Omega }{wQd\Omega }\] + +Using equation (\ref{eq:heat-flux}) results in, + +\[\int\limits_{\Omega }{w\rho C\frac{\partial u}{\partial t}d\Omega }+\int\limits_{\Omega }{k\nabla w\cdot \nabla ud\Omega }=\int\limits_{\Omega }{wQd\Omega }-\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }\] + +For steady-state analysis, we simply omit the first term. + +\[\int\limits_{\Omega }{k\nabla w\cdot \nabla ud\Omega }=\int\limits_{\Omega }{wQd\Omega }-\int\limits_{{{\Gamma }_{N}}}{w{{q}_{0}}d\Gamma }\quad\text{(steady state analysis)}\] + +## Discretization +The domain is discretized using iso-parametric elements. + +\[\begin{align} + & \mathbf{r}=\sum\limits_{a}{{{N}_{a}}{{\mathbf{r}}_{a}}} \\ + & u=\sum\limits_{a}{{{N}_{a}}{{u}_{a}}} \\ + & w=\sum\limits_{a}{{{N}_{a}}{{w}_{a}}} \\ +\end{align}\] + +Considering first the simpler case of steady state, the discretized equation becomes the following. + +\[\sum\limits_{a,b}{{{w}_{a}}\left( \int\limits_{\Omega }{k\nabla {{N}_{a}}\cdot \nabla {{N}_{b}}d\Omega } \right){{u}_{b}}}=\sum\limits_{a}{{{w}_{a}}\left( \int\limits_{\Omega }{{{N}_{a}}Qd\Omega } \right)}-\sum\limits_{a}{{{w}_{a}}\left( \int\limits_{{{\Gamma }_{N}}}{N_{a}^{'}{{q}_{0}}d\Gamma } \right)}\] + +Note that in the second term on the right-hand side, the shape functions are restricted to the boundary surface. + +This can be written more concisely. + +\[\mathbf{Ku}=\mathbf{F}\] + +Here, $\mathbf{K}$ is the stiffness matrix and is calculated from, + +\[{{K}_{ab}}=\int\limits_{\Omega }{k\nabla {{N}_{a}}\cdot \nabla {{N}_{b}}d\Omega }\] + +Note that the stiffness matrix is symmetric. + +The vector $\mathbf{u}$ collects all the unknown temperature values, and $\mathbf{F}$ collects all the known force loads, + +\[{{F}_{a}}=\int\limits_{\Omega }{{{N}_{a}}Qd\Omega }-\int\limits_{{{\Gamma }_{N}}}{N_{a}^{'}{{q}_{0}}d\Gamma }\] + +Note again the slight abuse of notation, since the second term is evaluated over a surface, using different shape functions than the first term. + +Returning now to the transient case, we still need to discretize the term, + +\[\int\limits_{\Omega }{w\rho C\dot{u}d\Omega }\] + +Here, we use $\dot{u}$ to denote the time derivative of the temperature. Introducing shape functions, similar as before, + +\[\dot{u}=\sum\limits_{a}{{{N}_{a}}{{{\dot{u}}}_{a}}}\] + +The discretized integral then becomes, + +\[\sum\limits_{a,b}{{{w}_{a}}\left( \int\limits_{\Omega }{\rho C{{N}_{a}}{{N}_{b}}d\Omega } \right)\,{{{\dot{u}}}_{b}}}=\sum\limits_{a,b}{{{w}_{a}}{{M}_{ab}}{{{\dot{u}}}_{b}}}\] + +where, + +\[{{M}_{ab}}=\int\limits_{\Omega }{\rho C{{N}_{a}}{{N}_{b}}d\Omega }\,\] + +is the _capacitance matrix_. + +Bringing all the terms together results in the following equation, + +\[\mathbf{M\dot{u}}+\mathbf{Ku}=\mathbf{F}\] + +This is the _semi-discrete equation_. Solving it will be discussed below. + + +## Time Discretization + +We now turn our attention to solving the semi-discrete equation, which is repeated here. + +\[\mathbf{M\dot{u}}+\mathbf{Ku}=\mathbf{F}\] + +We’ll solve this equation using the implicit-Euler method (also known as backward Euler method). This method is unconditionally stable and first-order accurate. + +First, we replace $\dot{\mathbf{u}}$ with a finite difference approximation. + +\[\mathbf{\dot{u}}=\frac{1}{\Delta t}\left( {{\mathbf{u}}^{n+1}}-{{\mathbf{u}}^{n}} \right)\] + +Then, we solve the semi-discrete equation at time point n+1. + +\[\left( \frac{1}{\Delta t}\mathbf{M}+\mathbf{K} \right){{\mathbf{u}}^{n+1}}=\mathbf{F}+\frac{1}{\Delta t}\mathbf{M}{{\mathbf{u}}^{n}}\] + +Once again, we recover a linear system of equations that can be solved. \ No newline at end of file diff --git a/docs/febioheat_tutorial1.md b/docs/febioheat_tutorial1.md new file mode 100644 index 0000000..16cda84 --- /dev/null +++ b/docs/febioheat_tutorial1.md @@ -0,0 +1,3 @@ +# Tutorial 1 + +(coming soon to a manual nearby!) diff --git a/docs/febioheat_user_manual.md b/docs/febioheat_user_manual.md new file mode 100644 index 0000000..d76fc05 --- /dev/null +++ b/docs/febioheat_user_manual.md @@ -0,0 +1,128 @@ +## Introduction + +This user manual discusses the feature of the FEBioHeat plugin in the context of the FEBio input file. + +This manual assumes the 4.0 FEBio format specification. + +## Heat Module + +The module type must be set to `heat` for solving a heat-transfer problem with this plugin. + +```xml + +``` + +## Heat transfer constitutive model + +Currently, only one constitutive model is supported for heat-transfer analysis, namely the `isotropic Fourier`, which defines the heat flux $\mathbf{q}$ as follows, + +\[\mathbf{q}=-k\,\nabla T\] + +Here, $T$ is the temperature and $k$ the _thermal conductivity_. + +This material has the following parameters as shown in the following table. + +| parameter | description |units (SI)| +|-----------|---------------------------|----------| +|`density` | The material density | kg/m^3 | +|`k` | The thermal conductivity | W/m.K | +|`c` | The specific heat | J/kg.K | + +The units are just given as an example of a consistent set of units for these variables. Users can use their preferred unit system instead. + +Note that density and specific heat are only used for transient analysis. + +An example of material definition follows below. + +```xml + + 1.0 + 0.4 + 1.0 + +``` + +## Heat Transfer Boundary Conditions +### Fixed and prescribed temperature + +To prescribe a temperature on a boundary, the “prescribed temperature” boundary condition can be used. A load controller can be assigned to the value parameter to make the prescribed temperature a function of time. For example, + +```xml + + 1.0 + 0 + +``` + +For prescribing a zero temperature on a boundary, a slightly more efficient boundary condition can be used instead, namely the `zero temperature` boundary condition. + +```xml + +``` + +This is more efficient since the degree of freedom for nodes assigned to a `zero temperature` boundary condition, is effectively ignored, resulting in smaller linear system of equations that needs to be solved. + +### Heat flux +The heat flux surface load prescribes the heat flux on a surface of the mesh. It is defined as a surface_load with the type attribute set to `heatflux`. It takes a single parameter, namely `flux`, which specifies the value of the heat flux. + +```xml + + 2.5 + +``` + +### Convective heat flux +A convective heat flux applies a flux boundary condition where the flux is proportional to the difference between the surface temperature $T$ and the ambient temperature $T_a$. + +\[{{q}_{c}}={{h}_{c}}\left( T-{{T}_{a}} \right)\] + +This boundary load is defined as a `surface_load` with the type attribute set to `convective_heatflux` and requires two parameters. + +| parameter | description | units (SI)| +|-----------|------------------------------|-----------| +| `hc` | The proportionality constant | W/m2.K | +| `Ta` | The ambient temperature | K | + +An example is given below. + +```xml + + 60.0 + 25 + +``` + +## Heat source + +A heat source can be added by defining a body_load with the type set to `heat_source`. Only one parameter is required, `Q`, to define the heat source value. + +| parameter | description | units (SI)| +|-----------|--------------------------------|-----------| +| `Q` | The constant heat source value | W/m^3 | + +```xml + + 13.5 + +``` + +## Output variables + +The heat transfer plugin adds several output variables to both the plot file and the log file. + +### Plot file variables +The following table shows a list of all new variables defined by the plugin. + +| variable | description | +|-------------|---------------------------------------| +|`temperature`|The nodal temperature at the mesh nodes| +|`heat flux` |The average heat flux over each element| + +### Log file variables + +The following table shows a list a log file variables. + +|variable | description |class | +|---------|-------------------|---------| +| `T` | nodal temperature |node_data| + diff --git a/docs/figs/FEBioHeat.png b/docs/figs/FEBioHeat.png new file mode 100644 index 0000000..be6b999 Binary files /dev/null and b/docs/figs/FEBioHeat.png differ diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..fd303f8 --- /dev/null +++ b/docs/index.md @@ -0,0 +1,14 @@ +--- +date: + created: 2026-1-18 +--- + +# FEBioHeat + +This is the documentation for **FEBioHeat**, an FEBio plugin that allows users to solve steady-state and transient heat transfer problems. + +Please see the [Tutorial](febioheat_tutorial1.md) page for a practical guide on using the FEBioHeat plugin with FEBioStudio. + +The [User Manual](febioheat_user_manual.md) page describes the features of the plugin in the context of the FEBio Input file. + +For some background on the theory, please see the [Theory](febioheat_theory.md) page. diff --git a/docs/js/mathjax_config.js b/docs/js/mathjax_config.js new file mode 100644 index 0000000..7309da9 --- /dev/null +++ b/docs/js/mathjax_config.js @@ -0,0 +1,7 @@ +window.MathJax = { + tex: { + tags: 'ams', + inlineMath: [['$', '$'], ['\\(', '\\)']], + displayMath: [['$$', '$$'], ['\\[', '\\]']] + } +}; diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 0000000..166f097 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,51 @@ +site_name: FEBioHeat +site_description: Documentation site for FEBioHeat +site_author: Steve Maas +theme: + name: material + logo: figs/FEBioHeat.png + palette: + primary: indigo + accent: indigo + font: + text: 'Roboto' + code: 'Roboto Mono' + features: + - navigation.tabs + - navigation.top + - search.highlight + - search.suggest + - toc.integrate + - content.external.links + - content.code.annotate + - content.code.copy +markdown_extensions: + - admonition + - codehilite + - footnotes + - attr_list + - md_in_html + - toc: + permalink: true + - pymdownx.arithmatex: + generic: true + - pymdownx.details + - pymdownx.superfences + - pymdownx.blocks.caption + - def_list + - pymdownx.tasklist: + custom_checkbox: true + - pymdownx.emoji: + emoji_index: !!python/name:material.extensions.emoji.twemoji + emoji_generator: !!python/name:material.extensions.emoji.to_svg + +extra_javascript: + - js/mathjax_config.js + - https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js + +nav: + - Home: index.md + - Tutorials: + - Tutorial 1: febioheat_tutorial1.md + - User Manual: febioheat_user_manual.md + - Theory: febioheat_theory.md \ No newline at end of file