Решение нелинейной FEM в MATLAB

Question

Я пытаюсь решить проблему тепловой диффузии на тетраэдрических конечных элементах с узловыми источниками тепла, которые зависят от вектора решения, в MATLAB. Нелинейная система уравнений выглядит следующим образом:

В U»+ А U = Q (T)

с Б будучи capactiy матрицы тепла, А является матрицей проводимости, д быть терминами источника и U будучи температурой. Я использую схему предсказателя-корректора Adams-Bashforth/Trapezoid Rule с итерацией Пикара, за которой следует контроль по времени. Температура для исходных терминов оценивается точно между температурой последнего временного шага и температурой предсказателя. Вот упрощенная версия кода-предиктор-корректор. Вычисление источников - это функция.

% predictor 
    K0 = t(n)-t(n-1);    
    Upre(dirichlet) = u_d_t(coordinates(dirichlet,:)); 
    Upre(FreeNodes) = U(FreeNodes,n) + (dt/2)*((2+dt/K0)*U_dt(FreeNodes,3) - (dt/K0)*U_dt(FreeNodes,1));  % predictor step 
    Uguess = Upre;  % save as initial guess for Picard iteration 

    % corrector with picard iteration 
    while res >= picard_tolerance 
     T_theta = Uguess*theta + (1-theta)*U(:,n); 
     b = q(T_theta); 
     % Building right-hand side vector (without Dirichlet boundary conditions yet) 
     rhs = ((2/dt)*B*U(:,n) + B*U_dt(:,1))+b; 

     % Applying Dirichlet Boundary Conditions to the Solution Vector 
     Ucor(dirichlet) = u_d_t(coordinates(dirichlet,:)); 
     rhs = rhs - ((2/dt)*B+A)*Ucor; 

     % Solving the linearized system using the backslash operator 
     % P*U(n+1) = f(Un) => U(n+1) = P\f(Un) 
     Ucor(FreeNodes) = ((2/dt)*B(FreeNodes,FreeNodes)+A(FreeNodes,FreeNodes))\rhs(FreeNodes); 
     res = norm(Uguess-Ucor); 
     Uguess = Ucor; 
     U(:,n+1) = Ucor; 
    end 

Как вы можете видеть, я использую оператор обратной косой черты для решения системы. Нелинейность системы не должна быть плохой. Однако с увеличением временных шагов метод пикарда сходится медленнее и, в конечном итоге, перестает сходиться. Мне нужно гораздо больше времени, хотя, поэтому я поставил весь корректирующий шаг в функцию и попытался решить его с помощью fsolve, чтобы убедиться, что я добился более быстрой конвергенции. К сожалению, fsolve, похоже, никогда не заканчивает первый шаг. Я полагаю, что я не настроил параметры fsolve правильно. Может ли кто-нибудь сказать мне, как настроить fsolve для больших разреженных нелинейных систем (мы говорим о тысячах до сотен тысяч уравнений). Или может быть лучшее решение, чем fsolve для этой проблемы? Помощь и - поскольку я не специалист или инженер-вычислитель, явный совет очень ценен!

Answer 1

В моем опыте нелинейные уравнения решаются путем линеаризации для решения приращения температуры и итерации к конвергенции с использованием чего-то вроде решения Ньютона Рафсона. Поэтому, если вы используете неявную схему интеграции, у вас есть внешнее временное решение с внутренним нелинейным решением для температурного шага по временному шагу.