# Modeling with Non-constant Coefficients in FEATool Multiphysics

**FEATool Multiphysics** has been designed to be as easy to use as
possible, and conveniently features built-in
parsing and evaluation of non-constant variables and coefficients.
This allows users to quickly and easily enter modeling expressions
formulated much as one would write on paper, all without having to
develop, write, and compile external user defined functions as other
more complex simulation software often requires.

The following applies to all physics coefficients such as point
constraints, boundary conditions, equation coefficients, initial
conditions and postprocessing expressions, as well as command line
functions such as `evalexpr`

,
`intsubd`

,
`intbdr`

,
`minmaxsubd`

,
`minmaxbdr`

, and
`postplot`

. A coefficient may be
composed as a string expression including any number of the various
coefficient types listed below, combined with the operators
`+-*/^`

, parentheses pairs `()`

, as well as
numeric constants and scalars.

## Space coordinates

The space coordinates `x`

, `y`

, and
`z`

for Cartesian coordinate systems, and `r`

and `z`

for axi-symmetric and cylindrical coordinate
systems are valid coefficient expressions. (Note that the space
coordinate names can be re-assigned from their default definitions
when starting a *New Model*).

The y-coordinate is for example used to define a parabolic velocity
inlet profile as *u _{inlet}* =

`4*0.3/0.41^2*y*(0.41-y)`

in the
flow around a cylinder CFD tutorial example.## Dependent variables

Solution unknowns or dependent variables can also be used in
coefficients, with their corresponding names given by the included
physics modes, such as for example `T`

for temperature. (As
with the space coordinate names, dependent variable names may also be
re-defined from the default names when a new physics mode is
added. Also note that adding dependent variables in coefficients often
make a problem more non-linear and harder to solve.)

The introductory
heat transfer multiphysics model for
example shows how the temperature `T`

is coupled to the
fluid flow via the source term
`alpha*g*rho*(T-Tc)`

, where *alpha*, *g*, *rho*,
and *Tc*, are model constants, while the flow velocities
`u`

and `v`

in turn are coupled back and driving
the temperature field through the convective terms.

## Space derivatives of dependent variables

By appending a space coordinate postfix to a dependent variable, for
example `Tx`

, the derivatives in the corresponding
direction can be evaluated (in the example here *Tx* is equivalent to
$\frac{\partial T}{\partial x}$). The capacitance in the
spherical capacitor model tutorial
is for example calculated by integrating the expression
`eps0*epsr*(Vr^2+Vz^2)*pi*r`

, where *V* is
the computed potential.

Second order derivatives, as for example `Txy`

equivalent
to $\frac{\partial T}{\partial x\partial y}$, can also be evaluated
for dependent variables discretized by the second order Lagrange
finite element basis functions.

Note that derivatives of dependent variables are evaluated by direct differentiation of the underlying FEM shape function polynomials and some overall accuracy may be lost in this process. For example, a variable discretized with a first order linear basis function, will feature a piecewise constant first derivative on each element, and zero second derivative (even if a global second derivative would exist analytically). For higher order accuracy with derivatives a projection or gradient reconstruction technique is necessary.

## Time

The variable name `t`

is reserved for the current time in
instationary simulations (evaluates to *0* for stationary
problems). For example, the
instationary flow around a cylinder example
defines a parabolic and time dependent inlet velocity as
*u _{inlet}* =

`6*sin(pi*t/8)*(y*(0.41-y))/0.41^2`

.## Normals

Boundary coefficients may use the external normals by prefixing a
space dimension name with the character *n*, as for example
`nx`

for the normal in the x-direction. Boundary normals
are evaluated and computed as the outward pointing unit vector from
the center of each external cell edge or face.

## Discontinuous functions

Logical switch expressions such as `T>0`

and
`T<=T0`

are also valid syntax and will evaluate to *1* when
true and *0* everywhere else. Valid expressions may include the
characters `<`

, `>`

, `<=`

,
`>=`

, `&`

, and `|`

(corresponding to
the *less than*, *greater than*, *less than or equal*, *greater than
or equal*, *and*, and *or* operators).

These types of switch expressions can be used to build quite complex
expressions and coefficients, for example in the
flow over a backwards facing step model
the recirculation region is highlighted and visualized with the
postprocessing expression
`x/hstep*(u<0)*(y<0)`

. This expression
effectively isolates the region in *y < 0* where the x-velocity *u* is
negative, and overlays it with the scaled coordinate *x/hstep*,
resulting in the following plot

Abruptly varying and discontinuous functions like these can sometimes introduce numerical instabilities if used in equation and boundary conditions. From the perspective of the solver it might then be preferable to use a regularized or smoothed Heaviside function, such as

`(abs(w)<1)*(1-(0.5*(1+w+1/pi*sin(pi*w)))) + (w<-1)`

or a smoothed Dirac delta function

`(w>-1 & w<1)*15/16*(1-2*w^2+w^4)/h_grid`

where *w* is a weighting function defining the transition and width of
the smoothed region. For example, with *w =*
`(sqrt(x^2+y^2)-0.5)/(2*h_grid)`

a circle with radius *0.5*
and transition region width *2*h_grid* is defined as shown below.

Smoothed step and delta functions are commonly used in models with moving and immersed boundaries which cannot be resolved by the mesh. Step functions are used to define the discontinuous materials, and delta functions to define surface tension forces and reactions. For example see the multiphase flow models ex_multiphase1.m.

## Built-in functions

All of the common built-in
mathematical functions and constants in MATLAB,
such as `pi`

, `eps`

, `sqrt`

,
`sin`

, `cos`

, `log`

,
`exp`

, `abs`

etc., may also be used to compose
expressions.

## Custom and user-defined functions

If a combination of the above methods is insufficient to implement an expression, or the expression can not be described analytically one can use a custom MATLAB m-file function. This could be the case if one wants to use tabulated or experimental data as input which must be interpolated, or complex coefficients with memory effects such as hysteresis. For example, if a coefficient is described as

`myfun(x,y,t,T)`

FEATool will first evaluate the inner function arguments (in this case
*x*, *y*, *t*, and *T*) and then call the function *myfun* with these
arguments. The input arguments can be any combination of the valid
coefficients described above. This assumes a function file *myfun.m*
can be found on any of the MATLAB search paths and directories. The
function must return a numeric array with the same size and dimensions
as the input arguments. A simplified function m-file could for example
look like the following

```
function [ result ] = myfun( x, y, t, T )
persistent mydata T0
if( isempty(mydata) )
load( 'path_to_mydata_file', 'data_x', 'data_y', 'data_values' );
T0 = 300;
end
interpolated_data = interp2( data_x, data_y, data_values, x, y );
result = T0 + (1-cos(2*pi*t))*T.^2.*interpolated_data;
index_limit = find( result < 0 );
result(index_limit) = 0;
```

In a custom function one can then use the usual MATLAB script
functions and techniques such as *load* (to load external data),
*save*, *interp*, *interp2*, *persistent*, *varargin*, *system*
etc. (To see information and syntax about a MATLAB command one can
enter `help command`

in the command line interface.)

Note that depending on the grid size the function may be called several times per evaluation step for different grid cell blocks (for efficiency reasons the default setting is to assembly block sizes of 50000 grid cells or less per function call).

## Reserved coefficient names

The coefficient name `h_grid`

is reserved for the mean grid
cell size or diameter, and is primarily used in the artificial and
convective stabilization techniques.

Internal and local variables may be prefixed with double underscores
`__varname`

, and it is therefore not recommended to name
variables with underscores to prevent potential name collisions.

## Tips

An easy way to test and check an expression for correctness is to simply plot and visualize it in postprocessing mode.

For two-dimensional surface plots in the FEATool GUI one can simply click on a point in the plot to evaluate the visualized expression in the point.

When working in the FEATool GUI it is often convenient to use the Model Constants and Expressions dialog box to define and store modeling expressions. The named expressions are then available everywhere in the GUI from equation coefficients to postprocessing expressions.

The more non-linearities that are introduced the harder it will be for the solver to converge to a solution. Highly non-linear problems typically require more relaxation (which can be achieved by decreasing the

*Non-linear relaxation parameter*in the Solver Settings) and a very good initial guess. One can try to introduce a scalar weighting coefficient to the non-linear coefficients, and step by step increasing it from*0*(no contributions) to*1*(fully non-linear contributions) while using the previous solution as initial guess. Lastly, for stationary problems one can also employ pseudo time stepping with the Backward-Euler scheme to try to use the time dependent solver to reach a steady state.