Page 132 - Process Modelling and Simulation With Finite Element Methods

Multiphysics 119

Now for the FEMLAB implementation. Before you launch FEMLAB from
MATLAB, make sure that you change the current directory to the one with your
watrdens.m m-file. To make sure that it is available, try out some value between
273K and 373K:

>> watrdens (330)
ans =
-0,0151

Now when you launch FEMLAB, it will inherit this current directory and have
the waterdens m-file function at its disposal. Load the saved version of
freeconv.mat from your distribution, and then Pull down the options menu and
select Add/Edit constants. Replace the Rayleigh number entry with the gravity
group Cr, and set it initially to Cr=O. Now edit the NS subdomain settings
and set
Fx=-Gr*watrdens (T)
Now onto solving. Click on the solver button (=) on the toolbar. If FEMLAB
hasn’t already popped up the message “No differentiation rule for function
watrdens”, it will now. By default, FEMLAB computes symbolic derivatives of
just about everything in sight in assembling the stiffness matrices, constraint
matrices, and load vectors. So it naturally is annoyed at us for not telling it how
to differentiate watrdens(T). FEMLAB has a place in its FEM structure for
differentiation rules if you provide a function that can be differentiated
analytically (femmles) which is used by the femdiff FEMLAB function. There
is a handle on fem.rules (Options menu, Differentiation Rules) in the
FEMLAB GUI, but we will use the MATLAB programming language.
I tried the following: Pull down the Solver menu, select Parameters, and
uncheck the tick box F in the automatic differentiation section. Now click on the
solve button. The solution progress window should now manifest as the nonlinear
solver whirls away. Eventually it reports many NaNs and Infs in the solution,
which should be interpreted as utter failure. So Plan B was necessary. I created a
second MATLAB m-file function for the numerical derivative of watrdens(T):

function a=dwatrden!ttemp)
%DWATRDENS Interpolates the expansivity of water between 273 and
373 deg K
temp=[O 3.98 5 10 15 18 20 25 30 35 38 40 45 50 55 60 65 70 75 80
85 90 95 1001;
dens=[0.99987 1. 0.99999 0.99973 0.99913 0.99862 0.99823 0,99707
0.99567 0.99406 0,99299 0.99224 0.99025 ...
0.98807 0.98573 0.98324 0.98059 0.97781 0.97489 0.97183 0.96865
0.96534 0.96192 0.958381;
temp=temp+273;
dens= (dens-dens (1) )/dens (1) ;
pp=spline(temp,dens); % Piecewise polynomial form of a cubic spline
[br,co,npy,ncol=unmkpp(pp); Breaks apart the pp form
%